From 098c680024b52a1efd19d1f7c312bed48b67b87a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ruben=20Calj=C3=A9?= Date: Thu, 30 Nov 2023 12:41:02 +0100 Subject: [PATCH] Improve output (#297) * Add add_nans-argument to nlmod.gwf.output.get_heads_da * Add parameters for hdry and hnoflo And add tests * Add test for a transport model and geotop and a minor change in gis.py * Add test for nlmod.plot.map_array * Add modpath test (and add support for structured grids) * Test if docs complete without geotop notebook * Fix notebooks * Revert "Test if docs complete without geotop notebook" This reverts commit b130a1f51d461093e9f89e38f6e0f0a15cbf1bba. * Only refine Lek at the boundary * Increase speed of notebooks * Update path of model of notebook 17 (so it is in gitignore) * Change hhnk url, which failed in test * Handle comments from @dbrakenhoff Remove add_nans again, and just use the value for hdry or hnoflo to not add NaNs --- docs/examples/02_surface_water.ipynb | 14 +- docs/examples/09_schoonhoven.ipynb | 8 +- docs/examples/11_grid_rotation.ipynb | 14 +- docs/examples/14_stromingen_example.ipynb | 2 +- docs/examples/15_geotop.ipynb | 6 +- docs/examples/16_groundwater_transport.ipynb | 10 +- docs/examples/17_unsaturated_zone_flow.ipynb | 3 +- nlmod/dims/grid.py | 59 +++++--- nlmod/gis.py | 2 +- nlmod/mfoutput/mfoutput.py | 29 ++-- nlmod/modpath/modpath.py | 49 ++++--- nlmod/read/waterboard.py | 2 +- tests/test_001_model.py | 2 +- tests/test_002_regis_geotop.py | 17 +++ tests/test_012_plot.py | 2 + tests/test_015_gwf_output.py | 7 +- tests/test_022_gwt.py | 133 ++++++++++++++++++- tests/test_025_modpath.py | 45 +++++++ 18 files changed, 318 insertions(+), 86 deletions(-) create mode 100644 tests/test_025_modpath.py diff --git a/docs/examples/02_surface_water.ipynb b/docs/examples/02_surface_water.ipynb index 143d4eb7..a42b10dc 100644 --- a/docs/examples/02_surface_water.ipynb +++ b/docs/examples/02_surface_water.ipynb @@ -223,7 +223,7 @@ "metadata": {}, "source": [ "#### Save the data to use in other notebooks as well\n", - "We save the bgt-data to a geojson file, so we can use the data in other notebooks with surface water as well" + "We save the bgt-data to a GeoPackage file, so we can use the data in other notebooks with surface water as well" ] }, { @@ -233,8 +233,8 @@ "metadata": {}, "outputs": [], "source": [ - "fname_bgt = os.path.join(cachedir, \"bgt.geojson\")\n", - "bgt.to_file(fname_bgt, driver=\"GeoJSON\")" + "fname_bgt = os.path.join(cachedir, \"bgt.gpkg\")\n", + "bgt.to_file(fname_bgt)" ] }, { @@ -787,14 +787,6 @@ "cbar = fig.colorbar(qm, shrink=1.0)\n", "cbar.set_label(\"head [m NAP]\")" ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6e5833d6", - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { diff --git a/docs/examples/09_schoonhoven.ipynb b/docs/examples/09_schoonhoven.ipynb index d915b077..aceedd44 100644 --- a/docs/examples/09_schoonhoven.ipynb +++ b/docs/examples/09_schoonhoven.ipynb @@ -91,7 +91,7 @@ "metadata": {}, "outputs": [], "source": [ - "fname_bgt = os.path.join(cachedir, \"bgt.geojson\")\n", + "fname_bgt = os.path.join(cachedir, \"bgt.gpkg\")\n", "if not os.path.isfile(fname_bgt):\n", " raise (\n", " Exception(\n", @@ -207,7 +207,7 @@ "metadata": {}, "outputs": [], "source": [ - "refinement_features = [(bgt[bgt[\"bronhouder\"] == \"L0002\"], 2)]\n", + "refinement_features = [(bgt[bgt[\"bronhouder\"] == \"L0002\"].dissolve().boundary, 2)]\n", "ds = nlmod.grid.refine(ds, refinement_features=refinement_features)" ] }, @@ -503,7 +503,7 @@ " colorbar_label=\"head [m NAP]\",\n", " title=\"mean head\",\n", ")\n", - "bgt.plot(ax=pc.axes, edgecolor=\"k\", facecolor=\"none\");" + "bgt.dissolve().plot(ax=pc.axes, edgecolor=\"k\", facecolor=\"none\");" ] }, { @@ -745,7 +745,7 @@ "line = ax.add_collection(lc)\n", "nlmod.plot.colorbar_inside(line, label=\"Travel time (years)\")\n", "\n", - "bgt.plot(ax=ax, edgecolor=\"k\", facecolor=\"none\")" + "bgt.dissolve().plot(ax=ax, edgecolor=\"k\", facecolor=\"none\")" ] }, { diff --git a/docs/examples/11_grid_rotation.ipynb b/docs/examples/11_grid_rotation.ipynb index 3f4d9132..7d60bdc0 100644 --- a/docs/examples/11_grid_rotation.ipynb +++ b/docs/examples/11_grid_rotation.ipynb @@ -165,17 +165,9 @@ "metadata": {}, "outputs": [], "source": [ - "bgt = nlmod.gwf.surface_water.get_gdf(ds)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "e023dbca", - "metadata": {}, - "outputs": [], - "source": [ - "bgt.plot()" + "bgt = nlmod.read.bgt.get_bgt(extent)\n", + "bgt = nlmod.gwf.surface_water.add_stages_from_waterboards(bgt, extent=extent)\n", + "bgt = nlmod.grid.gdf_to_grid(bgt, ds).set_index(\"cellid\")" ] }, { diff --git a/docs/examples/14_stromingen_example.ipynb b/docs/examples/14_stromingen_example.ipynb index 0dea2056..10da127f 100644 --- a/docs/examples/14_stromingen_example.ipynb +++ b/docs/examples/14_stromingen_example.ipynb @@ -106,7 +106,7 @@ "# wells = nlmod.read.get_wells(extent) # no well data is available just yet\n", "\n", "# surface water features and levels\n", - "fname_bgt = os.path.join(cachedir, \"bgt.geojson\")\n", + "fname_bgt = os.path.join(cachedir, \"bgt.gpkg\")\n", "if not os.path.isfile(fname_bgt):\n", " raise (\n", " Exception(\n", diff --git a/docs/examples/15_geotop.ipynb b/docs/examples/15_geotop.ipynb index 25dc5321..13500e10 100644 --- a/docs/examples/15_geotop.ipynb +++ b/docs/examples/15_geotop.ipynb @@ -40,7 +40,7 @@ "outputs": [], "source": [ "# Define an extent and a line from southwest to northeast through this extent\n", - "extent = [95000.0, 105000.0, 494000.0, 500000.0]\n", + "extent = [100000.0, 105000.0, 499800.0, 500000.0]\n", "line = LineString([(extent[0], extent[2]), (extent[1], extent[3])])" ] }, @@ -386,13 +386,13 @@ "var = \"kh\"\n", "norm = matplotlib.colors.Normalize(0.0, 40.0)\n", "\n", - "f, axes = nlmod.plot.get_map(extent, ncols=2)\n", + "f, axes = nlmod.plot.get_map(extent, nrows=2, figsize=20)\n", "pc = nlmod.plot.data_array(regis[var].loc[layer], ax=axes[0], norm=norm)\n", "nlmod.plot.colorbar_inside(pc, bounds=[0.02, 0.05, 0.02, 0.9], ax=axes[0])\n", "nlmod.plot.title_inside(\"REGIS\", ax=axes[0])\n", "pc = nlmod.plot.data_array(regis[f\"{var}_gt\"].loc[layer], ax=axes[1], norm=norm)\n", "nlmod.plot.title_inside(\"GeoTOP\", ax=axes[1])\n", - "nlmod.plot.colorbar_inside(pc, bounds=[0.02, 0.05, 0.02, 0.9], ax=axes[1])" + "nlmod.plot.colorbar_inside(pc, bounds=[0.02, 0.05, 0.02, 0.9], ax=axes[1]);" ] }, { diff --git a/docs/examples/16_groundwater_transport.ipynb b/docs/examples/16_groundwater_transport.ipynb index 6b1c8d93..34d70dc7 100644 --- a/docs/examples/16_groundwater_transport.ipynb +++ b/docs/examples/16_groundwater_transport.ipynb @@ -67,7 +67,7 @@ "start_time = \"2010-1-1\"\n", "starting_head = 1.0\n", "\n", - "municipalities = nlmod.read.boundaries.get_municipalities(extent=extent_hbossche)" + "municipalities = nlmod.read.administrative.get_municipalities(extent=extent_hbossche)" ] }, { @@ -600,14 +600,8 @@ } ], "metadata": { - "kernelspec": { - "display_name": "artesia", - "language": "python", - "name": "python3" - }, "language_info": { - "name": "python", - "version": "3.10.12" + "name": "python" } }, "nbformat": 4, diff --git a/docs/examples/17_unsaturated_zone_flow.ipynb b/docs/examples/17_unsaturated_zone_flow.ipynb index 6dcd1256..b9500d4a 100644 --- a/docs/examples/17_unsaturated_zone_flow.ipynb +++ b/docs/examples/17_unsaturated_zone_flow.ipynb @@ -66,8 +66,7 @@ "metadata": {}, "outputs": [], "source": [ - "model_name = \"test_uzf\"\n", - "model_ws = os.path.join(\"data\", model_name)\n", + "model_name = model_ws = \"model17\"\n", "figdir, cachedir = nlmod.util.get_model_dirs(model_ws)" ] }, diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py index 42e4e3f3..13bc8715 100644 --- a/nlmod/dims/grid.py +++ b/nlmod/dims/grid.py @@ -107,12 +107,37 @@ def xy_to_icell2d(xy, ds): icell2d : int number of the icell2d value of a cell containing the xy point. """ - + msg = "xy_to_icell2d can only be applied to a vertex grid" + assert ds.gridtype == "vertex", msg icell2d = (np.abs(ds.x.data - xy[0]) + np.abs(ds.y.data - xy[1])).argmin().item() return icell2d +def xy_to_row_col(xy, ds): + """get the row and column values of a point defined by its x and y coordinates. + + Parameters + ---------- + xy : list, tuple + coordinates of ta point. + ds : xarary dataset + model dataset. + + Returns + ------- + row : int + number of the row value of a cell containing the xy point. + col : int + number of the column value of a cell containing the xy point. + """ + msg = "xy_to_row_col can only be applied to a structured grid" + assert ds.gridtype == "structured", msg + row = np.abs(ds.y.data - xy[1]).argmin() + col = np.abs(ds.x.data - xy[0]).argmin() + return row, col + + def xyz_to_cid(xyz, ds=None, modelgrid=None): """get the icell2d value of a point defined by its x and y coordinates. @@ -1142,11 +1167,10 @@ def gdf_to_data_array_struc( def gdf_to_da( - gdf, ds, column, agg_method=None, fill_value=np.NaN, min_total_overlap=0.0, - ix=None + gdf, ds, column, agg_method=None, fill_value=np.NaN, min_total_overlap=0.0, ix=None ): """Project vector data on a grid. Aggregate data if multiple - geometries are in a single cell. Supports structured and vertex grids. + geometries are in a single cell. Supports structured and vertex grids. This method replaces gdf_to_data_array_struc. Parameters @@ -1385,12 +1409,17 @@ def aggregate_vector_per_cell(gdf, fields_methods, modelgrid=None): gdf_copy["area_times_field"] = gdf_copy["area"] * gdf_copy[field] # skipna is not implemented by groupby therefore we use min_count=1 - celldata[field] = ( - gdf_copy.groupby(by="cellid")["area_times_field"].sum(min_count=1) / - gdf_copy.groupby(by="cellid")["area"].sum(min_count=1) - ) - - elif method in ("nearest", "length_weighted", "max_length", "max_area", "center_grid"): + celldata[field] = gdf_copy.groupby(by="cellid")["area_times_field"].sum( + min_count=1 + ) / gdf_copy.groupby(by="cellid")["area"].sum(min_count=1) + + elif method in ( + "nearest", + "length_weighted", + "max_length", + "max_area", + "center_grid", + ): for cid, group in tqdm(gr, desc="Aggregate vector data"): agg_dic = _get_aggregates_values(group, {field: method}, modelgrid) @@ -1403,7 +1432,7 @@ def aggregate_vector_per_cell(gdf, fields_methods, modelgrid=None): return celldata -def gdf_to_bool_da(gdf, ds, ix=None, buffer=0., **kwargs): +def gdf_to_bool_da(gdf, ds, ix=None, buffer=0.0, **kwargs): """convert a GeoDataFrame with polygon geometries into a data array corresponding to the modelgrid in which each cell is 1 (True) if one or more geometries are (partly) in that cell. @@ -1430,7 +1459,7 @@ def gdf_to_bool_da(gdf, ds, ix=None, buffer=0., **kwargs): return da -def gdf_to_bool_ds(gdf, ds, da_name, keep_coords=None, ix=None, buffer=0., **kwargs): +def gdf_to_bool_ds(gdf, ds, da_name, keep_coords=None, ix=None, buffer=0.0, **kwargs): """convert a GeoDataFrame with polygon geometries into a model dataset with a data_array named 'da_name' in which each cell is 1 (True) if one or more geometries are (partly) in that cell. @@ -1465,7 +1494,7 @@ def gdf_to_bool_ds(gdf, ds, da_name, keep_coords=None, ix=None, buffer=0., **kwa return ds_out -def gdf_to_count_da(gdf, ds, ix=None, buffer=0., **kwargs): +def gdf_to_count_da(gdf, ds, ix=None, buffer=0.0, **kwargs): """Counts in how many polygons a coordinate of ds appears. Parameters @@ -1509,7 +1538,7 @@ def gdf_to_count_da(gdf, ds, ix=None, buffer=0., **kwargs): geoms = [gdf] for geom in geoms: - if buffer > 0.: + if buffer > 0.0: cids = ix.intersects(geom.buffer(buffer), **kwargs)["cellids"] else: cids = ix.intersects(geom, **kwargs)["cellids"] @@ -1526,7 +1555,7 @@ def gdf_to_count_da(gdf, ds, ix=None, buffer=0., **kwargs): return da -def gdf_to_count_ds(gdf, ds, da_name, keep_coords=None, ix=None, buffer=0., **kwargs): +def gdf_to_count_ds(gdf, ds, da_name, keep_coords=None, ix=None, buffer=0.0, **kwargs): """Counts in how many polygons a coordinate of ds appears. Parameters diff --git a/nlmod/gis.py b/nlmod/gis.py index 90ecfb28..d9ee9ce3 100644 --- a/nlmod/gis.py +++ b/nlmod/gis.py @@ -396,7 +396,7 @@ def ds_to_ugrid_nc_file( # Convert boolean layers to integer and int64 to int32 for var in variables: if np.issubdtype(ds[var].dtype, bool): - ds[var].encoding["dtype"] = np.int + ds[var].encoding["dtype"] = np.int32 elif np.issubdtype(ds[var].dtype, str) or np.issubdtype(ds[var].dtype, object): # convert the string to an index of unique strings index = np.unique(ds[var], return_inverse=True)[1] diff --git a/nlmod/mfoutput/mfoutput.py b/nlmod/mfoutput/mfoutput.py index 8d889941..64158e3d 100644 --- a/nlmod/mfoutput/mfoutput.py +++ b/nlmod/mfoutput/mfoutput.py @@ -77,7 +77,7 @@ def _get_time_index(fobj, ds=None, gwf_or_gwt=None): return tindex -def _create_da(arr, modelgrid, times): +def _create_da(arr, modelgrid, times, hdry=-1e30, hnoflo=1e30): """Create data array based on array, modelgrid, and time array. Parameters @@ -88,6 +88,12 @@ def _create_da(arr, modelgrid, times): flopy modelgrid object times : list or array list or array containing times as floats (usually in days) + hdry : float, optional + The value of dry cells, which will be replaced by NaNs. If hdry is None, the + values of dry cells will not be replaced by NaNs. The default is -1e30. + hnoflo : float, optional + The value of no-flow cells, which will be replaced by NaNs. If hnoflo is None, + the values of no-flow cells will not be replaced by NaNs. The default is 1e30. Returns ------- @@ -99,10 +105,15 @@ def _create_da(arr, modelgrid, times): dims, coords = get_dims_coords_from_modelgrid(modelgrid) da = xr.DataArray(data=arr, dims=("time",) + dims, coords=coords) - # set dry/no-flow to nan - hdry = -1e30 - hnoflo = 1e30 - da = da.where((da != hdry) & (da != hnoflo)) + if hdry is not None or hnoflo is not None: + # set dry/no-flow to nan + if hdry is None: + mask = da != hnoflo + elif hnoflo is None: + mask = da != hdry + else: + mask = (da != hdry) & (da != hnoflo) + da = da.where(mask) # set local time coordinates da.coords["time"] = ds_time_idx(times) @@ -168,7 +179,7 @@ def _get_heads_da( stacked_arr = stacked_arr[:, :, 0, :] # create data array - da = _create_da(stacked_arr, modelgrid, hobj.get_times()) + da = _create_da(stacked_arr, modelgrid, hobj.get_times(), **kwargs) return da @@ -263,10 +274,12 @@ def _get_flopy_data_object(var, ds=None, gwml=None, fname=None, grbfile=None): extension = "_gwt.ucn" else: raise (ValueError(f"Unknown variable {var}")) - msg = f"Load the {var}s using either ds, {ml_name} or fname" - assert ((ds is not None) + (gwml is not None) + (fname is not None)) == 1, msg + if fname is None: if ds is None: + if gwml is None: + msg = f"Load the {var}s using either ds, {ml_name} or fname" + raise (ValueError(msg)) # return gwf.output.head(), gwf.output.budget() or gwt.output.concentration() return getattr(gwml.output, var)() fname = os.path.join(ds.model_ws, ds.model_name + extension) diff --git a/nlmod/modpath/modpath.py b/nlmod/modpath/modpath.py index 20db899e..125c3aa3 100644 --- a/nlmod/modpath/modpath.py +++ b/nlmod/modpath/modpath.py @@ -10,7 +10,7 @@ from packaging.version import parse as parse_version from .. import util -from ..dims.grid import xy_to_icell2d +from ..dims.grid import xy_to_icell2d, xy_to_row_col logger = logging.getLogger(__name__) @@ -82,14 +82,26 @@ def xy_to_nodes(xy_list, mpf, ds, layer=0): nodes = [] for i, xy in enumerate(xy_list): - icell2d = xy_to_icell2d(xy, ds) - if mpf.ib[layer[i], icell2d] > 0: - node = layer[i] * mpf.ib.shape[1] + icell2d - nodes.append(node) + if len(mpf.ib.shape) == 3: + row, col = xy_to_row_col(xy, ds) + if mpf.ib[layer[i], row, col] > 0: + nodes.append(get_node_structured(layer[i], row, col, mpf.ib.shape)) + else: + icell2d = xy_to_icell2d(xy, ds) + if mpf.ib[layer[i], icell2d] > 0: + nodes.append(get_node_vertex(layer[i], icell2d, mpf.ib.shape)) return nodes +def get_node_structured(lay, row, col, shape): + return lay * shape[1] * shape[2] + row * shape[2] + col + + +def get_node_vertex(lay, icell2d, shape): + return lay * shape[1] + icell2d + + def package_to_nodes(gwf, package_name, mpf): """Return a list of nodes from the cells with certain boundary conditions. @@ -123,10 +135,12 @@ def package_to_nodes(gwf, package_name, mpf): ) nodes = [] for cid in pkg_cid: - if mpf.ib[cid[0], cid[1]] > 0: - node = cid[0] * mpf.ib.shape[1] + cid[1] - nodes.append(node) - + if len(mpf.ib.shape) == 3: + if mpf.ib[cid[0], cid[1], cid[2]] > 0: + nodes.append(get_node_structured(cid[0], cid[1], cid[2], mpf.ib.shape)) + else: + if mpf.ib[cid[0], cid[1]] > 0: + nodes.append(get_node_vertex(cid[0], cid[1], mpf.ib.shape)) return nodes @@ -149,15 +163,16 @@ def layer_to_nodes(mpf, modellayer): if not isinstance(modellayer, (list, tuple)): modellayer = [modellayer] nodes = [] - node = 0 - for lay in range(mpf.ib.shape[0]): - for icell2d in range(mpf.ib.shape[1]): - # only add specific layers - if lay in modellayer: + for lay in modellayer: + if len(mpf.ib.shape) == 3: + for row in range(mpf.ib.shape[1]): + for col in range(mpf.ib.shape[2]): + if mpf.ib[lay, row, col] > 0: + nodes.append(get_node_structured(lay, row, col, mpf.ib.shape)) + else: + for icell2d in range(mpf.ib.shape[1]): if mpf.ib[lay, icell2d] > 0: - nodes.append(node) - node += 1 - + nodes.append(get_node_vertex(lay, icell2d, mpf.ib.shape)) return nodes diff --git a/nlmod/read/waterboard.py b/nlmod/read/waterboard.py index ff30970d..597d662b 100644 --- a/nlmod/read/waterboard.py +++ b/nlmod/read/waterboard.py @@ -196,7 +196,7 @@ def get_configuration(): config["Hollands Noorderkwartier"] = { "bgt_code": "W0651", "watercourses": { - "url": "https://kaarten.hhnk.nl/arcgis/rest/services/od_legger/od_legger_wateren_2022_oppervlaktewateren_ti/MapServer", + "url": "https://kaarten.hhnk.nl/arcgis/rest/services/od_legger/od_legger_wateren_2023_oppervlaktewateren_vg/MapServer/", "bottom_height": "WS_BODEMHOOGTE", }, "level_areas": { diff --git a/tests/test_001_model.py b/tests/test_001_model.py index d73130a0..3349357a 100644 --- a/tests/test_001_model.py +++ b/tests/test_001_model.py @@ -181,7 +181,7 @@ def test_create_sea_model(tmpdir): _ = nlmod.gwf.dis(ds, gwf) # create node property flow - _ = nlmod.gwf.npf(ds, gwf) + _ = nlmod.gwf.npf(ds, gwf, save_flows=True) # Create the initial conditions package _ = nlmod.gwf.ic(ds, gwf, starting_head=1.0) diff --git a/tests/test_002_regis_geotop.py b/tests/test_002_regis_geotop.py index 81bf1bce..77d93dc3 100644 --- a/tests/test_002_regis_geotop.py +++ b/tests/test_002_regis_geotop.py @@ -48,3 +48,20 @@ def test_get_regis_geotop_keep_all_layers( extent, use_regis=True, use_geotop=True, remove_nan_layers=False ) assert regis_geotop_ds.dims["layer"] == 137 + + +def test_add_kh_and_kv(): + gt = nlmod.read.geotop.get_geotop( + [118200, 118300, 439800, 439900], probabilities=True + ) + + # test with a value for kh for each lithoclass + df = nlmod.read.geotop.get_lithok_props() + # stochastic = None/False is allready tested in methods above, so we onlt test stochastic=True + gt = nlmod.read.geotop.add_kh_and_kv(gt, df, stochastic=True) + + # test with a value for kh and kv for each combination of lithoclass and stratigraphic unit + df = nlmod.read.geotop.get_kh_kv_table() + gt = nlmod.read.geotop.add_kh_and_kv(gt, df) + # again, but using the stochastic method + gt = nlmod.read.geotop.add_kh_and_kv(gt, df, stochastic=True) diff --git a/tests/test_012_plot.py b/tests/test_012_plot.py index 7dd08e80..9eb3279d 100644 --- a/tests/test_012_plot.py +++ b/tests/test_012_plot.py @@ -33,6 +33,8 @@ def test_plot_get_map(): def test_map_array(): ds = util.get_ds_structured() + nlmod.plot.map_array(ds["kh"], ds, ilay=0) + gwf = util.get_gwf(ds) nlmod.plot.flopy.map_array(ds["kh"], gwf, ilay=0) diff --git a/tests/test_015_gwf_output.py b/tests/test_015_gwf_output.py index 8083bee5..4bbe224c 100644 --- a/tests/test_015_gwf_output.py +++ b/tests/test_015_gwf_output.py @@ -190,7 +190,12 @@ def test_get_budget_da_from_file_vertex_with_grb(): nlmod.gwf.output.get_budget_da("CHD", fname=fname_cbc, grbfile=grbfile) -def test_gxg(): +def test_postprocess_head(): ds = test_001_model.get_ds_from_cache("basic_sea_model") head = nlmod.gwf.get_heads_da(ds) + nlmod.gwf.calculate_gxg(head) + + nlmod.gwf.get_gwl_from_wet_cells(head, botm=ds["botm"]) + + nlmod.gwf.get_head_at_point(head, float(ds.x.mean()), float(ds.y.mean()), ds=ds) diff --git a/tests/test_022_gwt.py b/tests/test_022_gwt.py index aaa017b2..3f864c4f 100644 --- a/tests/test_022_gwt.py +++ b/tests/test_022_gwt.py @@ -1,8 +1,137 @@ +import tempfile +import os +import pandas as pd +import xarray as xr import nlmod -def test_prepare(): - ds = nlmod.get_ds([0, 1000, 0, 2000]) +def test_gwt_model(): + extent = [103700, 106700, 527500, 528500] + + tmpdir = tempfile.gettempdir() + model_name = "trnsprt_tst" + model_ws = os.path.join(tmpdir, model_name) + + layer_model = nlmod.read.get_regis(extent, botm_layer="MSz1") + # create a model ds + ds = nlmod.to_model_ds( + layer_model, + model_name, + model_ws, + delr=100.0, + delc=100.0, + transport=True, + ) + + # add time discretisation + time = pd.date_range("2000", "2023", freq="YS") + ds = nlmod.time.set_ds_time(ds, start=time[0], time=time[1:]) + ds = nlmod.gwt.prepare.set_default_transport_parameters( ds, transport_type="chloride" ) + + # We download the digital terrain model (AHN4) + ahn = nlmod.read.ahn.get_ahn4(ds.extent) + # calculate the average surface level in each cell + ds["ahn"] = nlmod.resample.structured_da_to_ds(ahn, ds, method="average") + + # we then determine the part of each cell that is covered by sea from the original fine ahn + ds["sea"] = nlmod.read.rws.calculate_sea_coverage(ahn, ds=ds, method="average") + + # download knmi recharge data + knmi_ds = nlmod.read.knmi.get_recharge(ds) + + # update model dataset + ds.update(knmi_ds) + + # create simulation + sim = nlmod.sim.sim(ds) + + # create time discretisation + nlmod.sim.tdis(ds, sim) + + # create ims + ims = nlmod.sim.ims(sim, complexity="MODERATE") + + # create groundwater flow model + gwf = nlmod.gwf.gwf(ds, sim) + + # Create discretization + nlmod.gwf.dis(ds, gwf) + + # create node property flow + nlmod.gwf.npf(ds, gwf) + + # create storage + nlmod.gwf.sto(ds, gwf) + + # Create the initial conditions package + nlmod.gwf.ic(ds, gwf, starting_head=1.0) + + # Create the output control package + nlmod.gwf.oc(ds, gwf) + + # build ghb package + nlmod.gwf.ghb( + ds, + gwf, + bhead=0.0, + cond=ds["sea"] * ds["area"] / 1.0, + auxiliary=18_000.0, + ) + + # build surface level drain package + elev = ds["ahn"].where(ds["sea"] == 0) + nlmod.gwf.surface_drain_from_ds(ds, gwf, elev=elev, resistance=10.0) + + # create recharge package + nlmod.gwf.rch(ds, gwf, mask=ds["sea"] == 0) + + # BUY: buoyancy package for GWF model + nlmod.gwf.buy(ds, gwf) + + # GWT: groundwater transport model + gwt = nlmod.gwt.gwt(ds, sim) + + # add IMS for GWT model and register it + ims = nlmod.sim.ims(sim, pname="ims_gwt", filename=f"{gwt.name}.ims") + nlmod.sim.register_ims_package(sim, gwt, ims) + + # DIS: discretization package + nlmod.gwt.dis(ds, gwt) + + # IC: initial conditions package + nlmod.gwt.ic(ds, gwt, strt=18_000) + + # ADV: advection package + nlmod.gwt.adv(ds, gwt) + + # DSP: dispersion package + nlmod.gwt.dsp(ds, gwt) + + # MST: mass transfer package + nlmod.gwt.mst(ds, gwt) + + # SSM: source-sink mixing package + nlmod.gwt.ssm(ds, gwt) + + # OC: output control + nlmod.gwt.oc(ds, gwt) + + # GWF-GWT Exchange + nlmod.gwt.gwfgwt(ds, sim) + + nlmod.sim.write_and_run(sim, ds) + + h = nlmod.gwf.output.get_heads_da(ds) + c = nlmod.gwt.output.get_concentration_da(ds) + + # calculate concentration at groundwater surface + nlmod.gwt.get_concentration_at_gw_surface(c) + + # Convert calculated heads to equivalent freshwater heads, and vice versa + hf = nlmod.gwt.output.freshwater_head(ds, h, c) + hp = nlmod.gwt.output.pointwater_head(ds, hf, c) + + xr.testing.assert_allclose(h, hp) diff --git a/tests/test_025_modpath.py b/tests/test_025_modpath.py new file mode 100644 index 00000000..6e22c1a6 --- /dev/null +++ b/tests/test_025_modpath.py @@ -0,0 +1,45 @@ +import os +import xarray as xr +import flopy +import nlmod + + +def test_modpath(): + # start with runned model from test_001_model.test_create_sea_model + model_ws = os.path.join(os.path.dirname(os.path.realpath(__file__)), "data") + ds = xr.open_dataset(os.path.join(model_ws, "basic_sea_model.nc")) + + sim = flopy.mf6.MFSimulation.load("mfsim.nam", sim_ws=ds.model_ws) + gwf = sim.get_model(ds.model_name) + + xy_start = [ + (float(ds.x.mean()), float(ds.y.mean())), + (float(ds.x.mean()) + 100, float(ds.y.mean())), + ] + + # create a modpath model + mpf = nlmod.modpath.mpf(gwf) + + # create the basic modpath package + nlmod.modpath.bas(mpf) + + # find the nodes for given xy + nodes = nlmod.modpath.xy_to_nodes(xy_start, mpf, ds, layer=5) + + # create a particle tracking group at the cell faces + pg = nlmod.modpath.pg_from_fdt(nodes) + + # create the modpath simulation file + nlmod.modpath.sim(mpf, pg, "backward", gwf=gwf) + + # run modpath model + nlmod.modpath.write_and_run(mpf) + + # test reading of pathline file + nlmod.modpath.load_pathline_data(mpf) + + # get the nodes from a package + nodes = nlmod.modpath.package_to_nodes(gwf, "GHB", mpf) + + # get nodes of all cells in the top modellayer + nodes = nlmod.modpath.layer_to_nodes(mpf, 0)