diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py index 8405650d..0272d702 100644 --- a/nlmod/dims/grid.py +++ b/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 diff --git a/nlmod/dims/layers.py b/nlmod/dims/layers.py index 9da3d8c7..89d667f7 100644 --- a/nlmod/dims/layers.py +++ b/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): diff --git a/tests/test_006_caching.py b/tests/test_006_caching.py index 214186b3..d82f68c3 100644 --- a/tests/test_006_caching.py +++ b/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.""" diff --git a/tests/test_009_layers.py b/tests/test_009_layers.py index e93e2a8a..00d37a39 100644 --- a/tests/test_009_layers.py +++ b/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) diff --git a/tests/test_013_surface_water.py b/tests/test_013_surface_water.py index cfc77507..519fc933 100644 --- a/tests/test_013_surface_water.py +++ b/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,6 +56,10 @@ 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) @@ -62,7 +67,7 @@ def test_gdf_lake(): 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( @@ -70,14 +75,20 @@ def test_gdf_lake(): "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") diff --git a/tests/test_026_grid.py b/tests/test_026_grid.py index b5dc4d4d..c27a3c03 100644 --- a/tests/test_026_grid.py +++ b/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()