Add tests for layer methods

nlmod/dims/grid.py

@@ -1316,7 +1316,7 @@ def gdf_to_da(
     agg_method : str, optional
         aggregation method to handle multiple geometries in one cell, options
         are:
-        - max, min, mean,
+        - max, min, mean, sum
         - length_weighted (lines), max_length (lines),
         - area_weighted (polygon), max_area (polygon).
         The default is 'max'.
@@ -1353,13 +1353,12 @@ def gdf_to_da(
             raise ValueError(
                 "multiple geometries in one cell please define aggregation method"
             )
+        modelgrid = None
         if agg_method in ["nearest"]:
             modelgrid = modelgrid_from_ds(ds)
-            gdf_agg = aggregate_vector_per_cell(
-                gdf_cellid, {column: agg_method}, modelgrid
-            )
-        else:
-            gdf_agg = aggregate_vector_per_cell(gdf_cellid, {column: agg_method})
+        gdf_agg = aggregate_vector_per_cell(
+            gdf_cellid, {column: agg_method}, modelgrid=modelgrid
+        )
     else:
         # aggregation not neccesary
         gdf_agg = gdf_cellid[[column]]
@@ -1512,18 +1511,18 @@ def aggregate_vector_per_cell(gdf, fields_methods, modelgrid=None):
     # check geometry types
     geom_types = gdf.geometry.type.unique()
     if len(geom_types) > 1:
-        if (
-            len(geom_types) == 2
-            and ("Polygon" in geom_types)
-            and ("MultiPolygon" in geom_types)
+        if len(geom_types) == 2 and (
+            set(geom_types) == set(["LineString", "MultiLineString"])
+            or set(geom_types) == set(["Polygon", "MultiPolygon"])
         ):
             pass
         else:
             raise TypeError("cannot aggregate geometries of different types")
     if bool({"length_weighted", "max_length"} & set(fields_methods.values())):
-        assert (
-            geom_types[0] == "LineString"
-        ), "can only use length methods with line geometries"
+        if ("LineString" in geom_types) or ("MultiLineString" in geom_types):
+            pass
+        else:
+            raise TypeError("can only use length methods with line geometries")
     if bool({"area_weighted", "max_area"} & set(fields_methods.values())):
         if ("Polygon" in geom_types) or ("MultiPolygon" in geom_types):
             pass

nlmod/dims/layers.py

@@ -1470,6 +1470,8 @@ def get_first_active_layer_from_idomain(idomain, nodata=-999):
 
     first_active_layer = xr.where(idomain[0] == 1, 0, nodata)
     for i in range(1, idomain.shape[0]):
+        if not (first_active_layer == nodata).any():
+            break
         first_active_layer = xr.where(
             (first_active_layer == nodata) & (idomain[i] == 1),
             i,
@@ -1479,6 +1481,26 @@ def get_first_active_layer_from_idomain(idomain, nodata=-999):
     return first_active_layer
 
 
+def get_last_active_layer(ds, **kwargs):
+    """Get the last active layer in each cell from a model ds.
+
+    Parameters
+    ----------
+    ds : xr.DataSet
+        Model Dataset
+    **kwargs : dict
+        Kwargs are passed on to get_last_active_layer_from_idomain.
+
+    Returns
+    -------
+    last_active_layer : xr.DataArray
+        raster in which each cell has the zero based number of the last
+        active layer. Shape can be (y, x) or (icell2d)
+    """
+    idomain = get_idomain(ds)
+    return get_last_active_layer_from_idomain(idomain, **kwargs)
+
+
 def get_last_active_layer_from_idomain(idomain, nodata=-999):
     """Get the last (bottom) active layer in each cell from the idomain.
 
@@ -1498,8 +1520,10 @@ def get_last_active_layer_from_idomain(idomain, nodata=-999):
     """
     logger.debug("get last active modellayer for each cell in idomain")
 
-    last_active_layer = xr.where(idomain[-1] == 1, 0, nodata)
+    last_active_layer = xr.where(idomain[-1] == 1, idomain.shape[0] - 1, nodata)
     for i in range(idomain.shape[0] - 2, -1, -1):
+        if not (last_active_layer == nodata).any():
+            break
         last_active_layer = xr.where(
             (last_active_layer == nodata) & (idomain[i] == 1),
             i,
@@ -1690,7 +1714,7 @@ def check_elevations_consistency(ds):
     thickness = calculate_thickness(ds)
     mask = thickness < 0.0
     if mask.any():
-        logger.warning(f"Thickness of layers is negative in {mask.sum()} cells.")
+        logger.warning(f"Thickness of layers is negative in {int(mask.sum())} cells.")
 
 
 def insert_layer(ds, name, top, bot, kh=None, kv=None, copy=True):

tests/test_006_caching.py

@@ -45,6 +45,9 @@ def test_cache_ahn_data_array():
             modification_time1 != modification_time3
         ), "Cache should have been rewritten"
 
+        # clear cache again
+        nlmod.cache.clear_cache(tmpdir, prompt=False)
+
 
 def test_cache_northsea_data_array():
     """Test caching of AHN data array. Does have dataset as argument."""

tests/test_009_layers.py

@@ -9,14 +9,13 @@
 from nlmod.plot import DatasetCrossSection
 
 
-def get_regis_horstermeer():
+def get_regis_horstermeer(cachedir=None, cachename="regis_horstermeer"):
     extent = [131000, 136800, 471500, 475700]
-    cachedir = os.path.join(os.path.dirname(os.path.realpath(__file__)), "data")
+    if cachedir is None:
+        cachedir = os.path.join(os.path.dirname(os.path.realpath(__file__)), "data")
     if not os.path.isdir(cachedir):
         os.makedirs(cachedir)
-    regis = nlmod.read.get_regis(
-        extent, cachedir=cachedir, cachename="regis_horstermeer"
-    )
+    regis = nlmod.read.get_regis(extent, cachedir=cachedir, cachename=cachename)
     return regis
 
 
@@ -159,6 +158,68 @@ def test_set_minimum_layer_thickness(plot=False):
         plot_test(regis, ds_new)
 
 
+def test_calculate_transmissivity():
+    regis = get_regis_horstermeer()
+    nlmod.layers.calculate_transmissivity(regis)
+
+
+def test_calculate_resistance():
+    regis = get_regis_horstermeer()
+    # with the default value of between_layers=True
+    nlmod.layers.calculate_resistance(regis)
+    # and also with between_layers=False
+    nlmod.layers.calculate_resistance(regis, between_layers=False)
+
+
+def test_get_layer_of_z():
+    regis = get_regis_horstermeer()
+    z = -100
+
+    layer = nlmod.layers.get_layer_of_z(regis, z)
+
+    assert (regis["botm"].isel(layer=layer) < z).all()
+    top = regis["botm"] + nlmod.layers.calculate_thickness(regis)
+    assert (top.isel(layer=layer) > z).all()
+
+
+def test_aggregate_by_weighted_mean_to_ds():
+    regis = get_regis_horstermeer()
+    regis2 = regis.copy(deep=True)
+
+    # botm needs to have the name "bottom'
+    regis2["bottom"] = regis2["botm"]
+    # top needs to be 3d
+    regis2["top"] = regis2["botm"] + nlmod.layers.calculate_thickness(regis2)
+    kh_new = nlmod.layers.aggregate_by_weighted_mean_to_ds(regis, regis2, "kh")
+    assert np.abs(kh_new - regis["kh"]).max() < 1e-5
+    # assert (kh_new.isnull() == regis["kh"].isnull()).all() # does not assert to True...
+
+
+def test_check_elevations_consistency(caplog):
+    regis = get_regis_horstermeer()
+    # there are no inconsistencies in this dataset, let's check for that:
+    nlmod.layers.check_elevations_consistency(regis)
+    assert len(caplog.text) == 0
+
+    # add an inconsistency by lowering the top of the model in part of the model domain
+    regis["top"][10:20, 20:25] = -5
+    nlmod.layers.check_elevations_consistency(regis)
+    assert "check_elevations_consistency" not in caplog.text
+    assert len(caplog.text) > 0
+    assert "Thickness of layers is negative in 50 cells" in caplog.text
+
+
+def test_get_first_and_last_active_layer():
+    regis = get_regis_horstermeer()
+    thickness = nlmod.layers.calculate_thickness(regis)
+
+    fal = nlmod.layers.get_first_active_layer(regis)
+    assert (thickness[fal] > 0).all()
+
+    lal = nlmod.layers.get_last_active_layer(regis)
+    assert (thickness[lal] > 0).all()
+
+
 def test_set_model_top(plot=False):
     regis = get_regis_horstermeer()
     ds_new = nlmod.layers.set_model_top(regis.copy(deep=True), 5.0)

tests/test_013_surface_water.py

@@ -1,5 +1,6 @@
 import os
 
+import numpy as np
 import geopandas as gpd
 import pandas as pd
 import flopy
@@ -55,29 +56,39 @@ def test_gdf_lake():
     )
     ds = nlmod.time.set_ds_time(ds, time=[1], start=pd.Timestamp.today())
     ds = nlmod.dims.refine(ds)
+    dims = ("time", "icell2d")
+    size = (len(ds.time), len(ds.icell2d))
+    ds["recharge"] = dims, np.full(size, 0.002)
+    ds["evaporation"] = dims, np.full(size, 0.001)
 
     sim = nlmod.sim.sim(ds)
     nlmod.sim.tdis(ds, sim)
     nlmod.sim.ims(sim)
     gwf = nlmod.gwf.gwf(ds, sim)
     nlmod.gwf.dis(ds, gwf)
 
-    ds["evap"] = (("time",), [0.0004])
+    ds["lake_evap"] = (("time",), [0.0004])
 
     # add lake with outlet and evaporation
     gdf_lake = gpd.GeoDataFrame(
         {
             "name": ["lake_0", "lake_0", "lake_1"],
             "strt": [1.0, 1.0, 2.0],
             "clake": [10.0, 10.0, 10.0],
-            "EVAPORATION": ["evap", "evap", "evap"],
+            "EVAPORATION": ["lake_evap", "lake_evap", "lake_evap"],
             "lakeout": ["lake_1", "lake_1", None],
             "outlet_invert": ["use_elevation", "use_elevation", None],
         },
         index=[14, 15, 16],
     )
 
-    nlmod.gwf.lake_from_gdf(gwf, gdf_lake, ds, boundname_column="name")
+    rainfall, evaporation = nlmod.gwf.clip_meteorological_data_from_ds(
+        gdf_lake, ds, boundname_column="name"
+    )
+    # do not pass evaporation to lake_from_gdf, as we have specified it in gdf_lake
+    nlmod.gwf.lake_from_gdf(
+        gwf, gdf_lake, ds, boundname_column="name", rainfall=rainfall
+    )
 
     # remove lake package
     gwf.remove_package("LAK_0")

tests/test_026_grid.py

@@ -134,6 +134,48 @@ def test_fillnan_da():
     assert not top[mask].isnull().any()
 
 
+def test_interpolate_gdf_to_array():
+    bgt = get_bgt()
+    bgt.geometry = bgt.centroid
+    bgt["values"] = range(len(bgt))
+
+    regis = get_regis()
+    ds = nlmod.to_model_ds(regis, model_ws=model_ws)
+    sim = nlmod.sim.sim(ds)
+    gwf = nlmod.gwf.gwf(ds, sim)
+    nlmod.gwf.dis(ds, gwf)
+
+    nlmod.grid.interpolate_gdf_to_array(bgt, gwf, field="values", method="linear")
+
+
+def test_gdf_to_da_methods():
+    bgt = get_bgt()
+    regis = get_regis()
+    ds = nlmod.to_model_ds(regis)
+    bgt["values"] = range(len(bgt))
+
+    bgt_line = bgt.copy()
+    bgt_line.geometry = bgt.boundary
+
+    for agg_method in [
+        "nearest",
+        "area_weighted",
+        "max_area",
+        "length_weighted",
+        "max_length",
+        # "center_grid",
+        "max",
+        "min",
+        "mean",
+        "sum",
+    ]:
+        if agg_method in ["length_weighted", "max_length"]:
+            gdf = bgt_line
+        else:
+            gdf = bgt
+        nlmod.grid.gdf_to_da(gdf, ds, column="values", agg_method=agg_method)
+
+
 def test_gdf_to_bool_da():
     bgt = get_bgt()
 
@@ -158,6 +200,30 @@ def test_gdf_to_bool_da():
     assert da.any()
 
 
+def test_gdf_to_count_da():
+    bgt = get_bgt()
+
+    # test for a structured grid
+    ds = get_structured_model_ds()
+    da = nlmod.grid.gdf_to_count_da(bgt, ds)
+    assert da.any()
+
+    # test for a vertex grid
+    ds = get_vertex_model_ds()
+    da = nlmod.grid.gdf_to_count_da(bgt, ds)
+    assert da.any()
+
+    # tets for a slightly rotated structured grid
+    ds = get_structured_model_ds_rotated()
+    da = nlmod.grid.gdf_to_count_da(bgt, ds)
+    assert da.any()
+
+    # test for a rotated vertex grid
+    ds = get_vertex_model_ds_rotated()
+    da = nlmod.grid.gdf_to_count_da(bgt, ds)
+    assert da.any()
+
+
 def test_gdf_to_da():
     bgt = get_bgt()