From 66800d939409be723e34b62609b2950c06760480 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 8 Jul 2026 16:32:23 +0200 Subject: [PATCH 01/38] Boyscouting: Get rid of calls to .values --- imod/mf6/wel.py | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py index e39af906a..5ebec41de 100644 --- a/imod/mf6/wel.py +++ b/imod/mf6/wel.py @@ -61,11 +61,11 @@ def _assign_dims(arg: Any) -> tuple[Any, ...] | xr.DataArray: if arg.dims[0] != "time": arg = arg.transpose() da = xr.DataArray( - data=arg.values, coords={"time": arg["time"]}, dims=["time", "index"] + data=arg.to_numpy(), coords={"time": arg["time"]}, dims=["time", "index"] ) return da elif is_da: - return "index", arg.values + return "index", arg.to_numpy() else: return "index", arg @@ -347,11 +347,11 @@ class GridAgnosticWell(BoundaryCondition, IPointDataPackage, abc.ABC): @property def x(self) -> npt.NDArray[np.float64]: - return self.dataset["x"].values + return self.dataset["x"].to_numpy() @property def y(self) -> npt.NDArray[np.float64]: - return self.dataset["y"].values + return self.dataset["y"].to_numpy() @classmethod def _is_grid_agnostic_package(cls) -> bool: @@ -394,7 +394,9 @@ def _create_dataset_vars( # Carefully rename the dimension and set coordinates d_rename = {"index": "ncellid"} ds_vars = ds_vars.rename_dims(**d_rename).rename_vars(**d_rename) - ds_vars = ds_vars.assign_coords(**{"ncellid": cellid.coords["ncellid"].values}) + ds_vars = ds_vars.assign_coords( + **{"ncellid": cellid.coords["ncellid"].to_numpy()} + ) return ds_vars @@ -525,9 +527,9 @@ def _to_mf6_pkg( ds = ds.assign(**data_vars_dict) # type: ignore[arg-type] ds = remove_inactive(ds, idomain) - ds["save_flows"] = self["save_flows"].values[()] - ds["print_flows"] = self["print_flows"].values[()] - ds["print_input"] = self["print_input"].values[()] + ds["save_flows"] = enforce_scalar(self["save_flows"]) + ds["print_flows"] = enforce_scalar(self["print_flows"]) + ds["print_input"] = enforce_scalar(self["print_input"]) filtered_final_well_ids = self._gather_filtered_well_ids(ds, wells_df) if len(filtered_final_well_ids) > 0: @@ -1060,8 +1062,8 @@ def _find_well_value_at_layer( if (value is not None) and is_spatial_grid(value): value = imod.select.points_values( value, - x=well_dataset["x"].values, - y=well_dataset["y"].values, + x=well_dataset["x"].to_numpy(), + y=well_dataset["y"].to_numpy(), out_of_bounds="ignore", ) in_bounds = np.full(well_dataset.sizes["index"], False) From 0473920565aac87aab777b236e82d5dd2b8a51a0 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Thu, 9 Jul 2026 15:17:59 +0200 Subject: [PATCH 02/38] Pass id on to Mf6Wel --- imod/mf6/mf6_wel_adapter.py | 4 +++- imod/mf6/wel.py | 2 -- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/imod/mf6/mf6_wel_adapter.py b/imod/mf6/mf6_wel_adapter.py index c536e6612..23866bf0e 100644 --- a/imod/mf6/mf6_wel_adapter.py +++ b/imod/mf6/mf6_wel_adapter.py @@ -143,6 +143,7 @@ def __init__( self, cellid, rate, + id, concentration=None, concentration_boundary_type="aux", save_flows: Optional[bool] = None, @@ -153,6 +154,7 @@ def __init__( dict_dataset = { "cellid": cellid, "rate": rate, + "id": id, "concentration": concentration, "concentration_boundary_type": concentration_boundary_type, "save_flows": save_flows, @@ -169,7 +171,7 @@ def _ds_to_arrdict(self, ds): arrdict: Dict[str, Any] = {} arrdict["data_vars"] = [ - var_name for var_name in ds.data_vars if var_name != "cellid" + var_name for var_name in ds.data_vars if var_name not in ("cellid", "id") ] dsvar = {} diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py index 5ebec41de..30a201233 100644 --- a/imod/mf6/wel.py +++ b/imod/mf6/wel.py @@ -539,8 +539,6 @@ def _to_mf6_pkg( ) logger.log(loglevel=LogLevel.WARNING, message=message_end) - ds = ds.drop_vars("id") - data_vars_dict = {str(k): v for k, v in ds.data_vars.items()} return Mf6Wel(**data_vars_dict) # type: ignore[arg-type] From 7040e132f38d1386008f8a5ded7e0494c083434f Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Mon, 13 Jul 2026 14:29:18 +0200 Subject: [PATCH 03/38] Support IPF wells for LayeredWell.from_imod5_cap_data --- imod/mf6/utilities/imod5_converter.py | 20 ++++++++++++++++---- imod/mf6/wel.py | 13 ++++++++++++- 2 files changed, 28 insertions(+), 5 deletions(-) diff --git a/imod/mf6/utilities/imod5_converter.py b/imod/mf6/utilities/imod5_converter.py index adee28a6e..e27c21581 100644 --- a/imod/mf6/utilities/imod5_converter.py +++ b/imod/mf6/utilities/imod5_converter.py @@ -58,9 +58,16 @@ def fill_missing_layers( def _well_from_imod5_cap_point_data(cap_data: GridDataDict) -> dict[str, np.ndarray]: - raise NotImplementedError( - "Assigning sprinkling wells with an IPF file is not supported, please specify them as IDF." - ) + df_points = cap_data["artificial_recharge_layer"] + data = {} + # Order of columns is x, y, layer, the other columns are irrelevant here. + data["x"] = df_points.iloc[:, 0].to_numpy().astype(float) + data["y"] = df_points.iloc[:, 1].to_numpy().astype(float) + data["layer"] = df_points.iloc[:, 2].to_numpy().astype(int) + data["rate"] = np.zeros_like(data["x"], dtype=float) + data["id"] = df_points.index.to_numpy() + + return data def _well_from_imod5_cap_grid_data(cap_data: GridDataDict) -> dict[str, np.ndarray]: @@ -85,7 +92,7 @@ def _well_from_imod5_cap_grid_data(cap_data: GridDataDict) -> dict[str, np.ndarr def well_from_imod5_cap_data( imod5_data: Imod5DataDict, - target_dis: IRegridPackage, + target_dis: Optional[IRegridPackage], regridder_types: DataclassType, regrid_cache: RegridderWeightsCache, ) -> dict[str, np.ndarray]: @@ -121,6 +128,11 @@ def well_from_imod5_cap_data( if has_ipf_well: return _well_from_imod5_cap_point_data(cap_data) else: + if target_dis is None: + raise ValueError( + "target_dis must be provided when converting iMOD5 cap data " + "from grids (IDF)" + ) cap_data_regridded = regrid_imod5_cap_data( imod5_data, target_dis, regridder_types, regrid_cache )["cap"] diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py index 30a201233..abd280628 100644 --- a/imod/mf6/wel.py +++ b/imod/mf6/wel.py @@ -1449,7 +1449,7 @@ def _validate_imod5_depth_information( def from_imod5_cap_data( cls, imod5_data: Imod5DataDict, - target_dis: StructuredDiscretization, + target_dis: Optional[StructuredDiscretization] = None, regridder_types: CapDataWellRegridMethod = CapDataWellRegridMethod(), regrid_cache: RegridderWeightsCache = RegridderWeightsCache(), ): @@ -1489,6 +1489,17 @@ def from_imod5_cap_data( xarray datasets, under the key of the package type to which it belongs, as returned by :func:`imod.formats.prj.open_projectfile_data`. + target_dis: Optional[StructuredDiscretization] + The target discretization to which the data should be regridded. + Only necessary when "artificial_recharge_layer" is an IDF grid, + otherwise ignored. + regridder_types: CapDataWellRegridMethod + The regridder type to use for the regridding of the "artificial_recharge_layer" + and "artificial_recharge_capacity" grids. Only necessary when + "artificial_recharge_layer" is an IDF grid, otherwise ignored. + regrid_cache: RegridderWeightsCache + Cache for storing intermediate regridding results. Only necessary when + "artificial_recharge_layer" is an IDF grid, otherwise ignored. """ data = well_from_imod5_cap_data( imod5_data, target_dis, regridder_types, regrid_cache From aa1ab088e0acff2014e0056f61db48488610355d Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 15 Jul 2026 16:12:59 +0200 Subject: [PATCH 04/38] Temp add scratch script to try out data reworking --- imod/tests/_scratch.py | 120 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 120 insertions(+) create mode 100644 imod/tests/_scratch.py diff --git a/imod/tests/_scratch.py b/imod/tests/_scratch.py new file mode 100644 index 000000000..8d3e986b4 --- /dev/null +++ b/imod/tests/_scratch.py @@ -0,0 +1,120 @@ +# %% +import imod +import xarray as xr +from imod.util.dims import drop_layer_dim_cap_data + + +# %% +SUBUNIT = 0 + +df_points = imod.ipf.read(r"c:\Users\engelen\projects_wdir\imod-python\imod5_converter\NHI_sprint\Peelvenen\BASIS7\METASWAP\grid\Sprinkling\BEREGEN_LOC.IPF") +art_grid = imod.idf.open(r"c:\Users\engelen\projects_wdir\imod-python\imod5_converter\NHI_sprint\Peelvenen\BASIS7\METASWAP\grid\Sprinkling\BEREGENINGS_LOCATIES.IDF").compute().astype(int) + +imod5_data = {"cap": {"artificial_recharge": art_grid, "artificial_recharge_layer": df_points}} + +well = imod.mf6.LayeredWell.from_imod5_cap_data( + imod5_data, target_dis=None, regridder_types=None, regrid_cache=None +) +# %% +prj_data, period_data = imod.formats.prj.open_projectfile_data(r"c:\Users\engelen\projects_wdir\imod-python\imod5_converter\NHI_sprint\Peelvenen\prjfiles\Peelvenen_relative_paths_dis_npf.PRJ") + +# %% +dis_pkg = imod.mf6.StructuredDiscretization.from_imod5_data(prj_data, validate=False) +dis_pkg["idomain"] = dis_pkg["idomain"].clip(min=0) +npf_pkg = imod.mf6.NodePropertyFlow.from_imod5_data(prj_data, dis_pkg.dataset["idomain"]) + +prj_data = drop_layer_dim_cap_data(prj_data) +griddata, msw_active = imod.msw.GridData.from_imod5_data(prj_data, dis_pkg) +# %% +# Convert to args of Sprinkling._render() +mf6_well = well.to_mf6_pkg(dis_pkg["idomain"], dis_pkg["top"], dis_pkg["bottom"], npf_pkg["k"]) +isactive_1d, svat = griddata.generate_isactive_svat_arrays() + +# %% +# In from_imod5_cap_data +arl_points = df_points.iloc[:, :5] +arl_points.columns = ["x_p", "y_p", "layer_p", "id2grid_p", "capacity"] +# Enforce dtypes +arl_points = arl_points.astype({"x_p": float, "y_p": float, "layer_p": int, "id2grid_p": int, "capacity": float}) +arl_points["id"] = arl_points.index.astype(str) +arl_points = arl_points.set_index("id") + +points_ds = arl_points.to_xarray() +art_grid = art_grid.rename("id_msw") +dataset = xr.merge([art_grid, points_ds]) + +# %% +# Flatten id_msw grid → (y, x, id_msw) table, drop cells with no well +# TODO: Verify with Peter that only the first subunit (landuse: agriculture) is +# relevant for the SVAT mapping. +grids = xr.merge([art_grid, svat.sel(subunit=SUBUNIT, drop=True)]) + +art_df = ( + grids + .to_dataframe() + .reset_index() + .query("(id_msw > 0) & (svat > 0)") +) +# Drop unnecessary columns. We preserve the x, y coords as they might prove +# useful for debugging. +art_df = art_df.drop(["dx", "dy"], axis=1) + +# Join: each SVAT cell gets the matching well row(s) from arl_points +msw_merged_df = art_df.merge( + arl_points.reset_index(), # brings id back as a column + left_on="id_msw", + right_on="id2grid_p", + how="inner", + validate="many_to_one" +) + +# %% +# Promote id to dim to join datasets +mf6_well_ds = mf6_well.dataset.set_coords("id").swap_dims({"ncellid": "id"}) +# Convert the cellid DataArray to a broad table for easier manipulation. +mf6_cellid_df = mf6_well_ds["cellid"].to_dataset("dim_cellid").to_dataframe() +# Select only the cellid columns we need and reset index to promote id to column +# for merging +dim_cellid = ["layer", "row", "column"] +mf6_cellid_df = mf6_cellid_df.loc[:, dim_cellid].reset_index() +# Merge again to confine to wells actually used in the modflow6 model. +# This drops points that are not in art_grid. +msw_mf6_merged_df = msw_merged_df.merge( + mf6_cellid_df, + on="id", + how="inner", + validate="many_to_one" +) + +# %% +# Derive the SVAT subunit for each well from the SVAT grid based on the +# well's row/col location. +indexer = msw_mf6_merged_df.loc[:, dim_cellid].to_numpy() - 1 +# We need to align the SVAT grid with the dis_pkg grid as the SVAT grid might be +# smaller. +idomain_flat = dis_pkg.dataset["idomain"].isel(layer=0, drop=True) +_, svat_aligned = xr.align(idomain_flat, svat, join="left") +svat_groundwater = svat_aligned.data[SUBUNIT, indexer[:, 1], indexer[:, 2]] + +msw_mf6_merged_df["svat_groundwater"] = svat_groundwater.astype(int) + +# %% +# Deal with edge case: wells that are outside art_grid, but in model domain +outside_art = ~mf6_cellid_df["id"].isin(msw_mf6_merged_df["id"].unique()) +indexer_outside = mf6_cellid_df.loc[outside_art, dim_cellid].to_numpy() - 1 +svat_outside = svat_aligned.data[SUBUNIT, indexer_outside[:, 1], indexer_outside[:, 2]] + +# %% +# Select columns that need to be written to scap_svat.inp +dataframe = msw_mf6_merged_df[["svat", "layer", "svat_groundwater"]] +dataframe["svat"] = dataframe["svat"].astype(int) +capacity = msw_mf6_merged_df["capacity"] +is_gw_extraction = msw_mf6_merged_df["layer"] > 0 +dataframe["max_abstraction_groundwater"] = capacity.where(is_gw_extraction, 0.0) +dataframe["max_abstraction_surfacewater"] = capacity.where(~is_gw_extraction, 0.0) + +# Cases to catch: +# 1. Well is outside the model domain → surfacewater +# 2. Well is in layer 0 → surfacewater + + From 6405be447b34fcd41e1184809ada799f8383b0ce Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Thu, 16 Jul 2026 15:15:02 +0200 Subject: [PATCH 05/38] Change order of merging to first merge tables, then to grid and deal with edge case --- imod/tests/_scratch.py | 77 ++++++++++++++++++++++++------------------ 1 file changed, 45 insertions(+), 32 deletions(-) diff --git a/imod/tests/_scratch.py b/imod/tests/_scratch.py index 8d3e986b4..56b3dba8e 100644 --- a/imod/tests/_scratch.py +++ b/imod/tests/_scratch.py @@ -2,7 +2,7 @@ import imod import xarray as xr from imod.util.dims import drop_layer_dim_cap_data - +import pandas as pd # %% SUBUNIT = 0 @@ -16,16 +16,15 @@ imod5_data, target_dis=None, regridder_types=None, regrid_cache=None ) # %% +# Setup to get example args for _render() of SprinklingPoints prj_data, period_data = imod.formats.prj.open_projectfile_data(r"c:\Users\engelen\projects_wdir\imod-python\imod5_converter\NHI_sprint\Peelvenen\prjfiles\Peelvenen_relative_paths_dis_npf.PRJ") -# %% dis_pkg = imod.mf6.StructuredDiscretization.from_imod5_data(prj_data, validate=False) dis_pkg["idomain"] = dis_pkg["idomain"].clip(min=0) npf_pkg = imod.mf6.NodePropertyFlow.from_imod5_data(prj_data, dis_pkg.dataset["idomain"]) prj_data = drop_layer_dim_cap_data(prj_data) griddata, msw_active = imod.msw.GridData.from_imod5_data(prj_data, dis_pkg) -# %% # Convert to args of Sprinkling._render() mf6_well = well.to_mf6_pkg(dis_pkg["idomain"], dis_pkg["top"], dis_pkg["bottom"], npf_pkg["k"]) isactive_1d, svat = griddata.generate_isactive_svat_arrays() @@ -43,6 +42,23 @@ art_grid = art_grid.rename("id_msw") dataset = xr.merge([art_grid, points_ds]) +# %% +# Promote id to dim to join datasets +mf6_well_ds = mf6_well.dataset.set_coords("id").swap_dims({"ncellid": "id"}) +# Convert the cellid DataArray to a broad table for easier manipulation. +mf6_cellid_df = mf6_well_ds["cellid"].to_dataset("dim_cellid").to_dataframe() +# Select only the cellid columns we need and reset index to promote id to column +# for merging +dim_cellid = ["layer", "row", "column"] +mf6_cellid_df = mf6_cellid_df.loc[:, dim_cellid].reset_index() +# Merge again to confine to wells actually used in the modflow6 model. +# This drops points that are not in art_grid. +points_mf6_merged_df = arl_points.reset_index().merge( + mf6_cellid_df, + on="id", + how="right", + validate="many_to_one" +) # %% # Flatten id_msw grid → (y, x, id_msw) table, drop cells with no well # TODO: Verify with Peter that only the first subunit (landuse: agriculture) is @@ -60,32 +76,14 @@ art_df = art_df.drop(["dx", "dy"], axis=1) # Join: each SVAT cell gets the matching well row(s) from arl_points -msw_merged_df = art_df.merge( - arl_points.reset_index(), # brings id back as a column +msw_mf6_merged_df = art_df.merge( + points_mf6_merged_df, # brings id back as a column left_on="id_msw", right_on="id2grid_p", how="inner", validate="many_to_one" ) -# %% -# Promote id to dim to join datasets -mf6_well_ds = mf6_well.dataset.set_coords("id").swap_dims({"ncellid": "id"}) -# Convert the cellid DataArray to a broad table for easier manipulation. -mf6_cellid_df = mf6_well_ds["cellid"].to_dataset("dim_cellid").to_dataframe() -# Select only the cellid columns we need and reset index to promote id to column -# for merging -dim_cellid = ["layer", "row", "column"] -mf6_cellid_df = mf6_cellid_df.loc[:, dim_cellid].reset_index() -# Merge again to confine to wells actually used in the modflow6 model. -# This drops points that are not in art_grid. -msw_mf6_merged_df = msw_merged_df.merge( - mf6_cellid_df, - on="id", - how="inner", - validate="many_to_one" -) - # %% # Derive the SVAT subunit for each well from the SVAT grid based on the # well's row/col location. @@ -98,23 +96,38 @@ msw_mf6_merged_df["svat_groundwater"] = svat_groundwater.astype(int) -# %% -# Deal with edge case: wells that are outside art_grid, but in model domain -outside_art = ~mf6_cellid_df["id"].isin(msw_mf6_merged_df["id"].unique()) -indexer_outside = mf6_cellid_df.loc[outside_art, dim_cellid].to_numpy() - 1 -svat_outside = svat_aligned.data[SUBUNIT, indexer_outside[:, 1], indexer_outside[:, 2]] - # %% # Select columns that need to be written to scap_svat.inp dataframe = msw_mf6_merged_df[["svat", "layer", "svat_groundwater"]] dataframe["svat"] = dataframe["svat"].astype(int) capacity = msw_mf6_merged_df["capacity"] +# Set wells with layer > 0 to groundwater abstraction, and wells with layer = 0 +# to surfacewater abstraction. is_gw_extraction = msw_mf6_merged_df["layer"] > 0 dataframe["max_abstraction_groundwater"] = capacity.where(is_gw_extraction, 0.0) dataframe["max_abstraction_surfacewater"] = capacity.where(~is_gw_extraction, 0.0) -# Cases to catch: -# 1. Well is outside the model domain → surfacewater -# 2. Well is in layer 0 → surfacewater +# %% +# +# Deal with edge case: wells that are outside art_grid, but in model domain. +# These will be assigned to surfacewater abstraction. We can identify these by +# checking which wells in mf6_cellid_df are not in msw_mf6_merged_df. +is_outside_art = ~mf6_cellid_df["id"].isin(msw_mf6_merged_df["id"].unique()) +outside_df = points_mf6_merged_df.loc[is_outside_art, ["layer", "capacity"]] +# Select the SVAT subunit for these wells based on their row/col location. +indexer_outside = mf6_cellid_df.loc[is_outside_art, dim_cellid].to_numpy() - 1 +svat_outside = svat_aligned.data[SUBUNIT, indexer_outside[:, 1], indexer_outside[:, 2]] +outside_df["svat_groundwater"] = svat_outside.astype(int) +outside_df["svat"] = svat_outside.astype(int) +# Set capacity to surfacewater abstraction, and set groundwater abstraction to 0. +outside_df = outside_df.rename(columns={"capacity": "max_abstraction_surfacewater"}) +outside_df["max_abstraction_groundwater"] = 0.0 +outside_df = outside_df.query("svat > 0").reset_index(drop=True) +# %% +# +# Combine +dataframe_out = pd.concat([dataframe, outside_df], axis=0, ignore_index=True) +dataframe_out = dataframe_out.sort_values(by=["svat"]).reset_index(drop=True) +# %% From fe251dd0163580a135326c7c5f7ae0f3b3967c2b Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Thu, 16 Jul 2026 16:57:42 +0200 Subject: [PATCH 06/38] Make script work for both subunits --- imod/tests/_scratch.py | 69 ++++++++++++++++++++++++++++++++++-------- 1 file changed, 57 insertions(+), 12 deletions(-) diff --git a/imod/tests/_scratch.py b/imod/tests/_scratch.py index 56b3dba8e..048f3fc4a 100644 --- a/imod/tests/_scratch.py +++ b/imod/tests/_scratch.py @@ -4,9 +4,50 @@ from imod.util.dims import drop_layer_dim_cap_data import pandas as pd -# %% -SUBUNIT = 0 +def get_indexer(df: pd.DataFrame, columns: list[str]): + """ + Get the indexer for a dataframe of wells to select the SVAT subunit for each + well based on its row/col location in the model grid. + + Parameters + ---------- + df : pd.DataFrame + DataFrame containing the wells with columns "subunit", "row", + and "column" (1-based). + + Returns + ------- + np.ndarray + Indexer array for selecting SVAT subunits from svat_da. + """ + if not set(["row", "column"]).issubset(columns): + raise ValueError("columns must contain 'row' and 'column'") + df.loc[:, ["row", "column"]] -= 1 # Convert to 0-based indexing for xarray + + indexer = df.loc[:, columns].to_numpy() + return indexer.T + + +def double_length_df_subunit(df: pd.DataFrame, subunit_col: str = "subunit") -> pd.DataFrame: + """ + Duplicate the rows of a DataFrame for each subunit (0 and 1). + + Parameters + ---------- + df : pd.DataFrame + Input DataFrame to duplicate. + subunit_col : str, optional + Name of the column to assign subunit values, by default "subunit". + + Returns + ------- + pd.DataFrame + DataFrame with duplicated rows for each subunit. + """ + return pd.concat([df.assign(**{subunit_col: 0}), df.assign(**{subunit_col: 1})], ignore_index=True) + +# %% df_points = imod.ipf.read(r"c:\Users\engelen\projects_wdir\imod-python\imod5_converter\NHI_sprint\Peelvenen\BASIS7\METASWAP\grid\Sprinkling\BEREGEN_LOC.IPF") art_grid = imod.idf.open(r"c:\Users\engelen\projects_wdir\imod-python\imod5_converter\NHI_sprint\Peelvenen\BASIS7\METASWAP\grid\Sprinkling\BEREGENINGS_LOCATIES.IDF").compute().astype(int) @@ -61,10 +102,7 @@ ) # %% # Flatten id_msw grid → (y, x, id_msw) table, drop cells with no well -# TODO: Verify with Peter that only the first subunit (landuse: agriculture) is -# relevant for the SVAT mapping. -grids = xr.merge([art_grid, svat.sel(subunit=SUBUNIT, drop=True)]) - +grids = xr.merge([art_grid, svat]) art_df = ( grids .to_dataframe() @@ -84,15 +122,15 @@ validate="many_to_one" ) -# %% +# %% # Derive the SVAT subunit for each well from the SVAT grid based on the # well's row/col location. -indexer = msw_mf6_merged_df.loc[:, dim_cellid].to_numpy() - 1 +indexer = get_indexer(msw_mf6_merged_df, columns=["subunit", "row", "column"]) # We need to align the SVAT grid with the dis_pkg grid as the SVAT grid might be # smaller. idomain_flat = dis_pkg.dataset["idomain"].isel(layer=0, drop=True) _, svat_aligned = xr.align(idomain_flat, svat, join="left") -svat_groundwater = svat_aligned.data[SUBUNIT, indexer[:, 1], indexer[:, 2]] +svat_groundwater = svat_aligned.data[*indexer] msw_mf6_merged_df["svat_groundwater"] = svat_groundwater.astype(int) @@ -114,15 +152,20 @@ # checking which wells in mf6_cellid_df are not in msw_mf6_merged_df. is_outside_art = ~mf6_cellid_df["id"].isin(msw_mf6_merged_df["id"].unique()) outside_df = points_mf6_merged_df.loc[is_outside_art, ["layer", "capacity"]] +outside_df = double_length_df_subunit(outside_df) # Select the SVAT subunit for these wells based on their row/col location. -indexer_outside = mf6_cellid_df.loc[is_outside_art, dim_cellid].to_numpy() - 1 -svat_outside = svat_aligned.data[SUBUNIT, indexer_outside[:, 1], indexer_outside[:, 2]] +df_cellid_outside = mf6_cellid_df.loc[is_outside_art] +df_cellid_outside = double_length_df_subunit(df_cellid_outside) +indexer_outside = get_indexer(df_cellid_outside, columns=["subunit", "row", "column"]) +svat_outside = svat_aligned.data[indexer_outside[0], indexer_outside[1], indexer_outside[2]] outside_df["svat_groundwater"] = svat_outside.astype(int) outside_df["svat"] = svat_outside.astype(int) # Set capacity to surfacewater abstraction, and set groundwater abstraction to 0. outside_df = outside_df.rename(columns={"capacity": "max_abstraction_surfacewater"}) outside_df["max_abstraction_groundwater"] = 0.0 - +# drop subunit column as it is no longer needed +outside_df = outside_df.drop(columns=["subunit"]) +# drop wells that are outside the active metaswap model domain (svat = 0) outside_df = outside_df.query("svat > 0").reset_index(drop=True) # %% @@ -131,3 +174,5 @@ dataframe_out = pd.concat([dataframe, outside_df], axis=0, ignore_index=True) dataframe_out = dataframe_out.sort_values(by=["svat"]).reset_index(drop=True) # %% +# +# TODO: Verify if iMOD5 SVAT grid the same as the one derived in this script. \ No newline at end of file From 42948e5545e1a0735d8bf323fcfc3eddac10616c Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 22 Jul 2026 09:43:48 +0200 Subject: [PATCH 07/38] Prototype SprinklingPoints class --- imod/msw/regrid/regrid_schemes.py | 22 ++ imod/msw/sprinkling.py | 342 +++++++++++++++++++++++++++++- 2 files changed, 362 insertions(+), 2 deletions(-) diff --git a/imod/msw/regrid/regrid_schemes.py b/imod/msw/regrid/regrid_schemes.py index cf82018b1..084e78841 100644 --- a/imod/msw/regrid/regrid_schemes.py +++ b/imod/msw/regrid/regrid_schemes.py @@ -104,6 +104,28 @@ class SprinklingRegridMethod(DataclassType): max_abstraction_surfacewater: RegridVarType = (RegridderType.OVERLAP, "mean") +@dataclass(config=_CONFIG) +class SprinklingPointsRegridMethod(DataclassType): + """ + Object containing regridder methods for the + :class:`imod.msw.Sprinkling` package. This can be provided to the + ``regrid_like`` method to regrid with custom settings. + + Parameters + ---------- + art_grid: tuple, default (RegridderType.OVERLAP, "mode") + + Examples + -------- + Regrid with custom settings: + + >>> regrid_method = SprinklingPointsRegridMethod(art_grid=(RegridderType.OVERLAP,"min")) + >>> sprinking.regrid_like(target_grid, RegridderWeightsCache(), regrid_method) + """ + + art_grid: RegridVarType = (RegridderType.OVERLAP, "mode") + + @dataclass(config=_CONFIG) class MeteoGridRegridMethod(DataclassType): """ diff --git a/imod/msw/sprinkling.py b/imod/msw/sprinkling.py index 802aec550..6be532f94 100644 --- a/imod/msw/sprinkling.py +++ b/imod/msw/sprinkling.py @@ -9,7 +9,10 @@ from imod.mf6.mf6_wel_adapter import Mf6Wel from imod.msw.fixed_format import VariableMetaData from imod.msw.pkgbase import MetaSwapPackage -from imod.msw.regrid.regrid_schemes import SprinklingRegridMethod +from imod.msw.regrid.regrid_schemes import ( + SprinklingPointsRegridMethod, + SprinklingRegridMethod, +) from imod.msw.utilities.common import concat_imod5 from imod.msw.utilities.imod5_converter import ( get_cell_area_from_imod5_data, @@ -65,10 +68,181 @@ def _sprinkling_data_from_imod5_grid(cap_data: GridDataDict) -> GridDataDict: return data +def _extract_indexer_for_svat(df: pd.DataFrame, columns: list[str]): + """ + Get the indexer for a dataframe of wells to select the SVAT subunit for each + well based on its row/col location in the model grid. + + Parameters + ---------- + df : pd.DataFrame + DataFrame containing the wells with columns "subunit", "row", + and "column". "row" and "column" are 1-based indices. + columns : list[str] + List of column names to use for indexing. Must include "row" and "column". + + Returns + ------- + np.ndarray + Indexer array for selecting SVAT subunits from svat_da. + """ + if not set("row", "column").issubset(columns): + raise ValueError("columns must contain 'row' and 'column'") + df.loc[:, ["row", "column"]] -= 1 # Convert to 0-based indexing for xarray + + indexer = df.loc[:, columns].to_numpy() + return indexer.T + + +def _replicate_dataframe_by_subunit( + df: pd.DataFrame, subunit_col: str = "subunit" +) -> pd.DataFrame: + """ + Duplicate the rows of a DataFrame for each subunit (0 and 1). + + Parameters + ---------- + df : pd.DataFrame + Input DataFrame to duplicate. + subunit_col : str, optional + Name of the column to assign subunit values, by default "subunit". + + Returns + ------- + pd.DataFrame + DataFrame with duplicated rows for each subunit. + """ + subunit_nrs = [0, 1] + df_ls = [df.assign(**{subunit_col: subunit_nr}) for subunit_nr in subunit_nrs] + return pd.concat(df_ls, ignore_index=True) + + +def _get_mf6_cellid_dataframe(mf6_well: Mf6Wel) -> pd.DataFrame: + """ + Get cellids from the Mf6Wel objects dataset and convert to a dataframe for + easy merging with sprinkling data. + """ + # Promote id to dim to join datasets + mf6_well_ds = mf6_well.dataset.set_coords("id").swap_dims({"ncellid": "id"}) + # Convert the cellid DataArray to a broad table for easier manipulation. + mf6_cellid_df = mf6_well_ds["cellid"].to_dataset("dim_cellid").to_dataframe() + # Select only the cellid columns we need and reset index to promote id to column + # for merging + dim_cellid = ["layer", "row", "column"] + mf6_cellid_df = mf6_cellid_df.loc[:, dim_cellid].reset_index() + return mf6_cellid_df + + +def _make_sprinkling_well_points_dataframe( + sprinkling_dataset: xr.Dataset, mf6_cellid_df: pd.DataFrame +) -> pd.DataFrame: + """ + Create a dataframe of sprinkling well points from the sprinkling dataset and + merge it with the mf6_cellid_df to get the row/col of each well. + """ + # Get point data from sprinkling dataset and convert to dataframe for easy merging + points_keys = [ + key for key, da in sprinkling_dataset.data_vars.items() if "id" in da.dims + ] + sprinkling_points_df = ( + sprinkling_dataset[points_keys].drop_vars(["dx", "dy"]).to_dataframe() + ) + # Merge again to confine to wells actually used in the modflow6 model. + # This drops points that are outside model domain. + return sprinkling_points_df.reset_index().merge( + mf6_cellid_df, on="id", how="right", validate="many_to_one" + ) + + +def _merge_sprinkling_points_with_grids( + points_df: pd.DataFrame, svat: xr.DataArray, sprinkling_id_grid: xr.DataArray +) -> pd.DataFrame: + """ + Merge sprinkling points with SVAT grids. + """ + + # TODO: Rename "id_msw" and "id2grid_p" to something clearer like "id_sprinkling" + # Flatten id_msw grid → (y, x, id_msw) table, drop cells with no well + grids = xr.merge([sprinkling_id_grid, svat]) + art_df = grids.to_dataframe().reset_index().query("(id_msw > 0) & (svat > 0)") + # Drop unnecessary columns. We preserve the x, y coords as they might + # prove useful for debugging. + art_df = art_df.drop(["dx", "dy"], axis=1) + + # Join: each SVAT cell gets the matching well row(s) from arl_points + return art_df.merge( + points_df, # brings id back as a column + left_on="id_msw", + right_on="id2grid_p", + how="inner", + validate="many_to_one", + ) + + +def align_svat_with_dis( + svat: xr.DataArray, dis_pkg: StructuredDiscretization +) -> xr.DataArray: + """ + Align the SVAT grid with the dis_pkg grid as the SVAT grid might be smaller. + """ + idomain_flat = dis_pkg.dataset["idomain"].isel(layer=0, drop=True) + _, svat_aligned = xr.align(idomain_flat, svat, join="left") + return svat_aligned + + +def _get_svat_groundwater_for_wells( + msw_mf6_sprinkling_df: pd.DataFrame, svat_aligned: xr.DataArray +) -> np.ndarray: + """ + Get the SVAT subunit for each well from the SVAT grid based on the + well's row/col location. + """ + indexer = _extract_indexer_for_svat( + msw_mf6_sprinkling_df, columns=["subunit", "row", "column"] + ) + svat_groundwater = svat_aligned.data[*indexer] + return svat_groundwater.astype(int) + + +def _get_wells_outside_art_grid_dataframe( + mf6_cellid_df: pd.DataFrame, + points_df: pd.DataFrame, + msw_mf6_sprinkling_df: pd.DataFrame, + svat_aligned: xr.DataArray, +) -> pd.DataFrame: + """ + Deal with edge case: wells that are outside art_grid, but in model domain. + These will be assigned to surfacewater abstraction. We can identify these by + checking which wells in mf6_cellid_df are not in msw_mf6_merged_df. + """ + is_outside_art = ~mf6_cellid_df["id"].isin(msw_mf6_sprinkling_df["id"].unique()) + outside_df = points_df.loc[is_outside_art, ["layer", "capacity"]] + outside_df = _replicate_dataframe_by_subunit(outside_df) + # Select the SVAT subunit for these wells based on their row/col location. + cellid_outside_df = mf6_cellid_df.loc[is_outside_art] + cellid_outside_df = _replicate_dataframe_by_subunit(cellid_outside_df) + indexer_outside = _extract_indexer_for_svat( + cellid_outside_df, columns=["subunit", "row", "column"] + ) + svat_outside = svat_aligned.data[*indexer_outside] + outside_df["svat_groundwater"] = svat_outside.astype(int) + outside_df["svat"] = svat_outside.astype(int) + # Set capacity to surfacewater abstraction, and set groundwater abstraction to 0. + outside_df = outside_df.rename(columns={"capacity": "max_abstraction_surfacewater"}) + outside_df["max_abstraction_groundwater"] = 0.0 + # drop subunit column as it is no longer needed + outside_df = outside_df.drop(columns=["subunit"]) + # drop wells that are outside the active metaswap model domain (svat = 0) + return outside_df.query("svat > 0").reset_index(drop=True) + + class Sprinkling(MetaSwapPackage, IRegridPackage): """ This contains the sprinkling capacities of links between SVAT units and - groundwater/surface water locations. + groundwater/surface water locations. Input is provided as grids for the + maximum abstraction of groundwater and surfacewater to SVAT units. To + specify the sprinkling capacity as points, see + :class:`imod.msw.SprinklingPoints`. This class is responsible for the file `scap_svat.inp` @@ -221,3 +395,167 @@ def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "Sprinkling": data = _sprinkling_data_from_imod5_grid(cap_data) return cls(**data) + + +class SprinklingPoints(MetaSwapPackage, IRegridPackage): + """ + This contains the sprinkling capacities of links between SVAT units and + groundwater/surface water locations. This class is capable of handling point + data (IPF) for sprinkling wells, which is a mapping of grid cells to well + locations. To specify the sprinkling capacity as grid, see + :class:`imod.msw.Sprinkling`. + + This class is responsible for the file `scap_svat.inp` + + Parameters + ---------- + art_grid: xr.DataArray + Grid of the artificial recharge types, with subunit coordinate. + x_p: np.ndarray | list[float] + x-coordinates of the artificial recharge locations. + y_p: np.ndarray | list[float] + y-coordinates of the artificial recharge locations. + layer_p: np.ndarray | list[int] + layer indices of the artificial recharge locations. + id2grid_p: np.ndarray | list[int] + mapping of the artificial recharge locations to the grid cells. + capacity_p: np.ndarray | list[float] + abstraction capacities of the artificial recharge locations. + + """ + + _file_name = "scap_svat.inp" + _metadata_dict = { + "svat": VariableMetaData(10, 1, 99999999, int), + "max_abstraction_groundwater_mm_d": VariableMetaData(8, None, None, str), + "max_abstraction_surfacewater_mm_d": VariableMetaData(8, None, None, str), + "max_abstraction_groundwater": VariableMetaData(8, 0.0, 1e9, float), + "max_abstraction_surfacewater": VariableMetaData(8, 0.0, 1e9, float), + "svat_groundwater": VariableMetaData(10, 1, 99999999, int), + "layer": VariableMetaData(6, 1, 9999, int), + "trajectory": VariableMetaData(10, None, None, str), + } + + _with_subunit = ( + "max_abstraction_groundwater", + "max_abstraction_surfacewater", + ) + _without_subunit = () + + _to_fill = ( + "max_abstraction_groundwater_mm_d", + "max_abstraction_surfacewater_mm_d", + "trajectory", + ) + + _regrid_method = SprinklingPointsRegridMethod() + + def __init__( + self, + art_grid: xr.DataArray, + x_p: np.ndarray | list[float], + y_p: np.ndarray | list[float], + layer_p: np.ndarray | list[int], + id2grid_p: np.ndarray | list[int], + capacity_p: np.ndarray | list[float], + ): + super().__init__() + # Replicate well ids as they were also created in + # imod.mf6.LayeredWell.from_imod5_cap_data() + id_index = pd.Index(range(len(x_p)), name="id").astype(str) + points_ds = xr.Dataset( + { + "x_p": (("id",), x_p), + "y_p": (("id",), y_p), + "layer_p": (("id",), layer_p), + "id2grid_p": (("id",), id2grid_p), + "capacity_p": (("id",), capacity_p), + }, + coords={"id": id_index}, + ) + art_grid = art_grid.rename("id_msw") + self.dataset = xr.merge([art_grid, points_ds]) + + @classmethod + def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "SprinklingPoints": + cap_data = imod5_data["cap"] + art_grid = cap_data["artificial_recharge"] + df_points = cap_data["artificial_recharge_layer"] + + arl_points = df_points.iloc[:, :5] + arl_points.columns = ["x_p", "y_p", "layer_p", "id2grid_p", "capacity_p"] + # Enforce dtypes + arl_points = arl_points.astype( + { + "x_p": float, + "y_p": float, + "layer_p": int, + "id2grid_p": int, + "capacity_p": float, + } + ) + + return cls( + art_grid=art_grid, + x_p=arl_points["x_p"].to_numpy(), + y_p=arl_points["y_p"].to_numpy(), + layer_p=arl_points["layer_p"].to_numpy(), + id2grid_p=arl_points["id2grid_p"].to_numpy(), + capacity_p=arl_points["capacity_p"].to_numpy(), + ) + + def _render(self, file, index, svat, mf6_dis, mf6_well): + """ + Render the sprinkling points to the scap_svat.inp file. + + This method first merges the sprinkling points with the mf6_well cellids + to get the row/col of each well, then merges the sprinkling points with + the svat and id_msw grid. It then selects the columns that need to be + written to scap_svat.inp and sets wells with layer > 0 to groundwater + abstraction, and wells with layer = 0 to surfacewater abstraction. + Finally, it deals with edge cases for wells that are outside art_grid + but in the model domain, and writes the dataframe to the file. + """ + # Merge the sprinkling points with the mf6_well cellids to get the + # row/col of each well. + mf6_cellid_df = _get_mf6_cellid_dataframe(mf6_well) + points_df = _make_sprinkling_well_points_dataframe(self.dataset, mf6_cellid_df) + # Merge the sprinkling points with the svat and id_msw grid + msw_mf6_sprinkling_df = _merge_sprinkling_points_with_grids( + points_df, svat, self.dataset["id_msw"] + ) + svat_aligned = align_svat_with_dis(svat, mf6_dis) + msw_mf6_sprinkling_df["svat_groundwater"] = _get_svat_groundwater_for_wells( + msw_mf6_sprinkling_df, svat_aligned + ) + + # Select columns that need to be written to scap_svat.inp + dataframe = msw_mf6_sprinkling_df[["svat", "layer", "svat_groundwater"]] + dataframe["svat"] = dataframe["svat"].astype(int) + capacity = msw_mf6_sprinkling_df["capacity"] + # Set wells with layer > 0 to groundwater abstraction, and wells with layer = 0 + # to surfacewater abstraction. + is_gw_extraction = msw_mf6_sprinkling_df["layer"] > 0 + dataframe["max_abstraction_groundwater"] = capacity.where(is_gw_extraction, 0.0) + dataframe["max_abstraction_surfacewater"] = capacity.where( + ~is_gw_extraction, 0.0 + ) + + # TODO: Make sure wells in svats are all present in the dataframe. If + # these svats are 0 in the art_grid, they should get a 0.0 capacity. + + # Deal with edge case: wells that are outside art_grid, but in model domain. + # These will be assigned to surfacewater abstraction. + outside_df = _get_wells_outside_art_grid_dataframe( + mf6_cellid_df, points_df, msw_mf6_sprinkling_df, svat_aligned + ) + + dataframe_out = pd.concat([dataframe, outside_df], axis=0, ignore_index=True) + dataframe_out = dataframe_out.sort_values(by=["svat"]).reset_index(drop=True) + + for var in self._to_fill: + dataframe_out[var] = "" + + self._check_range(dataframe_out) + + return self._write_dataframe_fixed_width(file, dataframe_out) From 87a5e02a2a91efc4a6419972f2cbfdb42a5fed8b Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 22 Jul 2026 09:43:58 +0200 Subject: [PATCH 08/38] Format --- imod/tests/_scratch.py | 82 ++++++++++++++++++++++++++---------------- 1 file changed, 52 insertions(+), 30 deletions(-) diff --git a/imod/tests/_scratch.py b/imod/tests/_scratch.py index 048f3fc4a..be0f2dba2 100644 --- a/imod/tests/_scratch.py +++ b/imod/tests/_scratch.py @@ -1,8 +1,10 @@ # %% -import imod +import pandas as pd import xarray as xr + +import imod from imod.util.dims import drop_layer_dim_cap_data -import pandas as pd + def get_indexer(df: pd.DataFrame, columns: list[str]): """ @@ -20,7 +22,7 @@ def get_indexer(df: pd.DataFrame, columns: list[str]): np.ndarray Indexer array for selecting SVAT subunits from svat_da. """ - if not set(["row", "column"]).issubset(columns): + if not set("row", "column").issubset(columns): raise ValueError("columns must contain 'row' and 'column'") df.loc[:, ["row", "column"]] -= 1 # Convert to 0-based indexing for xarray @@ -28,7 +30,9 @@ def get_indexer(df: pd.DataFrame, columns: list[str]): return indexer.T -def double_length_df_subunit(df: pd.DataFrame, subunit_col: str = "subunit") -> pd.DataFrame: +def double_length_df_subunit( + df: pd.DataFrame, subunit_col: str = "subunit" +) -> pd.DataFrame: """ Duplicate the rows of a DataFrame for each subunit (0 and 1). @@ -44,30 +48,49 @@ def double_length_df_subunit(df: pd.DataFrame, subunit_col: str = "subunit") -> pd.DataFrame DataFrame with duplicated rows for each subunit. """ - return pd.concat([df.assign(**{subunit_col: 0}), df.assign(**{subunit_col: 1})], ignore_index=True) + return pd.concat( + [df.assign(**{subunit_col: 0}), df.assign(**{subunit_col: 1})], + ignore_index=True, + ) # %% -df_points = imod.ipf.read(r"c:\Users\engelen\projects_wdir\imod-python\imod5_converter\NHI_sprint\Peelvenen\BASIS7\METASWAP\grid\Sprinkling\BEREGEN_LOC.IPF") -art_grid = imod.idf.open(r"c:\Users\engelen\projects_wdir\imod-python\imod5_converter\NHI_sprint\Peelvenen\BASIS7\METASWAP\grid\Sprinkling\BEREGENINGS_LOCATIES.IDF").compute().astype(int) +df_points = imod.ipf.read( + r"c:\Users\engelen\projects_wdir\imod-python\imod5_converter\NHI_sprint\Peelvenen\BASIS7\METASWAP\grid\Sprinkling\BEREGEN_LOC.IPF" +) +art_grid = ( + imod.idf.open( + r"c:\Users\engelen\projects_wdir\imod-python\imod5_converter\NHI_sprint\Peelvenen\BASIS7\METASWAP\grid\Sprinkling\BEREGENINGS_LOCATIES.IDF" + ) + .compute() + .astype(int) +) -imod5_data = {"cap": {"artificial_recharge": art_grid, "artificial_recharge_layer": df_points}} +imod5_data = { + "cap": {"artificial_recharge": art_grid, "artificial_recharge_layer": df_points} +} well = imod.mf6.LayeredWell.from_imod5_cap_data( imod5_data, target_dis=None, regridder_types=None, regrid_cache=None ) # %% # Setup to get example args for _render() of SprinklingPoints -prj_data, period_data = imod.formats.prj.open_projectfile_data(r"c:\Users\engelen\projects_wdir\imod-python\imod5_converter\NHI_sprint\Peelvenen\prjfiles\Peelvenen_relative_paths_dis_npf.PRJ") +prj_data, period_data = imod.formats.prj.open_projectfile_data( + r"c:\Users\engelen\projects_wdir\imod-python\imod5_converter\NHI_sprint\Peelvenen\prjfiles\Peelvenen_relative_paths_dis_npf.PRJ" +) dis_pkg = imod.mf6.StructuredDiscretization.from_imod5_data(prj_data, validate=False) dis_pkg["idomain"] = dis_pkg["idomain"].clip(min=0) -npf_pkg = imod.mf6.NodePropertyFlow.from_imod5_data(prj_data, dis_pkg.dataset["idomain"]) +npf_pkg = imod.mf6.NodePropertyFlow.from_imod5_data( + prj_data, dis_pkg.dataset["idomain"] +) prj_data = drop_layer_dim_cap_data(prj_data) griddata, msw_active = imod.msw.GridData.from_imod5_data(prj_data, dis_pkg) # Convert to args of Sprinkling._render() -mf6_well = well.to_mf6_pkg(dis_pkg["idomain"], dis_pkg["top"], dis_pkg["bottom"], npf_pkg["k"]) +mf6_well = well.to_mf6_pkg( + dis_pkg["idomain"], dis_pkg["top"], dis_pkg["bottom"], npf_pkg["k"] +) isactive_1d, svat = griddata.generate_isactive_svat_arrays() # %% @@ -75,15 +98,20 @@ def double_length_df_subunit(df: pd.DataFrame, subunit_col: str = "subunit") -> arl_points = df_points.iloc[:, :5] arl_points.columns = ["x_p", "y_p", "layer_p", "id2grid_p", "capacity"] # Enforce dtypes -arl_points = arl_points.astype({"x_p": float, "y_p": float, "layer_p": int, "id2grid_p": int, "capacity": float}) +arl_points = arl_points.astype( + {"x_p": float, "y_p": float, "layer_p": int, "id2grid_p": int, "capacity": float} +) arl_points["id"] = arl_points.index.astype(str) arl_points = arl_points.set_index("id") points_ds = arl_points.to_xarray() + +# in def __init__ art_grid = art_grid.rename("id_msw") dataset = xr.merge([art_grid, points_ds]) # %% +# In render() # Promote id to dim to join datasets mf6_well_ds = mf6_well.dataset.set_coords("id").swap_dims({"ncellid": "id"}) # Convert the cellid DataArray to a broad table for easier manipulation. @@ -95,36 +123,28 @@ def double_length_df_subunit(df: pd.DataFrame, subunit_col: str = "subunit") -> # Merge again to confine to wells actually used in the modflow6 model. # This drops points that are not in art_grid. points_mf6_merged_df = arl_points.reset_index().merge( - mf6_cellid_df, - on="id", - how="right", - validate="many_to_one" + mf6_cellid_df, on="id", how="right", validate="many_to_one" ) # %% # Flatten id_msw grid → (y, x, id_msw) table, drop cells with no well grids = xr.merge([art_grid, svat]) -art_df = ( - grids - .to_dataframe() - .reset_index() - .query("(id_msw > 0) & (svat > 0)") -) +art_df = grids.to_dataframe().reset_index().query("(id_msw > 0) & (svat > 0)") # Drop unnecessary columns. We preserve the x, y coords as they might prove # useful for debugging. art_df = art_df.drop(["dx", "dy"], axis=1) # Join: each SVAT cell gets the matching well row(s) from arl_points msw_mf6_merged_df = art_df.merge( - points_mf6_merged_df, # brings id back as a column + points_mf6_merged_df, # brings id back as a column left_on="id_msw", right_on="id2grid_p", how="inner", - validate="many_to_one" + validate="many_to_one", ) # %% # Derive the SVAT subunit for each well from the SVAT grid based on the -# well's row/col location. +# well's row/col location. indexer = get_indexer(msw_mf6_merged_df, columns=["subunit", "row", "column"]) # We need to align the SVAT grid with the dis_pkg grid as the SVAT grid might be # smaller. @@ -141,12 +161,12 @@ def double_length_df_subunit(df: pd.DataFrame, subunit_col: str = "subunit") -> capacity = msw_mf6_merged_df["capacity"] # Set wells with layer > 0 to groundwater abstraction, and wells with layer = 0 # to surfacewater abstraction. -is_gw_extraction = msw_mf6_merged_df["layer"] > 0 +is_gw_extraction = msw_mf6_merged_df["layer"] > 0 dataframe["max_abstraction_groundwater"] = capacity.where(is_gw_extraction, 0.0) dataframe["max_abstraction_surfacewater"] = capacity.where(~is_gw_extraction, 0.0) -# %% -# +# %% +# # Deal with edge case: wells that are outside art_grid, but in model domain. # These will be assigned to surfacewater abstraction. We can identify these by # checking which wells in mf6_cellid_df are not in msw_mf6_merged_df. @@ -157,7 +177,9 @@ def double_length_df_subunit(df: pd.DataFrame, subunit_col: str = "subunit") -> df_cellid_outside = mf6_cellid_df.loc[is_outside_art] df_cellid_outside = double_length_df_subunit(df_cellid_outside) indexer_outside = get_indexer(df_cellid_outside, columns=["subunit", "row", "column"]) -svat_outside = svat_aligned.data[indexer_outside[0], indexer_outside[1], indexer_outside[2]] +svat_outside = svat_aligned.data[ + indexer_outside[0], indexer_outside[1], indexer_outside[2] +] outside_df["svat_groundwater"] = svat_outside.astype(int) outside_df["svat"] = svat_outside.astype(int) # Set capacity to surfacewater abstraction, and set groundwater abstraction to 0. @@ -175,4 +197,4 @@ def double_length_df_subunit(df: pd.DataFrame, subunit_col: str = "subunit") -> dataframe_out = dataframe_out.sort_values(by=["svat"]).reset_index(drop=True) # %% # -# TODO: Verify if iMOD5 SVAT grid the same as the one derived in this script. \ No newline at end of file +# TODO: Verify if iMOD5 SVAT grid the same as the one derived in this script. From f088ad279b4582936a709ccc404ffde61d0b4127 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 12 Aug 2026 13:41:40 +0200 Subject: [PATCH 09/38] Fix mypy errors --- imod/msw/sprinkling.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/imod/msw/sprinkling.py b/imod/msw/sprinkling.py index 6be532f94..148f16ebc 100644 --- a/imod/msw/sprinkling.py +++ b/imod/msw/sprinkling.py @@ -86,7 +86,7 @@ def _extract_indexer_for_svat(df: pd.DataFrame, columns: list[str]): np.ndarray Indexer array for selecting SVAT subunits from svat_da. """ - if not set("row", "column").issubset(columns): + if not set(["row", "column"]).issubset(columns): raise ValueError("columns must contain 'row' and 'column'") df.loc[:, ["row", "column"]] -= 1 # Convert to 0-based indexing for xarray @@ -186,7 +186,9 @@ def align_svat_with_dis( Align the SVAT grid with the dis_pkg grid as the SVAT grid might be smaller. """ idomain_flat = dis_pkg.dataset["idomain"].isel(layer=0, drop=True) - _, svat_aligned = xr.align(idomain_flat, svat, join="left") + # Assign to _dummy instead of _ to avoid MyPy 2.3 crashing on the next line. + # See https://github.com/python/mypy/issues/21824 + _dummy, svat_aligned = xr.align(idomain_flat, svat, join="left") return svat_aligned From eb8a9e2ae20effc3cba563694c4d8ce0a0e4a139 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 12 Aug 2026 13:45:48 +0200 Subject: [PATCH 10/38] Convert to set literal --- imod/msw/sprinkling.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/imod/msw/sprinkling.py b/imod/msw/sprinkling.py index 148f16ebc..6f065d99c 100644 --- a/imod/msw/sprinkling.py +++ b/imod/msw/sprinkling.py @@ -86,7 +86,7 @@ def _extract_indexer_for_svat(df: pd.DataFrame, columns: list[str]): np.ndarray Indexer array for selecting SVAT subunits from svat_da. """ - if not set(["row", "column"]).issubset(columns): + if not {"row", "column"}.issubset(columns): raise ValueError("columns must contain 'row' and 'column'") df.loc[:, ["row", "column"]] -= 1 # Convert to 0-based indexing for xarray From 308b1577727043744fb9ff173b609c6e22c2e910 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 12 Aug 2026 14:29:51 +0200 Subject: [PATCH 11/38] Clearer varname for inside_df --- imod/msw/sprinkling.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/imod/msw/sprinkling.py b/imod/msw/sprinkling.py index 6f065d99c..516100b94 100644 --- a/imod/msw/sprinkling.py +++ b/imod/msw/sprinkling.py @@ -532,18 +532,18 @@ def _render(self, file, index, svat, mf6_dis, mf6_well): ) # Select columns that need to be written to scap_svat.inp - dataframe = msw_mf6_sprinkling_df[["svat", "layer", "svat_groundwater"]] - dataframe["svat"] = dataframe["svat"].astype(int) + inside_df = msw_mf6_sprinkling_df[["svat", "layer", "svat_groundwater"]] + inside_df["svat"] = inside_df["svat"].astype(int) capacity = msw_mf6_sprinkling_df["capacity"] # Set wells with layer > 0 to groundwater abstraction, and wells with layer = 0 # to surfacewater abstraction. is_gw_extraction = msw_mf6_sprinkling_df["layer"] > 0 - dataframe["max_abstraction_groundwater"] = capacity.where(is_gw_extraction, 0.0) - dataframe["max_abstraction_surfacewater"] = capacity.where( + inside_df["max_abstraction_groundwater"] = capacity.where(is_gw_extraction, 0.0) + inside_df["max_abstraction_surfacewater"] = capacity.where( ~is_gw_extraction, 0.0 ) - # TODO: Make sure wells in svats are all present in the dataframe. If + # TODO: Make sure wells in svats are all present in the inside_df. If # these svats are 0 in the art_grid, they should get a 0.0 capacity. # Deal with edge case: wells that are outside art_grid, but in model domain. @@ -552,7 +552,7 @@ def _render(self, file, index, svat, mf6_dis, mf6_well): mf6_cellid_df, points_df, msw_mf6_sprinkling_df, svat_aligned ) - dataframe_out = pd.concat([dataframe, outside_df], axis=0, ignore_index=True) + dataframe_out = pd.concat([inside_df, outside_df], axis=0, ignore_index=True) dataframe_out = dataframe_out.sort_values(by=["svat"]).reset_index(drop=True) for var in self._to_fill: From 99dba8f5bbdda76a3aa239bb1829355f2084f91f Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 12 Aug 2026 15:27:03 +0200 Subject: [PATCH 12/38] Add comment --- imod/msw/sprinkling.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/imod/msw/sprinkling.py b/imod/msw/sprinkling.py index 516100b94..d4b6f6056 100644 --- a/imod/msw/sprinkling.py +++ b/imod/msw/sprinkling.py @@ -551,7 +551,7 @@ def _render(self, file, index, svat, mf6_dis, mf6_well): outside_df = _get_wells_outside_art_grid_dataframe( mf6_cellid_df, points_df, msw_mf6_sprinkling_df, svat_aligned ) - + # Prepare the final dataframe to be written to scap_svat.inp dataframe_out = pd.concat([inside_df, outside_df], axis=0, ignore_index=True) dataframe_out = dataframe_out.sort_values(by=["svat"]).reset_index(drop=True) From 157e91d98197524be3b6c7fd5782d696d69b4ca0 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 12 Aug 2026 15:27:23 +0200 Subject: [PATCH 13/38] Add SprinklingPoints to msw namespace --- imod/msw/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/imod/msw/__init__.py b/imod/msw/__init__.py index d8d053c49..6169964e6 100644 --- a/imod/msw/__init__.py +++ b/imod/msw/__init__.py @@ -16,5 +16,5 @@ from imod.msw.output_control import TimeOutputControl, VariableOutputControl from imod.msw.ponding import Ponding from imod.msw.scaling_factors import ScalingFactors -from imod.msw.sprinkling import Sprinkling +from imod.msw.sprinkling import Sprinkling, SprinklingPoints from imod.msw.vegetation import AnnualCropFactors From ab4815d7dcdbe11ccaa2d01cd257e0dd657aa1a9 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 12 Aug 2026 15:39:26 +0200 Subject: [PATCH 14/38] Put columns in right order --- imod/tests/fixtures/imod5_cap_data.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/imod/tests/fixtures/imod5_cap_data.py b/imod/tests/fixtures/imod5_cap_data.py index 88a2a5d7d..38c8e814d 100644 --- a/imod/tests/fixtures/imod5_cap_data.py +++ b/imod/tests/fixtures/imod5_cap_data.py @@ -146,11 +146,11 @@ def cap_data_sprinkling_points() -> Imod5DataDict: artificial_rch_type[:, 2] = 4000 data = { - "id": [3000, 4000], + "x": [1.5, 2.5], + "y": [1.5, 2.5], "layer": [2, 3], + "id": [3000, 4000], "capacity": [15.0, 30.0], - "y": [1.0, 2.0], - "x": [1.0, 2.0], } layer = pd.DataFrame(data=data) From e5af555565b155e1aed55029fd1fce002ac8eaed Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 12 Aug 2026 15:39:39 +0200 Subject: [PATCH 15/38] Start adding unittest --- imod/tests/test_msw/test_sprinkling.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/imod/tests/test_msw/test_sprinkling.py b/imod/tests/test_msw/test_sprinkling.py index 63812df9f..2350e9efc 100644 --- a/imod/tests/test_msw/test_sprinkling.py +++ b/imod/tests/test_msw/test_sprinkling.py @@ -400,3 +400,18 @@ def test_sprinkling_from_imod5_data__grid(cap_data_sprinkling_grid): np.testing.assert_array_equal( rural_ds["max_abstraction_surfacewater"].to_numpy(), expected_sw_abstraction ) + +@pytest.mark.unittest_jit +def test_sprinklingpoints_from_imod5_data__points(cap_data_sprinkling_points): + # Arrange + expected_vars = {'id2grid_p', 'capacity_p', 'layer_p', 'y_p', 'x_p', 'id_msw'} + rate = 25.0 * 0.25 * 1.0e-3 + + # Act + sprinkling = msw.SprinklingPoints.from_imod5_data(cap_data_sprinkling_points) + + # Assert + assert sprinkling.dataset.sizes == {"id": 2, "x": 3, "y": 3} + assert set(sprinkling.dataset.keys()) == expected_vars + + From 01f4a5e535e4e587a23c7e1e592965054f0154cf Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 12 Aug 2026 17:08:25 +0200 Subject: [PATCH 16/38] Start refactoring test_module: Separate cases from tests --- imod/tests/test_msw/test_sprinkling.py | 317 ++++++++++--------------- 1 file changed, 127 insertions(+), 190 deletions(-) diff --git a/imod/tests/test_msw/test_sprinkling.py b/imod/tests/test_msw/test_sprinkling.py index 2350e9efc..f217038b8 100644 --- a/imod/tests/test_msw/test_sprinkling.py +++ b/imod/tests/test_msw/test_sprinkling.py @@ -1,54 +1,28 @@ import tempfile +from dataclasses import dataclass from pathlib import Path +from typing import Callable, Optional import numpy as np import pytest import xarray as xr from numpy import nan from numpy.testing import assert_almost_equal, assert_equal +from pytest_cases import parametrize_with_cases from imod import msw from imod.mf6.mf6_wel_adapter import Mf6Wel from imod.mf6.wel import derive_cellid_from_points -def test_simple_model_all_svats(fixed_format_parser): +@pytest.fixture(scope="function") +def sprinkling_svat_index(): x = [1.0, 2.0, 3.0] y = [1.0, 2.0, 3.0] subunit = [0, 1] dx = 1.0 dy = 1.0 # fmt: off - max_abstraction_groundwater = xr.DataArray( - np.array( - [ - [[nan, 100.0, nan], - [nan, 200.0, nan], - [nan, 300.0, nan]], - [[nan, 100.0, nan], - [nan, 200.0, nan], - [nan, 300.0, nan]] - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} - ) - - max_abstraction_surfacewater = xr.DataArray( - np.array( - [ - [[nan, 100.0, nan], - [nan, 200.0, nan], - [nan, 300.0, nan]], - [[nan, 100.0, nan], - [nan, 200.0, nan], - [nan, 300.0, nan]] - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} - ) - svat = xr.DataArray( np.array( [ @@ -66,66 +40,60 @@ def test_simple_model_all_svats(fixed_format_parser): ) # fmt: on index = (svat != 0).values.ravel() - - # Well - well_layer = [3, 2, 1] - well_y = y - well_x = [2.0, 2.0, 2.0] - well_rate = [-5.0] * 3 - cellids = derive_cellid_from_points(svat, well_x, well_y, well_layer) - well = Mf6Wel(cellids, well_rate) - - sprinkling = msw.Sprinkling( - max_abstraction_groundwater, - max_abstraction_surfacewater, - ) - - with tempfile.TemporaryDirectory() as output_dir: - output_dir = Path(output_dir) - sprinkling.write(output_dir, index, svat, None, well) - - results = fixed_format_parser( - output_dir / msw.Sprinkling._file_name, - msw.Sprinkling._metadata_dict, + return svat, index + + +@dataclass +class AbstractionCaseData: + max_abstraction_groundwater: Optional[xr.DataArray] = None + max_abstraction_surfacewater: Optional[xr.DataArray] = None + expected_abs_gw: Optional[np.ndarray] = None + expected_abs_sw: Optional[np.ndarray] = None + expected_layer: Optional[np.ndarray] = None + expected_svat_gw: Optional[np.ndarray] = None + + +class AbstractionGrids: + def case_all_svats(self, sprinkling_svat_index) -> AbstractionCaseData: + svat, _ = sprinkling_svat_index + case_data = AbstractionCaseData() + case_data.max_abstraction_groundwater = xr.full_like(svat, 0.0) + case_data.max_abstraction_surfacewater = xr.full_like(svat, 0.0) + # fmt: off + case_data.max_abstraction_groundwater.data = np.array( + [ + [[nan, 100.0, nan], + [nan, 200.0, nan], + [nan, 300.0, nan]], + [[nan, 100.0, nan], + [nan, 200.0, nan], + [nan, 300.0, nan]] + ] ) - - assert_equal(results["svat"], np.array([1, 2, 3, 4])) - assert_almost_equal( - results["max_abstraction_groundwater"], - np.array([100.0, 300.0, 100.0, 200.0]), - ) - assert_almost_equal( - results["max_abstraction_surfacewater"], - np.array([100.0, 300.0, 100.0, 200.0]), - ) - assert_equal(results["layer"], np.array([3, 1, 3, 2])) - assert_equal(results["svat_groundwater"], np.array([1, 2, 3, 4])) - - -def test_simple_model_some_svats(fixed_format_parser): - x = [1.0, 2.0, 3.0] - y = [1.0, 2.0, 3.0] - subunit = [0, 1] - dx = 1.0 - dy = 1.0 - # fmt: off - max_abstraction_groundwater = xr.DataArray( - np.array( + case_data.max_abstraction_surfacewater.data = np.array( [ [[nan, 100.0, nan], [nan, 200.0, nan], [nan, 300.0, nan]], - [[nan, nan, nan], + [[nan, 100.0, nan], [nan, 200.0, nan], - [nan, nan, nan]] + [nan, 300.0, nan]] ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} - ) - - max_abstraction_surfacewater = xr.DataArray( - np.array( + ) + # fmt: on + case_data.expected_abs_gw = np.array([100.0, 300.0, 100.0, 200.0]) + case_data.expected_abs_sw = np.array([100.0, 300.0, 100.0, 200.0]) + case_data.expected_layer = np.array([3, 1, 3, 2]) + case_data.expected_svat_gw = np.array([1, 2, 3, 4]) + return case_data + + def case_some_svats(self, sprinkling_svat_index): + svat, _ = sprinkling_svat_index + case_data = AbstractionCaseData() + case_data.max_abstraction_groundwater = xr.full_like(svat, 0.0) + case_data.max_abstraction_surfacewater = xr.full_like(svat, 0.0) + # fmt: off + case_data.max_abstraction_groundwater.data = np.array( [ [[nan, 100.0, nan], [nan, 200.0, nan], @@ -134,73 +102,32 @@ def test_simple_model_some_svats(fixed_format_parser): [nan, 200.0, nan], [nan, nan, nan]] ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} - ) - - svat = xr.DataArray( - np.array( + ) + case_data.max_abstraction_surfacewater.data = np.array( [ - [[0, 1, 0], - [0, 0, 0], - [0, 2, 0]], - - [[0, 3, 0], - [0, 4, 0], - [0, 0, 0]], + [[nan, 100.0, nan], + [nan, 200.0, nan], + [nan, 300.0, nan]], + [[nan, nan, nan], + [nan, 200.0, nan], + [nan, nan, nan]] ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} - ) - # fmt: on - index = (svat != 0).values.ravel() - - # Well - well_layer = [3, 2, 1] - well_y = [1.0, 2.0, 3.0] - well_x = [2.0, 2.0, 2.0] - well_rate = [-5.0] * 3 - cellids = derive_cellid_from_points(svat, well_x, well_y, well_layer) - well = Mf6Wel(cellids, well_rate) - - coupler_mapping = msw.Sprinkling( - max_abstraction_groundwater, - max_abstraction_surfacewater, - ) - - with tempfile.TemporaryDirectory() as output_dir: - output_dir = Path(output_dir) - coupler_mapping.write(output_dir, index, svat, None, well) - - results = fixed_format_parser( - output_dir / msw.Sprinkling._file_name, - msw.Sprinkling._metadata_dict, ) - - assert_equal(results["svat"], np.array([1, 2, 4])) - assert_almost_equal( - results["max_abstraction_groundwater"], - np.array([100.0, 300.0, 200.0]), - ) - assert_almost_equal( - results["max_abstraction_surfacewater"], - np.array([100.0, 300.0, 200.0]), - ) - assert_equal(results["layer"], np.array([3, 1, 2])) - assert_equal(results["svat_groundwater"], np.array([1, 2, 4])) - - -def test_simple_model_inconsistent_active_capacity(fixed_format_parser): - x = [1.0, 2.0, 3.0] - y = [1.0, 2.0, 3.0] - subunit = [0, 1] - dx = 1.0 - dy = 1.0 - # fmt: off - max_abstraction_groundwater = xr.DataArray( - np.array( + # fmt: on + case_data.expected_abs_gw = np.array([100.0, 300.0, 200.0]) + case_data.expected_abs_sw = np.array([100.0, 300.0, 200.0]) + case_data.expected_layer = np.array([3, 1, 2]) + case_data.expected_svat_gw = np.array([1, 2, 4]) + + return case_data + + def case_inconsistent_active_capacity(self, sprinkling_svat_index): + svat, _ = sprinkling_svat_index + case_data = AbstractionCaseData() + case_data.max_abstraction_groundwater = xr.full_like(svat, 0.0) + case_data.max_abstraction_surfacewater = xr.full_like(svat, 0.0) + # fmt: off + case_data.max_abstraction_groundwater.data = np.array( [ [[nan, 100.0, nan], [nan, 0.0, nan], @@ -209,13 +136,8 @@ def test_simple_model_inconsistent_active_capacity(fixed_format_parser): [nan, 200.0, nan], [nan, nan, nan]] ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} - ) - - max_abstraction_surfacewater = xr.DataArray( - np.array( + ) + case_data.max_abstraction_surfacewater.data = np.array( [ [[nan, 0.0, nan], [nan, 200.0, nan], @@ -224,62 +146,58 @@ def test_simple_model_inconsistent_active_capacity(fixed_format_parser): [nan, 200.0, nan], [nan, nan, nan]] ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} - ) + ) + # fmt: on + case_data.expected_abs_gw = np.array([100.0, 0.0, 200.0]) + case_data.expected_abs_sw = np.array([0.0, 300.0, 200.0]) + case_data.expected_layer = np.array([3, 1, 2]) + case_data.expected_svat_gw = np.array([1, 2, 4]) - svat = xr.DataArray( - np.array( - [ - [[0, 1, 0], - [0, 0, 0], - [0, 2, 0]], + return case_data - [[0, 3, 0], - [0, 4, 0], - [0, 0, 0]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} - ) - # fmt: on - index = (svat != 0).values.ravel() + +@parametrize_with_cases("case_data", cases=AbstractionGrids) +def test_simple_model( + fixed_format_parser: Callable, + sprinkling_svat_index: tuple[xr.DataArray, np.ndarray], + case_data: AbstractionCaseData, +): + svat, index = sprinkling_svat_index # Well well_layer = [3, 2, 1] well_y = [1.0, 2.0, 3.0] well_x = [2.0, 2.0, 2.0] well_rate = [-5.0] * 3 + well_id = ["a", "b", "c"] cellids = derive_cellid_from_points(svat, well_x, well_y, well_layer) - well = Mf6Wel(cellids, well_rate) + well = Mf6Wel(cellids, well_rate, well_id) - coupler_mapping = msw.Sprinkling( - max_abstraction_groundwater, - max_abstraction_surfacewater, + sprinkling = msw.Sprinkling( + case_data.max_abstraction_groundwater, + case_data.max_abstraction_surfacewater, ) with tempfile.TemporaryDirectory() as output_dir: output_dir = Path(output_dir) - coupler_mapping.write(output_dir, index, svat, None, well) + sprinkling.write(output_dir, index, svat, None, well) results = fixed_format_parser( output_dir / msw.Sprinkling._file_name, msw.Sprinkling._metadata_dict, ) - assert_equal(results["svat"], np.array([1, 2, 4])) + assert_equal(results["svat"], case_data.expected_svat_gw) assert_almost_equal( results["max_abstraction_groundwater"], - np.array([100.0, 0.0, 200.0]), + case_data.expected_abs_gw, ) assert_almost_equal( results["max_abstraction_surfacewater"], - np.array([0.0, 300.0, 200.0]), + case_data.expected_abs_sw, ) - assert_equal(results["layer"], np.array([3, 1, 2])) - assert_equal(results["svat_groundwater"], np.array([1, 2, 4])) + assert_equal(results["layer"], case_data.expected_layer) + assert_equal(results["svat_groundwater"], case_data.expected_svat_gw) def test_simple_model_1_subunit(fixed_format_parser): @@ -332,8 +250,9 @@ def test_simple_model_1_subunit(fixed_format_parser): well_y = [1.0, 3.0] well_x = [2.0, 2.0] well_rate = [-5.0] * 2 + well_id = ["a", "c"] cellids = derive_cellid_from_points(svat, well_x, well_y, well_layer) - well = Mf6Wel(cellids, well_rate) + well = Mf6Wel(cellids, well_rate, well_id) sprinkling = msw.Sprinkling( max_abstraction_groundwater, @@ -401,11 +320,11 @@ def test_sprinkling_from_imod5_data__grid(cap_data_sprinkling_grid): rural_ds["max_abstraction_surfacewater"].to_numpy(), expected_sw_abstraction ) + @pytest.mark.unittest_jit def test_sprinklingpoints_from_imod5_data__points(cap_data_sprinkling_points): # Arrange - expected_vars = {'id2grid_p', 'capacity_p', 'layer_p', 'y_p', 'x_p', 'id_msw'} - rate = 25.0 * 0.25 * 1.0e-3 + expected_vars = {"id2grid_p", "capacity_p", "layer_p", "y_p", "x_p", "id_msw"} # Act sprinkling = msw.SprinklingPoints.from_imod5_data(cap_data_sprinkling_points) @@ -413,5 +332,23 @@ def test_sprinklingpoints_from_imod5_data__points(cap_data_sprinkling_points): # Assert assert sprinkling.dataset.sizes == {"id": 2, "x": 3, "y": 3} assert set(sprinkling.dataset.keys()) == expected_vars - - + # No unit conversion is done in SprinklingPoints, as the capacity is already + # in m3/d + np.testing.assert_almost_equal(sprinkling.dataset["capacity_p"].max(), 100.0) + + +# @pytest.mark.unittest_jit +# def test_sprinklingpoints_write__points(cap_data_sprinkling_points, tmp_path): +# well_x = cap_data_sprinkling_points["cap"]["x"].values +# well_y = cap_data_sprinkling_points["cap"]["y"].values +# well_layer = cap_data_sprinkling_points["cap"]["layer"].values +# +# # Arrange +# # cellids = derive_cellid_from_points(svat, well_x, well_y, well_layer) +# # well = Mf6Wel(cellids, well_rate) +# +# # Act +# sprinkling = msw.SprinklingPoints.from_imod5_data(cap_data_sprinkling_points) +# +# # sprinkling.write() +# From 67058da9e8e4cee17240cdbb97840d3c5bdb2587 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Thu, 13 Aug 2026 11:18:33 +0200 Subject: [PATCH 17/38] Refactor into test for one subunit --- imod/tests/test_msw/test_sprinkling.py | 69 +++++++------------------- 1 file changed, 18 insertions(+), 51 deletions(-) diff --git a/imod/tests/test_msw/test_sprinkling.py b/imod/tests/test_msw/test_sprinkling.py index f217038b8..6e7845cc0 100644 --- a/imod/tests/test_msw/test_sprinkling.py +++ b/imod/tests/test_msw/test_sprinkling.py @@ -200,53 +200,19 @@ def test_simple_model( assert_equal(results["svat_groundwater"], case_data.expected_svat_gw) -def test_simple_model_1_subunit(fixed_format_parser): - x = [1.0, 2.0, 3.0] - y = [1.0, 2.0, 3.0] - subunit = [0] - dx = 1.0 - dy = 1.0 - # fmt: off - max_abstraction_groundwater = xr.DataArray( - np.array( - [ - [[nan, 100.0, nan], - [nan, 200.0, nan], - [nan, 300.0, nan]] - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} - ) - - max_abstraction_surfacewater = xr.DataArray( - np.array( - [ - [[nan, 100.0, nan], - [nan, 200.0, nan], - [nan, 300.0, nan]] - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} - ) +@parametrize_with_cases("case_data", cases=AbstractionGrids) +def test_simple_model_1_subunit( + fixed_format_parser: Callable, + sprinkling_svat_index: tuple[xr.DataArray, np.ndarray], + case_data: AbstractionCaseData, +): + svat, index = sprinkling_svat_index - svat = xr.DataArray( - np.array( - [ - [[0, 1, 0], - [0, 0, 0], - [0, 2, 0]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} - ) - # fmt: on - index = (svat != 0).values.ravel() + svat = svat.isel(subunit=[0]) + index = index[:9] # Only the first subunit # Well - well_layer = [3, 2] + well_layer = [3, 1] well_y = [1.0, 3.0] well_x = [2.0, 2.0] well_rate = [-5.0] * 2 @@ -255,8 +221,8 @@ def test_simple_model_1_subunit(fixed_format_parser): well = Mf6Wel(cellids, well_rate, well_id) sprinkling = msw.Sprinkling( - max_abstraction_groundwater, - max_abstraction_surfacewater, + case_data.max_abstraction_groundwater.isel(subunit=[0]), + case_data.max_abstraction_surfacewater.isel(subunit=[0]), ) with tempfile.TemporaryDirectory() as output_dir: @@ -268,17 +234,18 @@ def test_simple_model_1_subunit(fixed_format_parser): msw.Sprinkling._metadata_dict, ) - assert_equal(results["svat"], np.array([1, 2])) + assert_equal(results["svat"], case_data.expected_svat_gw[:2]) assert_almost_equal( results["max_abstraction_groundwater"], - np.array([100.0, 300.0]), + case_data.expected_abs_gw[:2], ) assert_almost_equal( results["max_abstraction_surfacewater"], - np.array([100.0, 300.0]), + case_data.expected_abs_sw[:2], ) - assert_equal(results["layer"], np.array([3, 2])) - assert_equal(results["svat_groundwater"], np.array([1, 2])) + assert_equal(results["layer"], case_data.expected_layer[:2]) + assert_equal(results["svat_groundwater"], case_data.expected_svat_gw[:2]) + @pytest.mark.unittest_jit From cb7458a6a6507741669e573a4f9a1be9e9226b57 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Thu, 13 Aug 2026 13:22:00 +0200 Subject: [PATCH 18/38] Add SprinklingPoints to the public API --- docs/api/msw.rst | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/docs/api/msw.rst b/docs/api/msw.rst index a14f60f59..2f7a1194a 100644 --- a/docs/api/msw.rst +++ b/docs/api/msw.rst @@ -56,6 +56,12 @@ Grid packages Sprinkling.from_imod5_data Sprinkling.get_regrid_methods Sprinkling.write + SprinklingPoints + SprinklingPoints.regrid_like + SprinklingPoints.clip_box + SprinklingPoints.from_imod5_data + SprinklingPoints.get_regrid_methods + SprinklingPoints.write Initial conditions ================== From a7b1d6dee904839f055dd4343d5a8cd0773224df Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Thu, 13 Aug 2026 13:49:52 +0200 Subject: [PATCH 19/38] Migrate from_imod5_ipf logic into designated function, update docstring, and throw TypeError when format not supported by class --- imod/msw/sprinkling.py | 153 ++++++++++++++++--------- imod/tests/test_msw/test_sprinkling.py | 7 +- 2 files changed, 105 insertions(+), 55 deletions(-) diff --git a/imod/msw/sprinkling.py b/imod/msw/sprinkling.py index d4b6f6056..698240378 100644 --- a/imod/msw/sprinkling.py +++ b/imod/msw/sprinkling.py @@ -1,3 +1,4 @@ +import textwrap from typing import TextIO import numpy as np @@ -29,8 +30,29 @@ def _ravel_per_subunit(da: xr.DataArray) -> np.ndarray: def _sprinkling_data_from_imod5_ipf(cap_data: GridDataDict) -> GridDataDict: - raise NotImplementedError( - "Assigning sprinkling wells with an IPF file is not supported, please specify them as IDF." + art_grid = cap_data["artificial_recharge"] + df_points = cap_data["artificial_recharge_layer"] + + arl_points = df_points.iloc[:, :5] + arl_points.columns = ["x_p", "y_p", "layer_p", "id2grid_p", "capacity_p"] + # Enforce dtypes + arl_points = arl_points.astype( + { + "x_p": float, + "y_p": float, + "layer_p": int, + "id2grid_p": int, + "capacity_p": float, + } + ) + + return dict( + art_grid=art_grid, + x_p=arl_points["x_p"].to_numpy(), + y_p=arl_points["y_p"].to_numpy(), + layer_p=arl_points["layer_p"].to_numpy(), + id2grid_p=arl_points["id2grid_p"].to_numpy(), + capacity_p=arl_points["capacity_p"].to_numpy(), ) @@ -351,32 +373,28 @@ def _render( @classmethod def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "Sprinkling": """ - Import sprinkling data from imod5 data. Abstraction data for sprinkling - is defined in iMOD5 either with grids (IDF) or points (IPF) combined - with a grid. Depending on the type, the method does different conversions: - - - grids (IDF) - The ``"artifical_recharge_layer"`` variable was defined as grid - (IDF), this grid defines in which layer a groundwater abstraction - well should be placed. The ``"artificial_recharge"`` grid contains - types which point to the type of abstraction: - - * 0: no abstraction - * 1: groundwater abstraction - * 2: surfacewater abstraction - - The ``"artificial_recharge_capacity"`` grid/constant defines the - capacity of each groundwater or surfacewater abstraction. This is an - ``1:1`` mapping: Each grid cell maps to a separate well. - - - points with grid (IPF & IDF) - The ``"artifical_recharge_layer"`` variable was defined as point - data (IPF), this table contains wellids with an abstraction capacity - and layer. The ``"artificial_recharge"`` grid contains a mapping of - grid cells to wellids in the point data. The - ``"artificial_recharge_capacity"`` is ignored as the abstraction - capacity is already defined in the point data. This is an ``n:1`` - mapping: multiple grid cells can map to one well. + Import sprinkling data from imod5 data artificial recharge grids. + Abstraction data for sprinkling is defined in iMOD5 either with grids + (IDF) or points (IPF) combined with a grid. This class can handle only + the purely grid (IDF) variant. For point data (IPF), use + :class:`imod.msw.SprinklingPoints.from_imod5_data()` instead. + + The iMOD5 data is expected to contain three grids for sprinkling: + + 1. The ``"artificial_recharge"`` grid contains types which point to the + type of abstraction: + + * **0**: no abstraction + * **1**: groundwater abstraction + * **2**: surfacewater abstraction + + 2. The ``"artificial_recharge_layer"`` defines in which layer a groundwater + abstraction well should be placed. + 3. The ``"artificial_recharge_capacity"`` grid/constant defines the + capacity of each groundwater or surfacewater abstraction. This is + converted from mm/d to m3/d using the cell area of the SVAT grid. + + This is an ``1:1`` mapping: Each grid cell maps to a separate well. Parameters ---------- @@ -392,9 +410,16 @@ def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "Sprinkling": """ cap_data = imod5_data["cap"] if isinstance(cap_data["artificial_recharge_layer"], pd.DataFrame): - data = _sprinkling_data_from_imod5_ipf(cap_data) - else: - data = _sprinkling_data_from_imod5_grid(cap_data) + msg = textwrap.dedent( + """ + Unsupported format for artificial_recharge_layer: expected a + grid (IDF) got a DataFrame for point data (IPF). Call + imod.msw.SprinklingPoints.from_imod5_data() instead. + """ + ) + raise TypeError(msg) + + data = _sprinkling_data_from_imod5_grid(cap_data) return cls(**data) @@ -480,31 +505,51 @@ def __init__( @classmethod def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "SprinklingPoints": - cap_data = imod5_data["cap"] - art_grid = cap_data["artificial_recharge"] - df_points = cap_data["artificial_recharge_layer"] + """ + Import sprinkling data from imod5 data artificial recharge grids. + Abstraction data for sprinkling is defined in iMOD5 either with grids + (IDF) or points (IPF) combined with a grid. This class can handle only + the point (IPF) variant. For grid data (IDF), use + :class:`imod.msw.Sprinkling.from_imod5_data()` instead. + + The iMOD5 data is expected to contain one grid (IDF) and one table with + points for sprinkling (IPF): + + 1. The ``"artificial_recharge"`` grid contains a mapping of + grid cells to wellids in the point data. + 2. The ``"artificial_recharge_layer"`` variable was defined as point + data (IPF), this table contains wellids with an abstraction capacity + and layer. + 3. The ``"artificial_recharge_capacity"`` is ignored as the abstraction + capacity is already defined in the point data. + + This is an ``n:1`` mapping: multiple grid cells can map to one well. - arl_points = df_points.iloc[:, :5] - arl_points.columns = ["x_p", "y_p", "layer_p", "id2grid_p", "capacity_p"] - # Enforce dtypes - arl_points = arl_points.astype( - { - "x_p": float, - "y_p": float, - "layer_p": int, - "id2grid_p": int, - "capacity_p": float, - } - ) + Parameters + ---------- + imod5_data: dict[str, dict[str, GridDataArray]] + dictionary containing the arrays mentioned in the project file as + xarray datasets, under the key of the package type to which it + belongs, as returned by + :func:`imod.formats.prj.open_projectfile_data`. - return cls( - art_grid=art_grid, - x_p=arl_points["x_p"].to_numpy(), - y_p=arl_points["y_p"].to_numpy(), - layer_p=arl_points["layer_p"].to_numpy(), - id2grid_p=arl_points["id2grid_p"].to_numpy(), - capacity_p=arl_points["capacity_p"].to_numpy(), - ) + Returns + ------- + SprinklingPoints package + """ + cap_data = imod5_data["cap"] + if isinstance(cap_data["artificial_recharge_layer"], pd.DataFrame): + data = _sprinkling_data_from_imod5_ipf(cap_data) + return cls(**data) + else: + msg = textwrap.dedent( + """ + Unsupported format for artificial_recharge_layer: expected a + DataFrame for point data (IPF), got a grid (IDF). Call + imod.msw.Sprinkling.from_imod5_data() instead. + """ + ) + raise TypeError(msg) def _render(self, file, index, svat, mf6_dis, mf6_well): """ diff --git a/imod/tests/test_msw/test_sprinkling.py b/imod/tests/test_msw/test_sprinkling.py index 6e7845cc0..7c039fef3 100644 --- a/imod/tests/test_msw/test_sprinkling.py +++ b/imod/tests/test_msw/test_sprinkling.py @@ -250,9 +250,14 @@ def test_simple_model_1_subunit( @pytest.mark.unittest_jit def test_sprinkling_from_imod5_data__points(cap_data_sprinkling_points): - with pytest.raises(NotImplementedError): + with pytest.raises(TypeError): msw.Sprinkling.from_imod5_data(cap_data_sprinkling_points) +@pytest.mark.unittest_jit +def test_sprinklingpoints_from_imod5_data__grid(cap_data_sprinkling_grid): + with pytest.raises(TypeError): + msw.SprinklingPoints.from_imod5_data(cap_data_sprinkling_grid) + @pytest.mark.unittest_jit def test_sprinkling_from_imod5_data__grid(cap_data_sprinkling_grid): From 2326477e97d9ca211b7c46219f35598a5c5afa5d Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Thu, 13 Aug 2026 14:32:34 +0200 Subject: [PATCH 20/38] Format and fix mypy issues by adding some specific typeddicts --- imod/msw/sprinkling.py | 81 +++++++++++++++++--------- imod/tests/test_msw/test_sprinkling.py | 2 +- 2 files changed, 54 insertions(+), 29 deletions(-) diff --git a/imod/msw/sprinkling.py b/imod/msw/sprinkling.py index 698240378..a394b9f1b 100644 --- a/imod/msw/sprinkling.py +++ b/imod/msw/sprinkling.py @@ -1,5 +1,5 @@ import textwrap -from typing import TextIO +from typing import TextIO, TypedDict, cast import numpy as np import pandas as pd @@ -18,10 +18,30 @@ from imod.msw.utilities.imod5_converter import ( get_cell_area_from_imod5_data, ) -from imod.typing import GridDataDict, Imod5DataDict, IntArray +from imod.typing import GridDataArray, GridDataDict, Imod5DataDict, IntArray from imod.typing.grid import zeros_like +# Some additional type aliases for sprinkling data, which is a bit more complex +# than other packages. +class CapSprinklingDataDict(TypedDict, total=False): + artificial_recharge: GridDataArray + artificial_recharge_layer: pd.DataFrame + artificial_recharge_capacity: GridDataArray + + +class SprinklingPointsDataDict(TypedDict, total=False): + x_p: np.ndarray | list[float] + y_p: np.ndarray | list[float] + layer_p: np.ndarray | list[int] + id2grid_p: np.ndarray | list[int] + capacity_p: np.ndarray | list[float] + + +class SprinklingPointsGridDataDict(SprinklingPointsDataDict, total=False): + art_grid: GridDataArray + + def _ravel_per_subunit(da: xr.DataArray) -> np.ndarray: # per defined well element, all subunits array_out = da.to_numpy().ravel() @@ -29,32 +49,37 @@ def _ravel_per_subunit(da: xr.DataArray) -> np.ndarray: return array_out[np.isfinite(array_out)] -def _sprinkling_data_from_imod5_ipf(cap_data: GridDataDict) -> GridDataDict: +def _sprinkling_data_from_imod5_ipf( + cap_data: CapSprinklingDataDict, +) -> SprinklingPointsGridDataDict: art_grid = cap_data["artificial_recharge"] df_points = cap_data["artificial_recharge_layer"] + # Select first 5 columns and enforce column names, iMOD5 expects columns in + # this order. The additional columns are metadata for the user and can be + # ignored. arl_points = df_points.iloc[:, :5] arl_points.columns = ["x_p", "y_p", "layer_p", "id2grid_p", "capacity_p"] # Enforce dtypes - arl_points = arl_points.astype( - { - "x_p": float, - "y_p": float, - "layer_p": int, - "id2grid_p": int, - "capacity_p": float, - } - ) + dtype_dict = { + "x_p": float, + "y_p": float, + "layer_p": int, + "id2grid_p": int, + "capacity_p": float, + } - return dict( - art_grid=art_grid, - x_p=arl_points["x_p"].to_numpy(), - y_p=arl_points["y_p"].to_numpy(), - layer_p=arl_points["layer_p"].to_numpy(), - id2grid_p=arl_points["id2grid_p"].to_numpy(), - capacity_p=arl_points["capacity_p"].to_numpy(), + arl_points = arl_points.astype(dtype_dict) + arl_point_dict = cast( + SprinklingPointsDataDict, + {key: arl_points[key].to_numpy() for key in dtype_dict.keys()}, ) + return { + "art_grid": art_grid, + **arl_point_dict, + } + def _sprinkling_data_from_imod5_grid(cap_data: GridDataDict) -> GridDataDict: # Convert units from mm/d to m3/d @@ -378,7 +403,7 @@ def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "Sprinkling": (IDF) or points (IPF) combined with a grid. This class can handle only the purely grid (IDF) variant. For point data (IPF), use :class:`imod.msw.SprinklingPoints.from_imod5_data()` instead. - + The iMOD5 data is expected to contain three grids for sprinkling: 1. The ``"artificial_recharge"`` grid contains types which point to the @@ -387,13 +412,13 @@ def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "Sprinkling": * **0**: no abstraction * **1**: groundwater abstraction * **2**: surfacewater abstraction - + 2. The ``"artificial_recharge_layer"`` defines in which layer a groundwater - abstraction well should be placed. + abstraction well should be placed. 3. The ``"artificial_recharge_capacity"`` grid/constant defines the capacity of each groundwater or surfacewater abstraction. This is converted from mm/d to m3/d using the cell area of the SVAT grid. - + This is an ``1:1`` mapping: Each grid cell maps to a separate well. Parameters @@ -418,7 +443,7 @@ def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "Sprinkling": """ ) raise TypeError(msg) - + data = _sprinkling_data_from_imod5_grid(cap_data) return cls(**data) @@ -514,15 +539,15 @@ def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "SprinklingPoints": The iMOD5 data is expected to contain one grid (IDF) and one table with points for sprinkling (IPF): - + 1. The ``"artificial_recharge"`` grid contains a mapping of grid cells to wellids in the point data. 2. The ``"artificial_recharge_layer"`` variable was defined as point data (IPF), this table contains wellids with an abstraction capacity and layer. 3. The ``"artificial_recharge_capacity"`` is ignored as the abstraction - capacity is already defined in the point data. - + capacity is already defined in the point data. + This is an ``n:1`` mapping: multiple grid cells can map to one well. Parameters @@ -537,7 +562,7 @@ def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "SprinklingPoints": ------- SprinklingPoints package """ - cap_data = imod5_data["cap"] + cap_data = cast(CapSprinklingDataDict, imod5_data["cap"]) if isinstance(cap_data["artificial_recharge_layer"], pd.DataFrame): data = _sprinkling_data_from_imod5_ipf(cap_data) return cls(**data) diff --git a/imod/tests/test_msw/test_sprinkling.py b/imod/tests/test_msw/test_sprinkling.py index 7c039fef3..765514f7c 100644 --- a/imod/tests/test_msw/test_sprinkling.py +++ b/imod/tests/test_msw/test_sprinkling.py @@ -247,12 +247,12 @@ def test_simple_model_1_subunit( assert_equal(results["svat_groundwater"], case_data.expected_svat_gw[:2]) - @pytest.mark.unittest_jit def test_sprinkling_from_imod5_data__points(cap_data_sprinkling_points): with pytest.raises(TypeError): msw.Sprinkling.from_imod5_data(cap_data_sprinkling_points) + @pytest.mark.unittest_jit def test_sprinklingpoints_from_imod5_data__grid(cap_data_sprinkling_grid): with pytest.raises(TypeError): From 709038cae1188635cbd758cab9531b6966dbe4ac Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Thu, 13 Aug 2026 17:20:53 +0200 Subject: [PATCH 21/38] Finalize renaming to "capacity_p", deal with extra edge case if svat_groundwater = 0, and properly order the columns of the dataframe --- imod/msw/sprinkling.py | 37 +++++++++++++++++++++++++++++-------- 1 file changed, 29 insertions(+), 8 deletions(-) diff --git a/imod/msw/sprinkling.py b/imod/msw/sprinkling.py index a394b9f1b..4731fda5f 100644 --- a/imod/msw/sprinkling.py +++ b/imod/msw/sprinkling.py @@ -265,7 +265,7 @@ def _get_wells_outside_art_grid_dataframe( checking which wells in mf6_cellid_df are not in msw_mf6_merged_df. """ is_outside_art = ~mf6_cellid_df["id"].isin(msw_mf6_sprinkling_df["id"].unique()) - outside_df = points_df.loc[is_outside_art, ["layer", "capacity"]] + outside_df = points_df.loc[is_outside_art, ["layer", "capacity_p"]] outside_df = _replicate_dataframe_by_subunit(outside_df) # Select the SVAT subunit for these wells based on their row/col location. cellid_outside_df = mf6_cellid_df.loc[is_outside_art] @@ -277,7 +277,9 @@ def _get_wells_outside_art_grid_dataframe( outside_df["svat_groundwater"] = svat_outside.astype(int) outside_df["svat"] = svat_outside.astype(int) # Set capacity to surfacewater abstraction, and set groundwater abstraction to 0. - outside_df = outside_df.rename(columns={"capacity": "max_abstraction_surfacewater"}) + outside_df = outside_df.rename( + columns={"capacity_p": "max_abstraction_surfacewater"} + ) outside_df["max_abstraction_groundwater"] = 0.0 # drop subunit column as it is no longer needed outside_df = outside_df.drop(columns=["subunit"]) @@ -604,7 +606,7 @@ def _render(self, file, index, svat, mf6_dis, mf6_well): # Select columns that need to be written to scap_svat.inp inside_df = msw_mf6_sprinkling_df[["svat", "layer", "svat_groundwater"]] inside_df["svat"] = inside_df["svat"].astype(int) - capacity = msw_mf6_sprinkling_df["capacity"] + capacity = msw_mf6_sprinkling_df["capacity_p"] # Set wells with layer > 0 to groundwater abstraction, and wells with layer = 0 # to surfacewater abstraction. is_gw_extraction = msw_mf6_sprinkling_df["layer"] > 0 @@ -612,21 +614,40 @@ def _render(self, file, index, svat, mf6_dis, mf6_well): inside_df["max_abstraction_surfacewater"] = capacity.where( ~is_gw_extraction, 0.0 ) + ############## + # EDGE CASES # + ############## + # Set wells located in inactive SVAT groundwater units (svat = 0) to + # surface water extraction. + well_in_inactive_cell = inside_df["svat_groundwater"] == 0 + inside_df.loc[well_in_inactive_cell, "max_abstraction_groundwater"] = 0.0 + inside_df.loc[well_in_inactive_cell, "max_abstraction_surfacewater"] = ( + capacity.where(well_in_inactive_cell, 0.0) + ) + inside_df.loc[well_in_inactive_cell, "svat_groundwater"] = inside_df.loc[ + well_in_inactive_cell, "svat" + ] - # TODO: Make sure wells in svats are all present in the inside_df. If - # these svats are 0 in the art_grid, they should get a 0.0 capacity. - - # Deal with edge case: wells that are outside art_grid, but in model domain. + # Wells that are outside art_grid, but in model domain. # These will be assigned to surfacewater abstraction. outside_df = _get_wells_outside_art_grid_dataframe( mf6_cellid_df, points_df, msw_mf6_sprinkling_df, svat_aligned ) + ############ + # FINALIZE # + ############ # Prepare the final dataframe to be written to scap_svat.inp dataframe_out = pd.concat([inside_df, outside_df], axis=0, ignore_index=True) + # Order rows by SVAT number to ensure consistent output for testing and + # debugging. dataframe_out = dataframe_out.sort_values(by=["svat"]).reset_index(drop=True) - + # Fill last columns with empty strings, as they are not used in the + # iMOD5 implementation but required by MetaSWAP. for var in self._to_fill: dataframe_out[var] = "" + # Order columns to match the metadata dict, which defines the order of + # columns in scap_svat.inp. + dataframe_out = dataframe_out[list(self._metadata_dict.keys())] self._check_range(dataframe_out) From 900de95d8463b2d402bfc5192bebd06b3af836a1 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Thu, 13 Aug 2026 17:21:10 +0200 Subject: [PATCH 22/38] Slightly better name --- imod/tests/fixtures/imod5_cap_data.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/imod/tests/fixtures/imod5_cap_data.py b/imod/tests/fixtures/imod5_cap_data.py index 38c8e814d..b670abdc1 100644 --- a/imod/tests/fixtures/imod5_cap_data.py +++ b/imod/tests/fixtures/imod5_cap_data.py @@ -153,13 +153,13 @@ def cap_data_sprinkling_points() -> Imod5DataDict: "capacity": [15.0, 30.0], } - layer = pd.DataFrame(data=data) + dataframe = pd.DataFrame(data=data) cap_data = { "boundary": boundary, "wetted_area": wetted_area, "urban_area": urban_area, "artificial_recharge": artificial_rch_type, - "artificial_recharge_layer": layer, + "artificial_recharge_layer": dataframe, "artificial_recharge_capacity": xr.DataArray(25.0), } From ab665d2b22bfb328dab82cdced3361ead2d205ae Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Thu, 13 Aug 2026 17:21:56 +0200 Subject: [PATCH 23/38] Add test to write from_imod5_data --- imod/tests/test_msw/test_sprinkling.py | 80 +++++++++++++++++++------- 1 file changed, 60 insertions(+), 20 deletions(-) diff --git a/imod/tests/test_msw/test_sprinkling.py b/imod/tests/test_msw/test_sprinkling.py index 765514f7c..4748864a3 100644 --- a/imod/tests/test_msw/test_sprinkling.py +++ b/imod/tests/test_msw/test_sprinkling.py @@ -21,7 +21,7 @@ def sprinkling_svat_index(): y = [1.0, 2.0, 3.0] subunit = [0, 1] dx = 1.0 - dy = 1.0 + dy = -1.0 # fmt: off svat = xr.DataArray( np.array( @@ -36,7 +36,8 @@ def sprinkling_svat_index(): ] ), dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} + coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy}, + name="svat", ) # fmt: on index = (svat != 0).values.ravel() @@ -306,21 +307,60 @@ def test_sprinklingpoints_from_imod5_data__points(cap_data_sprinkling_points): assert set(sprinkling.dataset.keys()) == expected_vars # No unit conversion is done in SprinklingPoints, as the capacity is already # in m3/d - np.testing.assert_almost_equal(sprinkling.dataset["capacity_p"].max(), 100.0) - - -# @pytest.mark.unittest_jit -# def test_sprinklingpoints_write__points(cap_data_sprinkling_points, tmp_path): -# well_x = cap_data_sprinkling_points["cap"]["x"].values -# well_y = cap_data_sprinkling_points["cap"]["y"].values -# well_layer = cap_data_sprinkling_points["cap"]["layer"].values -# -# # Arrange -# # cellids = derive_cellid_from_points(svat, well_x, well_y, well_layer) -# # well = Mf6Wel(cellids, well_rate) -# -# # Act -# sprinkling = msw.SprinklingPoints.from_imod5_data(cap_data_sprinkling_points) -# -# # sprinkling.write() -# + np.testing.assert_almost_equal(sprinkling.dataset["capacity_p"], [15.0, 30.0]) + + +# TODO: Create more test cases for SprinklingPoints.write() to test edge cases, +# such as wells in inactive SVATs, wells outside the model domain, etc. + + +@pytest.mark.unittest_jit +def test_sprinklingpoints_from_imod5_data_write__points( + sprinkling_svat_index, + fixed_format_parser, + cap_data_sprinkling_points, + cap_coupled_dis_grid, + tmp_path, +): + """ + Test with two wells: one inside the active SVAT area, and one outside the + active SVAT area but still in the model domain. Well nr 2. is not assigned + to anything. Well nr 1. is assigned and is located in the centre cell. In + subunit 1 this cell is inactive and the svats coupled to this well are + assigned to surface water, in subunit 2 this cell is active and this well is + coupled to the groundwater svat. + """ + # Arrange + svat, index = sprinkling_svat_index + df = cap_data_sprinkling_points["cap"]["artificial_recharge_layer"] + well_x = df["x"].to_numpy() + well_y = df["y"].to_numpy() + well_layer = df["layer"].to_numpy() + well_rate = xr.DataArray([0.0, 0.0], dims=("ncellid",)) + well_id = xr.DataArray(["0", "1"], dims=("ncellid",)) + cellids = derive_cellid_from_points(svat, well_x, well_y, well_layer) + mf6_well = Mf6Wel(cellids, well_rate, well_id) + mf6_dis = cap_coupled_dis_grid + directory = tmp_path / "sprinkling_points" + directory.mkdir(parents=True, exist_ok=True) + + # Act + sprinkling = msw.SprinklingPoints.from_imod5_data(cap_data_sprinkling_points) + sprinkling.write(directory, index, svat, mf6_dis, mf6_well) + + results = fixed_format_parser( + directory / msw.Sprinkling._file_name, + msw.Sprinkling._metadata_dict, + ) + + # Assert + # TODO: Check with Hendrik whether this is the appropriate behaviour. + np.testing.assert_equal(results["svat"], [1, 2, 3, 4]) + np.testing.assert_equal(results["svat_groundwater"], [1, 2, 4, 4]) + np.testing.assert_equal(results["layer"], [2, 2, 2, 2]) + np.testing.assert_equal( + results["max_abstraction_surfacewater"], [15.0, 15.0, 0.0, 0.0] + ) + np.testing.assert_equal( + results["max_abstraction_groundwater"], [0.0, 0.0, 15.0, 15.0] + ) From b7fce7830fde890558529b635b1b32ff318bb0ca Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Mon, 17 Aug 2026 10:17:10 +0200 Subject: [PATCH 24/38] Separate expected case data into separate dataclass --- imod/tests/test_msw/test_sprinkling.py | 122 +++++++++++++++---------- 1 file changed, 76 insertions(+), 46 deletions(-) diff --git a/imod/tests/test_msw/test_sprinkling.py b/imod/tests/test_msw/test_sprinkling.py index 4748864a3..2c5635feb 100644 --- a/imod/tests/test_msw/test_sprinkling.py +++ b/imod/tests/test_msw/test_sprinkling.py @@ -43,21 +43,34 @@ def sprinkling_svat_index(): index = (svat != 0).values.ravel() return svat, index +@dataclass +class ExpectedCaseData: + abs_gw: Optional[np.ndarray] = None + abs_sw: Optional[np.ndarray] = None + layer: Optional[np.ndarray] = None + svat: Optional[np.ndarray] = None + svat_gw: Optional[np.ndarray] = None + @dataclass -class AbstractionCaseData: +class SprinklingGridCaseData: max_abstraction_groundwater: Optional[xr.DataArray] = None max_abstraction_surfacewater: Optional[xr.DataArray] = None - expected_abs_gw: Optional[np.ndarray] = None - expected_abs_sw: Optional[np.ndarray] = None - expected_layer: Optional[np.ndarray] = None - expected_svat_gw: Optional[np.ndarray] = None +@dataclass +class SprinklingPointsCaseData: + art_grid: Optional[xr.DataArray] = None + x_p: Optional[xr.DataArray] = None + y_p: Optional[xr.DataArray] = None + layer_p: Optional[xr.DataArray] = None + id2grid_p: Optional[xr.DataArray] = None + capacity_p: Optional[xr.DataArray] = None -class AbstractionGrids: - def case_all_svats(self, sprinkling_svat_index) -> AbstractionCaseData: + +class SprinklingGridCases: + def case_all_svats(self, sprinkling_svat_index) -> tuple[SprinklingGridCaseData, ExpectedCaseData]: svat, _ = sprinkling_svat_index - case_data = AbstractionCaseData() + case_data = SprinklingGridCaseData() case_data.max_abstraction_groundwater = xr.full_like(svat, 0.0) case_data.max_abstraction_surfacewater = xr.full_like(svat, 0.0) # fmt: off @@ -82,15 +95,18 @@ def case_all_svats(self, sprinkling_svat_index) -> AbstractionCaseData: ] ) # fmt: on - case_data.expected_abs_gw = np.array([100.0, 300.0, 100.0, 200.0]) - case_data.expected_abs_sw = np.array([100.0, 300.0, 100.0, 200.0]) - case_data.expected_layer = np.array([3, 1, 3, 2]) - case_data.expected_svat_gw = np.array([1, 2, 3, 4]) - return case_data - - def case_some_svats(self, sprinkling_svat_index): + expected_data = ExpectedCaseData() + expected_data.abs_gw = np.array([100.0, 300.0, 100.0, 200.0]) + expected_data.abs_sw = np.array([100.0, 300.0, 100.0, 200.0]) + expected_data.layer = np.array([3, 1, 3, 2]) + expected_data.svat = np.array([1, 2, 3, 4]) + expected_data.svat_gw = np.array([1, 2, 3, 4]) + case_data.expected_data = expected_data + return case_data, expected_data + + def case_some_svats(self, sprinkling_svat_index) -> tuple[SprinklingGridCaseData, ExpectedCaseData]: svat, _ = sprinkling_svat_index - case_data = AbstractionCaseData() + case_data = SprinklingGridCaseData() case_data.max_abstraction_groundwater = xr.full_like(svat, 0.0) case_data.max_abstraction_surfacewater = xr.full_like(svat, 0.0) # fmt: off @@ -115,16 +131,18 @@ def case_some_svats(self, sprinkling_svat_index): ] ) # fmt: on - case_data.expected_abs_gw = np.array([100.0, 300.0, 200.0]) - case_data.expected_abs_sw = np.array([100.0, 300.0, 200.0]) - case_data.expected_layer = np.array([3, 1, 2]) - case_data.expected_svat_gw = np.array([1, 2, 4]) + expected_data = ExpectedCaseData() + expected_data.abs_gw = np.array([100.0, 300.0, 200.0]) + expected_data.abs_sw = np.array([100.0, 300.0, 200.0]) + expected_data.layer = np.array([3, 1, 2]) + expected_data.svat = np.array([1, 2, 4]) + expected_data.svat_gw = np.array([1, 2, 4]) - return case_data + return case_data, expected_data - def case_inconsistent_active_capacity(self, sprinkling_svat_index): + def case_inconsistent_active_capacity(self, sprinkling_svat_index) -> tuple[SprinklingGridCaseData, ExpectedCaseData]: svat, _ = sprinkling_svat_index - case_data = AbstractionCaseData() + case_data = SprinklingGridCaseData() case_data.max_abstraction_groundwater = xr.full_like(svat, 0.0) case_data.max_abstraction_surfacewater = xr.full_like(svat, 0.0) # fmt: off @@ -149,19 +167,29 @@ def case_inconsistent_active_capacity(self, sprinkling_svat_index): ] ) # fmt: on - case_data.expected_abs_gw = np.array([100.0, 0.0, 200.0]) - case_data.expected_abs_sw = np.array([0.0, 300.0, 200.0]) - case_data.expected_layer = np.array([3, 1, 2]) - case_data.expected_svat_gw = np.array([1, 2, 4]) + expected_data = ExpectedCaseData() + expected_data.abs_gw = np.array([100.0, 0.0, 200.0]) + expected_data.abs_sw = np.array([0.0, 300.0, 200.0]) + expected_data.layer = np.array([3, 1, 2]) + expected_data.svat = np.array([1, 2, 4]) + expected_data.svat_gw = np.array([1, 2, 4]) + + return case_data, expected_data + + +class SprinklingPointsCases: + def case_simple(self, sprinkling_svat_index): + svat, _ = sprinkling_svat_index + case_data = SprinklingPointsCaseData() - return case_data -@parametrize_with_cases("case_data", cases=AbstractionGrids) +@parametrize_with_cases("case_data, expected_data", cases=SprinklingGridCases) def test_simple_model( fixed_format_parser: Callable, sprinkling_svat_index: tuple[xr.DataArray, np.ndarray], - case_data: AbstractionCaseData, + case_data: SprinklingGridCaseData, + expected_data: ExpectedCaseData, ): svat, index = sprinkling_svat_index @@ -188,24 +216,25 @@ def test_simple_model( msw.Sprinkling._metadata_dict, ) - assert_equal(results["svat"], case_data.expected_svat_gw) + assert_equal(results["svat"], expected_data.svat) assert_almost_equal( results["max_abstraction_groundwater"], - case_data.expected_abs_gw, + expected_data.abs_gw, ) assert_almost_equal( results["max_abstraction_surfacewater"], - case_data.expected_abs_sw, + expected_data.abs_sw, ) - assert_equal(results["layer"], case_data.expected_layer) - assert_equal(results["svat_groundwater"], case_data.expected_svat_gw) + assert_equal(results["layer"], expected_data.layer) + assert_equal(results["svat_groundwater"], expected_data.svat_gw) -@parametrize_with_cases("case_data", cases=AbstractionGrids) +@parametrize_with_cases("case_data, expected_data", cases=SprinklingGridCases) def test_simple_model_1_subunit( fixed_format_parser: Callable, sprinkling_svat_index: tuple[xr.DataArray, np.ndarray], - case_data: AbstractionCaseData, + case_data: SprinklingGridCaseData, + expected_data: ExpectedCaseData, ): svat, index = sprinkling_svat_index @@ -235,17 +264,17 @@ def test_simple_model_1_subunit( msw.Sprinkling._metadata_dict, ) - assert_equal(results["svat"], case_data.expected_svat_gw[:2]) + assert_equal(results["svat"], expected_data.svat[:2]) assert_almost_equal( results["max_abstraction_groundwater"], - case_data.expected_abs_gw[:2], + expected_data.abs_gw[:2], ) assert_almost_equal( results["max_abstraction_surfacewater"], - case_data.expected_abs_sw[:2], + expected_data.abs_sw[:2], ) - assert_equal(results["layer"], case_data.expected_layer[:2]) - assert_equal(results["svat_groundwater"], case_data.expected_svat_gw[:2]) + assert_equal(results["layer"], expected_data.layer[:2]) + assert_equal(results["svat_groundwater"], expected_data.svat_gw[:2]) @pytest.mark.unittest_jit @@ -294,6 +323,11 @@ def test_sprinkling_from_imod5_data__grid(cap_data_sprinkling_grid): ) +# TODO: Create more test cases for SprinklingPoints.write() to test edge cases, +# such as wells in inactive SVATs, wells outside the model domain, etc. + + + @pytest.mark.unittest_jit def test_sprinklingpoints_from_imod5_data__points(cap_data_sprinkling_points): # Arrange @@ -310,10 +344,6 @@ def test_sprinklingpoints_from_imod5_data__points(cap_data_sprinkling_points): np.testing.assert_almost_equal(sprinkling.dataset["capacity_p"], [15.0, 30.0]) -# TODO: Create more test cases for SprinklingPoints.write() to test edge cases, -# such as wells in inactive SVATs, wells outside the model domain, etc. - - @pytest.mark.unittest_jit def test_sprinklingpoints_from_imod5_data_write__points( sprinkling_svat_index, From 4f39a69c3bed2a46fff473a7200ffd67fbd2b442 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Mon, 17 Aug 2026 17:31:33 +0200 Subject: [PATCH 25/38] Start adding unittest for SprinklingPoints with a bunch of test cases. Not all of them succeed yet. --- imod/tests/test_msw/test_sprinkling.py | 242 +++++++++++++++++++++++-- 1 file changed, 222 insertions(+), 20 deletions(-) diff --git a/imod/tests/test_msw/test_sprinkling.py b/imod/tests/test_msw/test_sprinkling.py index 2c5635feb..b549e6f16 100644 --- a/imod/tests/test_msw/test_sprinkling.py +++ b/imod/tests/test_msw/test_sprinkling.py @@ -13,15 +13,15 @@ from imod import msw from imod.mf6.mf6_wel_adapter import Mf6Wel from imod.mf6.wel import derive_cellid_from_points - +from imod.mf6.dis import StructuredDiscretization @pytest.fixture(scope="function") def sprinkling_svat_index(): x = [1.0, 2.0, 3.0] - y = [1.0, 2.0, 3.0] + y = [3.0, 2.0, 1.0] subunit = [0, 1] dx = 1.0 - dy = -1.0 + dy = 1.0 # fmt: off svat = xr.DataArray( np.array( @@ -45,6 +45,7 @@ def sprinkling_svat_index(): @dataclass class ExpectedCaseData: + xfail: Optional[str] = None abs_gw: Optional[np.ndarray] = None abs_sw: Optional[np.ndarray] = None layer: Optional[np.ndarray] = None @@ -60,11 +61,11 @@ class SprinklingGridCaseData: @dataclass class SprinklingPointsCaseData: art_grid: Optional[xr.DataArray] = None - x_p: Optional[xr.DataArray] = None - y_p: Optional[xr.DataArray] = None - layer_p: Optional[xr.DataArray] = None - id2grid_p: Optional[xr.DataArray] = None - capacity_p: Optional[xr.DataArray] = None + x_p: Optional[np.ndarray] = None + y_p: Optional[np.ndarray] = None + layer_p: Optional[np.ndarray] = None + id2grid_p: Optional[np.ndarray] = None + capacity_p: Optional[np.ndarray] = None class SprinklingGridCases: @@ -178,14 +179,159 @@ def case_inconsistent_active_capacity(self, sprinkling_svat_index) -> tuple[Spri class SprinklingPointsCases: - def case_simple(self, sprinkling_svat_index): + def case_simple(self, sprinkling_svat_index) -> tuple[SprinklingPointsCaseData, ExpectedCaseData]: + """ + Simple test case for sprinkling points. Each point is mapped to one svat. + """ + svat, _ = sprinkling_svat_index + case_data = SprinklingPointsCaseData() + case_data.art_grid = xr.full_like(svat.isel(subunit=0, drop=True), 0, dtype=int) + # fmt: off + case_data.art_grid.data = np.array( + [[0, 1, 0], + [0, 2, 0], + [0, 3, 0]] + ) + # fmt: on + case_data.x_p = [2.0, 2.0, 2.0] + case_data.y_p = [3.0, 2.0, 1.0] + case_data.layer_p = [1, 2, 3] + case_data.id2grid_p = [1, 2, 3] + case_data.capacity_p = [10.0, 20.0, 30.0] + + expected_data = ExpectedCaseData() + expected_data.svat = np.array([1, 2, 3, 4]) + expected_data.svat_gw = np.array([1, 2, 3, 4]) + expected_data.layer = np.array([1, 3, 1, 2]) + expected_data.abs_gw = np.array([10.0, 30.0, 10.0, 20.0]) + expected_data.abs_sw = np.array([0.0, 0.0, 0.0, 0.0]) + + return case_data, expected_data + + + def case_multi_point_one_art_cell(self, sprinkling_svat_index) -> tuple[SprinklingPointsCaseData, ExpectedCaseData]: + """ + Case where multiple points are assigned to the same SVAT. Not a common + usecase. Usually multiple cells coupled to one point. + """ + svat, _ = sprinkling_svat_index + case_data = SprinklingPointsCaseData() + case_data.art_grid = xr.full_like(svat.isel(subunit=0, drop=True), 0, dtype=int) + # fmt: off + case_data.art_grid.data = np.array( + [[0, 0, 0], + [0, 1, 0], + [0, 0, 0]] + ) + # fmt: on + case_data.x_p = [2.0, 2.0, 2.0] + case_data.y_p = [3.0, 2.0, 1.0] + case_data.layer_p = [1, 2, 3] + case_data.id2grid_p = [1, 1, 1] + case_data.capacity_p = [10.0, 20.0, 30.0] + + expected_data = ExpectedCaseData() + expected_data.xfail = "Multiple points cannot be connected to one grid cell" + return case_data, expected_data + + def case_one_point_multi_art_cell(self, sprinkling_svat_index) -> tuple[SprinklingPointsCaseData, ExpectedCaseData]: + """ + Case where one point is assigned to multiple art_grid cells. Quite a + common usecase. The point is located in the centre of the grid, where + there is only an svat in subunit 1. In subunit 0 this cell is not + active, therefore sprinkling capacity is assigned to surface water. + """ + svat, _ = sprinkling_svat_index + case_data = SprinklingPointsCaseData() + case_data.art_grid = xr.full_like(svat.isel(subunit=0, drop=True), 0, dtype=int) + # fmt: off + case_data.art_grid.data = np.array( + [[0, 1, 0], + [0, 1, 0], + [0, 1, 0]] + ) + # fmt: on + case_data.x_p = [2.0] + case_data.y_p = [2.0] + case_data.layer_p = [2] + case_data.id2grid_p = [1] + case_data.capacity_p = [10.0] + + expected_data = ExpectedCaseData() + expected_data.svat = np.array([1, 2, 3, 4]) + expected_data.svat_gw = np.array([1, 2, 4, 4]) + expected_data.layer = np.array([2, 2, 2, 2]) + expected_data.abs_gw = np.array([0.0, 0.0, 10.0, 10.0]) + expected_data.abs_sw = np.array([10.0, 10.0, 0.0, 0.0]) + + return case_data, expected_data + + + def case_art_grid_outside(self, sprinkling_svat_index) -> tuple[SprinklingPointsCaseData, ExpectedCaseData]: + """ + Case where art grid is located outside the active SVAT area, but still in + the model domain. The well should not be assigned. + """ + svat, _ = sprinkling_svat_index + case_data = SprinklingPointsCaseData() + case_data.art_grid = xr.full_like(svat.isel(subunit=0, drop=True), 0, dtype=int) + # fmt: off + case_data.art_grid.data = np.array( + [[0, 1, 4], + [0, 2, 0], + [0, 3, 0]] + ) + # fmt: on + case_data.x_p = [2.0, 2.0, 2.0, 2.0] + case_data.y_p = [3.0, 2.0, 1.0, 3.0] + case_data.layer_p = [1, 2, 3, 3] + case_data.id2grid_p = [1, 2, 3, 4] + case_data.capacity_p = [10.0, 20.0, 30.0, 40.0] + + expected_data = ExpectedCaseData() + expected_data.svat = np.array([1, 2, 3, 4]) + expected_data.svat_gw = np.array([1, 2, 3, 4]) + expected_data.layer = np.array([1, 3, 1, 2]) + expected_data.abs_gw = np.array([10.0, 30.0, 10.0, 20.0]) + expected_data.abs_sw = np.array([0.0, 0.0, 0.0, 0.0]) + + + return case_data, expected_data + + + def case_point_outside(self, sprinkling_svat_index) -> tuple[SprinklingPointsCaseData, ExpectedCaseData]: + """ + Case where one point is located outside the active SVAT area, but still in + the model domain. The well should be assigned as surface water extraction. + """ svat, _ = sprinkling_svat_index case_data = SprinklingPointsCaseData() + case_data.art_grid = xr.full_like(svat.isel(subunit=0, drop=True), 0, dtype=int) + # fmt: off + case_data.art_grid.data = np.array( + [[0, 1, 0], + [0, 2, 0], + [0, 3, 0]] + ) + # fmt: on + case_data.x_p = [2.0, 2.0, 3.0] + case_data.y_p = [1.0, 2.0, 1.0] + case_data.layer_p = [1, 2, 3] + case_data.id2grid_p = [1, 2, 1] + case_data.capacity_p = [10.0, 20.0, 30.0, 40.0] + + expected_data = ExpectedCaseData() + expected_data.svat = np.array([1, 1, 2, 3, 3, 4]) + expected_data.svat_gw = np.array([1, 1, 2, 3, 3, 4]) + expected_data.layer = np.array([1, 1, 2, 1, 1, 3]) + expected_data.abs_gw = np.array([10.0, 0.0, 20.0, 10.0, 0.0, 30.0]) + expected_data.abs_sw = np.array([0.0, 40.0, 0.0, 0.0, 40.0, 0.0]) + return case_data, expected_data @parametrize_with_cases("case_data, expected_data", cases=SprinklingGridCases) -def test_simple_model( +def test_grid_simple_model( fixed_format_parser: Callable, sprinkling_svat_index: tuple[xr.DataArray, np.ndarray], case_data: SprinklingGridCaseData, @@ -197,8 +343,10 @@ def test_simple_model( well_layer = [3, 2, 1] well_y = [1.0, 2.0, 3.0] well_x = [2.0, 2.0, 2.0] - well_rate = [-5.0] * 3 - well_id = ["a", "b", "c"] + well_rate_values = [-5.0] * 3 + well_rate = xr.DataArray(well_rate_values, dims=("ncellid",)) + well_id_values = ["0", "1", "2"] + well_id = xr.DataArray(well_id_values, dims=("ncellid",)) cellids = derive_cellid_from_points(svat, well_x, well_y, well_layer) well = Mf6Wel(cellids, well_rate, well_id) @@ -230,7 +378,7 @@ def test_simple_model( @parametrize_with_cases("case_data, expected_data", cases=SprinklingGridCases) -def test_simple_model_1_subunit( +def test_grid_simple_model_1_subunit( fixed_format_parser: Callable, sprinkling_svat_index: tuple[xr.DataArray, np.ndarray], case_data: SprinklingGridCaseData, @@ -245,8 +393,10 @@ def test_simple_model_1_subunit( well_layer = [3, 1] well_y = [1.0, 3.0] well_x = [2.0, 2.0] - well_rate = [-5.0] * 2 - well_id = ["a", "c"] + well_rate_values = [-5.0] * 2 + well_rate = xr.DataArray(well_rate_values, dims=("ncellid",)) + well_id_values = ["0", "2"] + well_id = xr.DataArray(well_id_values, dims=("ncellid",)) cellids = derive_cellid_from_points(svat, well_x, well_y, well_layer) well = Mf6Wel(cellids, well_rate, well_id) @@ -277,6 +427,63 @@ def test_simple_model_1_subunit( assert_equal(results["svat_groundwater"], expected_data.svat_gw[:2]) +@parametrize_with_cases("case_data, expected_data", cases=SprinklingPointsCases) +def test_points_simple_model( + fixed_format_parser: Callable, + sprinkling_svat_index: tuple[xr.DataArray, np.ndarray], + case_data: SprinklingPointsCaseData, + expected_data: ExpectedCaseData, +): + if expected_data.xfail: + pytest.xfail(expected_data.xfail) + svat, index = sprinkling_svat_index + + # Well + n_wells = len(case_data.x_p) + well_rate_values = [-5.0] * n_wells + well_rate = xr.DataArray(well_rate_values, dims=("ncellid",)) + well_id_values = [str(i) for i in np.arange(n_wells)] + well_id = xr.DataArray(well_id_values, dims=("ncellid",)) + + cellids = derive_cellid_from_points(svat, case_data.x_p, case_data.y_p, case_data.layer_p) + well = Mf6Wel(cellids, well_rate, well_id) + + layer_template = xr.DataArray([1.0, 2.0, 3.0], coords={"layer": [1, 2, 3]}, dims=("layer",)) + grid_2d_template = xr.ones_like(svat.isel(subunit=0, drop=True), dtype=float) + mf6_dis_template = layer_template * grid_2d_template + + dis = StructuredDiscretization(top=grid_2d_template, bottom=-mf6_dis_template, idomain=mf6_dis_template.astype(int)) + + sprinkling = msw.SprinklingPoints( + case_data.art_grid, + case_data.x_p, + case_data.y_p, + case_data.layer_p, + case_data.id2grid_p, + case_data.capacity_p + ) + + with tempfile.TemporaryDirectory() as output_dir: + output_dir = Path(output_dir) + sprinkling.write(output_dir, index, svat, dis, well) + + results = fixed_format_parser( + output_dir / msw.SprinklingPoints._file_name, + msw.SprinklingPoints._metadata_dict, + ) + + assert_equal(results["svat"], expected_data.svat) + assert_almost_equal( + results["max_abstraction_groundwater"], + expected_data.abs_gw, + ) + assert_almost_equal( + results["max_abstraction_surfacewater"], + expected_data.abs_sw, + ) + assert_equal(results["layer"], expected_data.layer) + assert_equal(results["svat_groundwater"], expected_data.svat_gw) + @pytest.mark.unittest_jit def test_sprinkling_from_imod5_data__points(cap_data_sprinkling_points): with pytest.raises(TypeError): @@ -323,11 +530,6 @@ def test_sprinkling_from_imod5_data__grid(cap_data_sprinkling_grid): ) -# TODO: Create more test cases for SprinklingPoints.write() to test edge cases, -# such as wells in inactive SVATs, wells outside the model domain, etc. - - - @pytest.mark.unittest_jit def test_sprinklingpoints_from_imod5_data__points(cap_data_sprinkling_points): # Arrange From 41885609f8a36a146509b4b9bc216c49406fd923 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Tue, 18 Aug 2026 08:11:34 +0200 Subject: [PATCH 26/38] Simplify case where the art grid has a mapping id in an inactive svat --- imod/tests/test_msw/test_sprinkling.py | 33 ++++++++++++++------------ 1 file changed, 18 insertions(+), 15 deletions(-) diff --git a/imod/tests/test_msw/test_sprinkling.py b/imod/tests/test_msw/test_sprinkling.py index b549e6f16..93c26dda8 100644 --- a/imod/tests/test_msw/test_sprinkling.py +++ b/imod/tests/test_msw/test_sprinkling.py @@ -269,31 +269,34 @@ def case_one_point_multi_art_cell(self, sprinkling_svat_index) -> tuple[Sprinkli def case_art_grid_outside(self, sprinkling_svat_index) -> tuple[SprinklingPointsCaseData, ExpectedCaseData]: """ - Case where art grid is located outside the active SVAT area, but still in - the model domain. The well should not be assigned. + Case where art grid is located outside the active SVAT area, but still + in the model domain. The well is inside the model domain. Sprinkling + should not be assigned. """ svat, _ = sprinkling_svat_index case_data = SprinklingPointsCaseData() case_data.art_grid = xr.full_like(svat.isel(subunit=0, drop=True), 0, dtype=int) # fmt: off case_data.art_grid.data = np.array( - [[0, 1, 4], - [0, 2, 0], - [0, 3, 0]] + [[0, 0, 4], + [0, 0, 0], + [0, 0, 0]] ) # fmt: on - case_data.x_p = [2.0, 2.0, 2.0, 2.0] - case_data.y_p = [3.0, 2.0, 1.0, 3.0] - case_data.layer_p = [1, 2, 3, 3] - case_data.id2grid_p = [1, 2, 3, 4] - case_data.capacity_p = [10.0, 20.0, 30.0, 40.0] + case_data.x_p = [2.0] + case_data.y_p = [2.0] + case_data.layer_p = [3] + case_data.id2grid_p = [4] + case_data.capacity_p = [40.0] + # TODO: This is wrong, this should be for the case where a well is + # outside the grid. expected_data = ExpectedCaseData() - expected_data.svat = np.array([1, 2, 3, 4]) - expected_data.svat_gw = np.array([1, 2, 3, 4]) - expected_data.layer = np.array([1, 3, 1, 2]) - expected_data.abs_gw = np.array([10.0, 30.0, 10.0, 20.0]) - expected_data.abs_sw = np.array([0.0, 0.0, 0.0, 0.0]) + expected_data.svat = np.array([4]) + expected_data.svat_gw = np.array([4]) + expected_data.layer = np.array([3]) + expected_data.abs_gw = np.array([0.0]) + expected_data.abs_sw = np.array([40.0]) return case_data, expected_data From a5f01aca26e7c745a8db85fcf2a55d308ae68810 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Tue, 18 Aug 2026 15:56:32 +0200 Subject: [PATCH 27/38] Fix edge case that wasn't working and reduce code a lot --- imod/msw/sprinkling.py | 69 +++++++++++------------------------------- 1 file changed, 17 insertions(+), 52 deletions(-) diff --git a/imod/msw/sprinkling.py b/imod/msw/sprinkling.py index 4731fda5f..4b71c94c5 100644 --- a/imod/msw/sprinkling.py +++ b/imod/msw/sprinkling.py @@ -253,40 +253,6 @@ def _get_svat_groundwater_for_wells( return svat_groundwater.astype(int) -def _get_wells_outside_art_grid_dataframe( - mf6_cellid_df: pd.DataFrame, - points_df: pd.DataFrame, - msw_mf6_sprinkling_df: pd.DataFrame, - svat_aligned: xr.DataArray, -) -> pd.DataFrame: - """ - Deal with edge case: wells that are outside art_grid, but in model domain. - These will be assigned to surfacewater abstraction. We can identify these by - checking which wells in mf6_cellid_df are not in msw_mf6_merged_df. - """ - is_outside_art = ~mf6_cellid_df["id"].isin(msw_mf6_sprinkling_df["id"].unique()) - outside_df = points_df.loc[is_outside_art, ["layer", "capacity_p"]] - outside_df = _replicate_dataframe_by_subunit(outside_df) - # Select the SVAT subunit for these wells based on their row/col location. - cellid_outside_df = mf6_cellid_df.loc[is_outside_art] - cellid_outside_df = _replicate_dataframe_by_subunit(cellid_outside_df) - indexer_outside = _extract_indexer_for_svat( - cellid_outside_df, columns=["subunit", "row", "column"] - ) - svat_outside = svat_aligned.data[*indexer_outside] - outside_df["svat_groundwater"] = svat_outside.astype(int) - outside_df["svat"] = svat_outside.astype(int) - # Set capacity to surfacewater abstraction, and set groundwater abstraction to 0. - outside_df = outside_df.rename( - columns={"capacity_p": "max_abstraction_surfacewater"} - ) - outside_df["max_abstraction_groundwater"] = 0.0 - # drop subunit column as it is no longer needed - outside_df = outside_df.drop(columns=["subunit"]) - # drop wells that are outside the active metaswap model domain (svat = 0) - return outside_df.query("svat > 0").reset_index(drop=True) - - class Sprinkling(MetaSwapPackage, IRegridPackage): """ This contains the sprinkling capacities of links between SVAT units and @@ -602,14 +568,16 @@ def _render(self, file, index, svat, mf6_dis, mf6_well): msw_mf6_sprinkling_df["svat_groundwater"] = _get_svat_groundwater_for_wells( msw_mf6_sprinkling_df, svat_aligned ) - + is_point_inside = msw_mf6_sprinkling_df["svat_groundwater"] > 0 # Select columns that need to be written to scap_svat.inp - inside_df = msw_mf6_sprinkling_df[["svat", "layer", "svat_groundwater"]] + inside_df = msw_mf6_sprinkling_df.loc[ + is_point_inside, ["svat", "layer", "svat_groundwater"] + ] inside_df["svat"] = inside_df["svat"].astype(int) - capacity = msw_mf6_sprinkling_df["capacity_p"] # Set wells with layer > 0 to groundwater abstraction, and wells with layer = 0 # to surfacewater abstraction. - is_gw_extraction = msw_mf6_sprinkling_df["layer"] > 0 + capacity = msw_mf6_sprinkling_df.loc[is_point_inside, "capacity_p"] + is_gw_extraction = inside_df["layer"] > 0 inside_df["max_abstraction_groundwater"] = capacity.where(is_gw_extraction, 0.0) inside_df["max_abstraction_surfacewater"] = capacity.where( ~is_gw_extraction, 0.0 @@ -617,22 +585,19 @@ def _render(self, file, index, svat, mf6_dis, mf6_well): ############## # EDGE CASES # ############## - # Set wells located in inactive SVAT groundwater units (svat = 0) to - # surface water extraction. - well_in_inactive_cell = inside_df["svat_groundwater"] == 0 - inside_df.loc[well_in_inactive_cell, "max_abstraction_groundwater"] = 0.0 - inside_df.loc[well_in_inactive_cell, "max_abstraction_surfacewater"] = ( - capacity.where(well_in_inactive_cell, 0.0) - ) - inside_df.loc[well_in_inactive_cell, "svat_groundwater"] = inside_df.loc[ - well_in_inactive_cell, "svat" - ] - - # Wells that are outside art_grid, but in model domain. + # 1. Wells that are outside art_grid, but in model domain. # These will be assigned to surfacewater abstraction. - outside_df = _get_wells_outside_art_grid_dataframe( - mf6_cellid_df, points_df, msw_mf6_sprinkling_df, svat_aligned + outside_df = msw_mf6_sprinkling_df.loc[ + ~is_point_inside, ["svat", "layer", "svat_groundwater", "capacity_p"] + ] + # Set capacity to surfacewater abstraction, and set groundwater abstraction to 0. + outside_df = outside_df.rename( + columns={"capacity_p": "max_abstraction_surfacewater"} ) + outside_df["max_abstraction_groundwater"] = 0.0 + # Set svat_groundwater to svat, as these wells are outside art_grid and + # will be assigned to surfacewater abstraction. + outside_df["svat_groundwater"] = outside_df["svat"] ############ # FINALIZE # ############ From 06f40fc20989d3031317748e6886314348a7b23b Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Tue, 18 Aug 2026 15:56:44 +0200 Subject: [PATCH 28/38] Finish unittest cases and format --- imod/tests/test_msw/test_sprinkling.py | 97 ++++++++++++++++---------- 1 file changed, 59 insertions(+), 38 deletions(-) diff --git a/imod/tests/test_msw/test_sprinkling.py b/imod/tests/test_msw/test_sprinkling.py index 93c26dda8..e10aa8e9e 100644 --- a/imod/tests/test_msw/test_sprinkling.py +++ b/imod/tests/test_msw/test_sprinkling.py @@ -11,9 +11,10 @@ from pytest_cases import parametrize_with_cases from imod import msw +from imod.mf6.dis import StructuredDiscretization from imod.mf6.mf6_wel_adapter import Mf6Wel from imod.mf6.wel import derive_cellid_from_points -from imod.mf6.dis import StructuredDiscretization + @pytest.fixture(scope="function") def sprinkling_svat_index(): @@ -43,6 +44,7 @@ def sprinkling_svat_index(): index = (svat != 0).values.ravel() return svat, index + @dataclass class ExpectedCaseData: xfail: Optional[str] = None @@ -58,6 +60,7 @@ class SprinklingGridCaseData: max_abstraction_groundwater: Optional[xr.DataArray] = None max_abstraction_surfacewater: Optional[xr.DataArray] = None + @dataclass class SprinklingPointsCaseData: art_grid: Optional[xr.DataArray] = None @@ -69,7 +72,9 @@ class SprinklingPointsCaseData: class SprinklingGridCases: - def case_all_svats(self, sprinkling_svat_index) -> tuple[SprinklingGridCaseData, ExpectedCaseData]: + def case_all_svats( + self, sprinkling_svat_index + ) -> tuple[SprinklingGridCaseData, ExpectedCaseData]: svat, _ = sprinkling_svat_index case_data = SprinklingGridCaseData() case_data.max_abstraction_groundwater = xr.full_like(svat, 0.0) @@ -105,7 +110,9 @@ def case_all_svats(self, sprinkling_svat_index) -> tuple[SprinklingGridCaseData, case_data.expected_data = expected_data return case_data, expected_data - def case_some_svats(self, sprinkling_svat_index) -> tuple[SprinklingGridCaseData, ExpectedCaseData]: + def case_some_svats( + self, sprinkling_svat_index + ) -> tuple[SprinklingGridCaseData, ExpectedCaseData]: svat, _ = sprinkling_svat_index case_data = SprinklingGridCaseData() case_data.max_abstraction_groundwater = xr.full_like(svat, 0.0) @@ -141,7 +148,9 @@ def case_some_svats(self, sprinkling_svat_index) -> tuple[SprinklingGridCaseData return case_data, expected_data - def case_inconsistent_active_capacity(self, sprinkling_svat_index) -> tuple[SprinklingGridCaseData, ExpectedCaseData]: + def case_inconsistent_active_capacity( + self, sprinkling_svat_index + ) -> tuple[SprinklingGridCaseData, ExpectedCaseData]: svat, _ = sprinkling_svat_index case_data = SprinklingGridCaseData() case_data.max_abstraction_groundwater = xr.full_like(svat, 0.0) @@ -179,7 +188,9 @@ def case_inconsistent_active_capacity(self, sprinkling_svat_index) -> tuple[Spri class SprinklingPointsCases: - def case_simple(self, sprinkling_svat_index) -> tuple[SprinklingPointsCaseData, ExpectedCaseData]: + def case_one_point_one_art_cell( + self, sprinkling_svat_index + ) -> tuple[SprinklingPointsCaseData, ExpectedCaseData]: """ Simple test case for sprinkling points. Each point is mapped to one svat. """ @@ -208,8 +219,9 @@ def case_simple(self, sprinkling_svat_index) -> tuple[SprinklingPointsCaseData, return case_data, expected_data - - def case_multi_point_one_art_cell(self, sprinkling_svat_index) -> tuple[SprinklingPointsCaseData, ExpectedCaseData]: + def case_multi_point_one_art_cell( + self, sprinkling_svat_index + ) -> tuple[SprinklingPointsCaseData, ExpectedCaseData]: """ Case where multiple points are assigned to the same SVAT. Not a common usecase. Usually multiple cells coupled to one point. @@ -234,7 +246,9 @@ def case_multi_point_one_art_cell(self, sprinkling_svat_index) -> tuple[Sprinkli expected_data.xfail = "Multiple points cannot be connected to one grid cell" return case_data, expected_data - def case_one_point_multi_art_cell(self, sprinkling_svat_index) -> tuple[SprinklingPointsCaseData, ExpectedCaseData]: + def case_one_point_multi_art_cell( + self, sprinkling_svat_index + ) -> tuple[SprinklingPointsCaseData, ExpectedCaseData]: """ Case where one point is assigned to multiple art_grid cells. Quite a common usecase. The point is located in the centre of the grid, where @@ -266,8 +280,9 @@ def case_one_point_multi_art_cell(self, sprinkling_svat_index) -> tuple[Sprinkli return case_data, expected_data - - def case_art_grid_outside(self, sprinkling_svat_index) -> tuple[SprinklingPointsCaseData, ExpectedCaseData]: + def case_art_grid_outside( + self, sprinkling_svat_index + ) -> tuple[SprinklingPointsCaseData, ExpectedCaseData]: """ Case where art grid is located outside the active SVAT area, but still in the model domain. The well is inside the model domain. Sprinkling @@ -289,20 +304,18 @@ def case_art_grid_outside(self, sprinkling_svat_index) -> tuple[SprinklingPoints case_data.id2grid_p = [4] case_data.capacity_p = [40.0] - # TODO: This is wrong, this should be for the case where a well is - # outside the grid. expected_data = ExpectedCaseData() - expected_data.svat = np.array([4]) - expected_data.svat_gw = np.array([4]) - expected_data.layer = np.array([3]) - expected_data.abs_gw = np.array([0.0]) - expected_data.abs_sw = np.array([40.0]) - + expected_data.svat = np.array([]) + expected_data.svat_gw = np.array([]) + expected_data.layer = np.array([]) + expected_data.abs_gw = np.array([]) + expected_data.abs_sw = np.array([]) return case_data, expected_data - - def case_point_outside(self, sprinkling_svat_index) -> tuple[SprinklingPointsCaseData, ExpectedCaseData]: + def case_point_outside( + self, sprinkling_svat_index + ) -> tuple[SprinklingPointsCaseData, ExpectedCaseData]: """ Case where one point is located outside the active SVAT area, but still in the model domain. The well should be assigned as surface water extraction. @@ -312,27 +325,26 @@ def case_point_outside(self, sprinkling_svat_index) -> tuple[SprinklingPointsCas case_data.art_grid = xr.full_like(svat.isel(subunit=0, drop=True), 0, dtype=int) # fmt: off case_data.art_grid.data = np.array( - [[0, 1, 0], - [0, 2, 0], - [0, 3, 0]] + [[0, 0, 0], + [0, 0, 0], + [0, 5, 0]] ) # fmt: on - case_data.x_p = [2.0, 2.0, 3.0] - case_data.y_p = [1.0, 2.0, 1.0] - case_data.layer_p = [1, 2, 3] - case_data.id2grid_p = [1, 2, 1] - case_data.capacity_p = [10.0, 20.0, 30.0, 40.0] + case_data.x_p = [3.0] + case_data.y_p = [1.0] + case_data.layer_p = [3] + case_data.id2grid_p = [5] + case_data.capacity_p = [40.0] expected_data = ExpectedCaseData() - expected_data.svat = np.array([1, 1, 2, 3, 3, 4]) - expected_data.svat_gw = np.array([1, 1, 2, 3, 3, 4]) - expected_data.layer = np.array([1, 1, 2, 1, 1, 3]) - expected_data.abs_gw = np.array([10.0, 0.0, 20.0, 10.0, 0.0, 30.0]) - expected_data.abs_sw = np.array([0.0, 40.0, 0.0, 0.0, 40.0, 0.0]) + expected_data.svat = np.array([2]) + expected_data.svat_gw = np.array([2]) + expected_data.layer = np.array([3]) + expected_data.abs_gw = np.array([0.0]) + expected_data.abs_sw = np.array([40.0]) return case_data, expected_data - @parametrize_with_cases("case_data, expected_data", cases=SprinklingGridCases) def test_grid_simple_model( fixed_format_parser: Callable, @@ -448,14 +460,22 @@ def test_points_simple_model( well_id_values = [str(i) for i in np.arange(n_wells)] well_id = xr.DataArray(well_id_values, dims=("ncellid",)) - cellids = derive_cellid_from_points(svat, case_data.x_p, case_data.y_p, case_data.layer_p) + cellids = derive_cellid_from_points( + svat, case_data.x_p, case_data.y_p, case_data.layer_p + ) well = Mf6Wel(cellids, well_rate, well_id) - layer_template = xr.DataArray([1.0, 2.0, 3.0], coords={"layer": [1, 2, 3]}, dims=("layer",)) + layer_template = xr.DataArray( + [1.0, 2.0, 3.0], coords={"layer": [1, 2, 3]}, dims=("layer",) + ) grid_2d_template = xr.ones_like(svat.isel(subunit=0, drop=True), dtype=float) mf6_dis_template = layer_template * grid_2d_template - dis = StructuredDiscretization(top=grid_2d_template, bottom=-mf6_dis_template, idomain=mf6_dis_template.astype(int)) + dis = StructuredDiscretization( + top=grid_2d_template, + bottom=-mf6_dis_template, + idomain=mf6_dis_template.astype(int), + ) sprinkling = msw.SprinklingPoints( case_data.art_grid, @@ -463,7 +483,7 @@ def test_points_simple_model( case_data.y_p, case_data.layer_p, case_data.id2grid_p, - case_data.capacity_p + case_data.capacity_p, ) with tempfile.TemporaryDirectory() as output_dir: @@ -487,6 +507,7 @@ def test_points_simple_model( assert_equal(results["layer"], expected_data.layer) assert_equal(results["svat_groundwater"], expected_data.svat_gw) + @pytest.mark.unittest_jit def test_sprinkling_from_imod5_data__points(cap_data_sprinkling_points): with pytest.raises(TypeError): From 919d438aeb771ac299eb5b1d19e0b9e765b415dc Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 19 Aug 2026 11:09:45 +0200 Subject: [PATCH 29/38] Ensure y-coords are properly oriented --- imod/tests/test_msw/test_sprinkling.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/imod/tests/test_msw/test_sprinkling.py b/imod/tests/test_msw/test_sprinkling.py index e10aa8e9e..1e1e229df 100644 --- a/imod/tests/test_msw/test_sprinkling.py +++ b/imod/tests/test_msw/test_sprinkling.py @@ -22,7 +22,7 @@ def sprinkling_svat_index(): y = [3.0, 2.0, 1.0] subunit = [0, 1] dx = 1.0 - dy = 1.0 + dy = -1.0 # fmt: off svat = xr.DataArray( np.array( @@ -356,7 +356,7 @@ def test_grid_simple_model( # Well well_layer = [3, 2, 1] - well_y = [1.0, 2.0, 3.0] + well_y = [3.0, 2.0, 1.0] well_x = [2.0, 2.0, 2.0] well_rate_values = [-5.0] * 3 well_rate = xr.DataArray(well_rate_values, dims=("ncellid",)) @@ -406,7 +406,7 @@ def test_grid_simple_model_1_subunit( # Well well_layer = [3, 1] - well_y = [1.0, 3.0] + well_y = [3.0, 1.0] well_x = [2.0, 2.0] well_rate_values = [-5.0] * 2 well_rate = xr.DataArray(well_rate_values, dims=("ncellid",)) @@ -610,7 +610,6 @@ def test_sprinklingpoints_from_imod5_data_write__points( ) # Assert - # TODO: Check with Hendrik whether this is the appropriate behaviour. np.testing.assert_equal(results["svat"], [1, 2, 3, 4]) np.testing.assert_equal(results["svat_groundwater"], [1, 2, 4, 4]) np.testing.assert_equal(results["layer"], [2, 2, 2, 2]) From 07353cafe144ed3b7c6a522df7d6455c8def7475 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 19 Aug 2026 11:25:14 +0200 Subject: [PATCH 30/38] Clearer name for the mapping id --- imod/msw/sprinkling.py | 29 +++++++++++++------------- imod/tests/_scratch.py | 14 ++++++------- imod/tests/test_msw/test_sprinkling.py | 16 +++++++------- 3 files changed, 29 insertions(+), 30 deletions(-) diff --git a/imod/msw/sprinkling.py b/imod/msw/sprinkling.py index 4b71c94c5..1f8bf1000 100644 --- a/imod/msw/sprinkling.py +++ b/imod/msw/sprinkling.py @@ -34,7 +34,7 @@ class SprinklingPointsDataDict(TypedDict, total=False): x_p: np.ndarray | list[float] y_p: np.ndarray | list[float] layer_p: np.ndarray | list[int] - id2grid_p: np.ndarray | list[int] + id_sprinkling_p: np.ndarray | list[int] capacity_p: np.ndarray | list[float] @@ -59,13 +59,13 @@ def _sprinkling_data_from_imod5_ipf( # this order. The additional columns are metadata for the user and can be # ignored. arl_points = df_points.iloc[:, :5] - arl_points.columns = ["x_p", "y_p", "layer_p", "id2grid_p", "capacity_p"] + arl_points.columns = ["x_p", "y_p", "layer_p", "id_sprinkling_p", "capacity_p"] # Enforce dtypes dtype_dict = { "x_p": float, "y_p": float, "layer_p": int, - "id2grid_p": int, + "id_sprinkling_p": int, "capacity_p": float, } @@ -208,10 +208,9 @@ def _merge_sprinkling_points_with_grids( Merge sprinkling points with SVAT grids. """ - # TODO: Rename "id_msw" and "id2grid_p" to something clearer like "id_sprinkling" - # Flatten id_msw grid → (y, x, id_msw) table, drop cells with no well + # Flatten id_sprinkling grid → (y, x, id_sprinkling) table, drop cells with no well grids = xr.merge([sprinkling_id_grid, svat]) - art_df = grids.to_dataframe().reset_index().query("(id_msw > 0) & (svat > 0)") + art_df = grids.to_dataframe().reset_index().query("(id_sprinkling > 0) & (svat > 0)") # Drop unnecessary columns. We preserve the x, y coords as they might # prove useful for debugging. art_df = art_df.drop(["dx", "dy"], axis=1) @@ -219,8 +218,8 @@ def _merge_sprinkling_points_with_grids( # Join: each SVAT cell gets the matching well row(s) from arl_points return art_df.merge( points_df, # brings id back as a column - left_on="id_msw", - right_on="id2grid_p", + left_on="id_sprinkling", + right_on="id_sprinkling_p", how="inner", validate="many_to_one", ) @@ -437,7 +436,7 @@ class SprinklingPoints(MetaSwapPackage, IRegridPackage): y-coordinates of the artificial recharge locations. layer_p: np.ndarray | list[int] layer indices of the artificial recharge locations. - id2grid_p: np.ndarray | list[int] + id_sprinkling_p: np.ndarray | list[int] mapping of the artificial recharge locations to the grid cells. capacity_p: np.ndarray | list[float] abstraction capacities of the artificial recharge locations. @@ -476,7 +475,7 @@ def __init__( x_p: np.ndarray | list[float], y_p: np.ndarray | list[float], layer_p: np.ndarray | list[int], - id2grid_p: np.ndarray | list[int], + id_sprinkling_p: np.ndarray | list[int], capacity_p: np.ndarray | list[float], ): super().__init__() @@ -488,12 +487,12 @@ def __init__( "x_p": (("id",), x_p), "y_p": (("id",), y_p), "layer_p": (("id",), layer_p), - "id2grid_p": (("id",), id2grid_p), + "id_sprinkling_p": (("id",), id_sprinkling_p), "capacity_p": (("id",), capacity_p), }, coords={"id": id_index}, ) - art_grid = art_grid.rename("id_msw") + art_grid = art_grid.rename("id_sprinkling") self.dataset = xr.merge([art_grid, points_ds]) @classmethod @@ -550,7 +549,7 @@ def _render(self, file, index, svat, mf6_dis, mf6_well): This method first merges the sprinkling points with the mf6_well cellids to get the row/col of each well, then merges the sprinkling points with - the svat and id_msw grid. It then selects the columns that need to be + the svat and id_sprinkling grid. It then selects the columns that need to be written to scap_svat.inp and sets wells with layer > 0 to groundwater abstraction, and wells with layer = 0 to surfacewater abstraction. Finally, it deals with edge cases for wells that are outside art_grid @@ -560,9 +559,9 @@ def _render(self, file, index, svat, mf6_dis, mf6_well): # row/col of each well. mf6_cellid_df = _get_mf6_cellid_dataframe(mf6_well) points_df = _make_sprinkling_well_points_dataframe(self.dataset, mf6_cellid_df) - # Merge the sprinkling points with the svat and id_msw grid + # Merge the sprinkling points with the svat and id_sprinkling grid msw_mf6_sprinkling_df = _merge_sprinkling_points_with_grids( - points_df, svat, self.dataset["id_msw"] + points_df, svat, self.dataset["id_sprinkling"] ) svat_aligned = align_svat_with_dis(svat, mf6_dis) msw_mf6_sprinkling_df["svat_groundwater"] = _get_svat_groundwater_for_wells( diff --git a/imod/tests/_scratch.py b/imod/tests/_scratch.py index be0f2dba2..912f6df7f 100644 --- a/imod/tests/_scratch.py +++ b/imod/tests/_scratch.py @@ -96,10 +96,10 @@ def double_length_df_subunit( # %% # In from_imod5_cap_data arl_points = df_points.iloc[:, :5] -arl_points.columns = ["x_p", "y_p", "layer_p", "id2grid_p", "capacity"] +arl_points.columns = ["x_p", "y_p", "layer_p", "id_sprinkling_p", "capacity"] # Enforce dtypes arl_points = arl_points.astype( - {"x_p": float, "y_p": float, "layer_p": int, "id2grid_p": int, "capacity": float} + {"x_p": float, "y_p": float, "layer_p": int, "id_sprinkling_p": int, "capacity": float} ) arl_points["id"] = arl_points.index.astype(str) arl_points = arl_points.set_index("id") @@ -107,7 +107,7 @@ def double_length_df_subunit( points_ds = arl_points.to_xarray() # in def __init__ -art_grid = art_grid.rename("id_msw") +art_grid = art_grid.rename("id_sprinkling") dataset = xr.merge([art_grid, points_ds]) # %% @@ -126,9 +126,9 @@ def double_length_df_subunit( mf6_cellid_df, on="id", how="right", validate="many_to_one" ) # %% -# Flatten id_msw grid → (y, x, id_msw) table, drop cells with no well +# Flatten id_sprinkling grid → (y, x, id_sprinkling) table, drop cells with no well grids = xr.merge([art_grid, svat]) -art_df = grids.to_dataframe().reset_index().query("(id_msw > 0) & (svat > 0)") +art_df = grids.to_dataframe().reset_index().query("(id_sprinkling > 0) & (svat > 0)") # Drop unnecessary columns. We preserve the x, y coords as they might prove # useful for debugging. art_df = art_df.drop(["dx", "dy"], axis=1) @@ -136,8 +136,8 @@ def double_length_df_subunit( # Join: each SVAT cell gets the matching well row(s) from arl_points msw_mf6_merged_df = art_df.merge( points_mf6_merged_df, # brings id back as a column - left_on="id_msw", - right_on="id2grid_p", + left_on="id_sprinkling", + right_on="id_sprinkling_p", how="inner", validate="many_to_one", ) diff --git a/imod/tests/test_msw/test_sprinkling.py b/imod/tests/test_msw/test_sprinkling.py index 1e1e229df..7e784d10c 100644 --- a/imod/tests/test_msw/test_sprinkling.py +++ b/imod/tests/test_msw/test_sprinkling.py @@ -67,7 +67,7 @@ class SprinklingPointsCaseData: x_p: Optional[np.ndarray] = None y_p: Optional[np.ndarray] = None layer_p: Optional[np.ndarray] = None - id2grid_p: Optional[np.ndarray] = None + id_sprinkling_p: Optional[np.ndarray] = None capacity_p: Optional[np.ndarray] = None @@ -207,7 +207,7 @@ def case_one_point_one_art_cell( case_data.x_p = [2.0, 2.0, 2.0] case_data.y_p = [3.0, 2.0, 1.0] case_data.layer_p = [1, 2, 3] - case_data.id2grid_p = [1, 2, 3] + case_data.id_sprinkling_p = [1, 2, 3] case_data.capacity_p = [10.0, 20.0, 30.0] expected_data = ExpectedCaseData() @@ -239,7 +239,7 @@ def case_multi_point_one_art_cell( case_data.x_p = [2.0, 2.0, 2.0] case_data.y_p = [3.0, 2.0, 1.0] case_data.layer_p = [1, 2, 3] - case_data.id2grid_p = [1, 1, 1] + case_data.id_sprinkling_p = [1, 1, 1] case_data.capacity_p = [10.0, 20.0, 30.0] expected_data = ExpectedCaseData() @@ -268,7 +268,7 @@ def case_one_point_multi_art_cell( case_data.x_p = [2.0] case_data.y_p = [2.0] case_data.layer_p = [2] - case_data.id2grid_p = [1] + case_data.id_sprinkling_p = [1] case_data.capacity_p = [10.0] expected_data = ExpectedCaseData() @@ -301,7 +301,7 @@ def case_art_grid_outside( case_data.x_p = [2.0] case_data.y_p = [2.0] case_data.layer_p = [3] - case_data.id2grid_p = [4] + case_data.id_sprinkling_p = [4] case_data.capacity_p = [40.0] expected_data = ExpectedCaseData() @@ -333,7 +333,7 @@ def case_point_outside( case_data.x_p = [3.0] case_data.y_p = [1.0] case_data.layer_p = [3] - case_data.id2grid_p = [5] + case_data.id_sprinkling_p = [5] case_data.capacity_p = [40.0] expected_data = ExpectedCaseData() @@ -482,7 +482,7 @@ def test_points_simple_model( case_data.x_p, case_data.y_p, case_data.layer_p, - case_data.id2grid_p, + case_data.id_sprinkling_p, case_data.capacity_p, ) @@ -557,7 +557,7 @@ def test_sprinkling_from_imod5_data__grid(cap_data_sprinkling_grid): @pytest.mark.unittest_jit def test_sprinklingpoints_from_imod5_data__points(cap_data_sprinkling_points): # Arrange - expected_vars = {"id2grid_p", "capacity_p", "layer_p", "y_p", "x_p", "id_msw"} + expected_vars = {"id_sprinkling_p", "capacity_p", "layer_p", "y_p", "x_p", "id_sprinkling"} # Act sprinkling = msw.SprinklingPoints.from_imod5_data(cap_data_sprinkling_points) From cfa5c63a76cc7a7c17fe2adbcb466bce34e4faba Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 19 Aug 2026 11:27:55 +0200 Subject: [PATCH 31/38] Format --- imod/msw/sprinkling.py | 4 +++- imod/tests/_scratch.py | 8 +++++++- imod/tests/test_msw/test_sprinkling.py | 9 ++++++++- 3 files changed, 18 insertions(+), 3 deletions(-) diff --git a/imod/msw/sprinkling.py b/imod/msw/sprinkling.py index 1f8bf1000..6fa7b2691 100644 --- a/imod/msw/sprinkling.py +++ b/imod/msw/sprinkling.py @@ -210,7 +210,9 @@ def _merge_sprinkling_points_with_grids( # Flatten id_sprinkling grid → (y, x, id_sprinkling) table, drop cells with no well grids = xr.merge([sprinkling_id_grid, svat]) - art_df = grids.to_dataframe().reset_index().query("(id_sprinkling > 0) & (svat > 0)") + art_df = ( + grids.to_dataframe().reset_index().query("(id_sprinkling > 0) & (svat > 0)") + ) # Drop unnecessary columns. We preserve the x, y coords as they might # prove useful for debugging. art_df = art_df.drop(["dx", "dy"], axis=1) diff --git a/imod/tests/_scratch.py b/imod/tests/_scratch.py index 912f6df7f..97aafa82a 100644 --- a/imod/tests/_scratch.py +++ b/imod/tests/_scratch.py @@ -99,7 +99,13 @@ def double_length_df_subunit( arl_points.columns = ["x_p", "y_p", "layer_p", "id_sprinkling_p", "capacity"] # Enforce dtypes arl_points = arl_points.astype( - {"x_p": float, "y_p": float, "layer_p": int, "id_sprinkling_p": int, "capacity": float} + { + "x_p": float, + "y_p": float, + "layer_p": int, + "id_sprinkling_p": int, + "capacity": float, + } ) arl_points["id"] = arl_points.index.astype(str) arl_points = arl_points.set_index("id") diff --git a/imod/tests/test_msw/test_sprinkling.py b/imod/tests/test_msw/test_sprinkling.py index 7e784d10c..0c43287e9 100644 --- a/imod/tests/test_msw/test_sprinkling.py +++ b/imod/tests/test_msw/test_sprinkling.py @@ -557,7 +557,14 @@ def test_sprinkling_from_imod5_data__grid(cap_data_sprinkling_grid): @pytest.mark.unittest_jit def test_sprinklingpoints_from_imod5_data__points(cap_data_sprinkling_points): # Arrange - expected_vars = {"id_sprinkling_p", "capacity_p", "layer_p", "y_p", "x_p", "id_sprinkling"} + expected_vars = { + "id_sprinkling_p", + "capacity_p", + "layer_p", + "y_p", + "x_p", + "id_sprinkling", + } # Act sprinkling = msw.SprinklingPoints.from_imod5_data(cap_data_sprinkling_points) From aa1dc245e03c26b1cdfaac2e48191c39358563ff Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 19 Aug 2026 15:51:52 +0200 Subject: [PATCH 32/38] Update unittest for LayeredWell.from_imod5_cap_data__points --- imod/tests/test_mf6/test_mf6_wel.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/imod/tests/test_mf6/test_mf6_wel.py b/imod/tests/test_mf6/test_mf6_wel.py index 94d240f87..6c975d745 100644 --- a/imod/tests/test_mf6/test_mf6_wel.py +++ b/imod/tests/test_mf6/test_mf6_wel.py @@ -1164,7 +1164,14 @@ def test_from_imod5_cap_data__big_grid( @pytest.mark.unittest_jit def test_from_imod5_cap_data__points(cap_data_sprinkling_points, cap_coupled_dis_grid): - with pytest.raises(NotImplementedError): - LayeredWell.from_imod5_cap_data( + # Act + well = LayeredWell.from_imod5_cap_data( cap_data_sprinkling_points, cap_coupled_dis_grid ) + # Assert + ds = well.dataset + np.testing.assert_allclose(ds["x"].to_numpy(), np.array([2.0, 2.0])) + np.testing.assert_allclose(ds["y"].to_numpy(), np.array([3.0, 2.0])) + np.testing.assert_equal(ds["layer"].to_numpy(), np.array([2, 3])) + np.testing.assert_allclose(ds["rate"].to_numpy(), np.array([0.0, 0.0])) + np.testing.assert_equal(ds["id"].to_numpy(), np.array(["0", "1"])) \ No newline at end of file From 31c781ac6dbe622b5fbea3697724e0fc151a2bac Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 19 Aug 2026 15:53:30 +0200 Subject: [PATCH 33/38] Implement subunits for art_grid --- imod/msw/sprinkling.py | 37 ++++--------- imod/tests/fixtures/imod5_cap_data.py | 6 +-- imod/tests/test_msw/test_sprinkling.py | 72 +++++++++++++++----------- 3 files changed, 54 insertions(+), 61 deletions(-) diff --git a/imod/msw/sprinkling.py b/imod/msw/sprinkling.py index 6fa7b2691..5ba5f384f 100644 --- a/imod/msw/sprinkling.py +++ b/imod/msw/sprinkling.py @@ -53,8 +53,13 @@ def _sprinkling_data_from_imod5_ipf( cap_data: CapSprinklingDataDict, ) -> SprinklingPointsGridDataDict: art_grid = cap_data["artificial_recharge"] - df_points = cap_data["artificial_recharge_layer"] + # Set urban landuse irrigation to 0, as sprinkling is not allowed for urban landuse. + subunit_template = xr.DataArray( + np.array([1, 0], dtype=int), dims="subunit", coords={"subunit": [0, 1]} + ) + art_grid = subunit_template * art_grid + df_points = cap_data["artificial_recharge_layer"] # Select first 5 columns and enforce column names, iMOD5 expects columns in # this order. The additional columns are metadata for the user and can be # ignored. @@ -141,29 +146,6 @@ def _extract_indexer_for_svat(df: pd.DataFrame, columns: list[str]): return indexer.T -def _replicate_dataframe_by_subunit( - df: pd.DataFrame, subunit_col: str = "subunit" -) -> pd.DataFrame: - """ - Duplicate the rows of a DataFrame for each subunit (0 and 1). - - Parameters - ---------- - df : pd.DataFrame - Input DataFrame to duplicate. - subunit_col : str, optional - Name of the column to assign subunit values, by default "subunit". - - Returns - ------- - pd.DataFrame - DataFrame with duplicated rows for each subunit. - """ - subunit_nrs = [0, 1] - df_ls = [df.assign(**{subunit_col: subunit_nr}) for subunit_nr in subunit_nrs] - return pd.concat(df_ls, ignore_index=True) - - def _get_mf6_cellid_dataframe(mf6_well: Mf6Wel) -> pd.DataFrame: """ Get cellids from the Mf6Wel objects dataset and convert to a dataframe for @@ -457,10 +439,7 @@ class SprinklingPoints(MetaSwapPackage, IRegridPackage): "trajectory": VariableMetaData(10, None, None, str), } - _with_subunit = ( - "max_abstraction_groundwater", - "max_abstraction_surfacewater", - ) + _with_subunit = ("id_sprinkling",) _without_subunit = () _to_fill = ( @@ -497,6 +476,8 @@ def __init__( art_grid = art_grid.rename("id_sprinkling") self.dataset = xr.merge([art_grid, points_ds]) + self._pkgcheck() + @classmethod def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "SprinklingPoints": """ diff --git a/imod/tests/fixtures/imod5_cap_data.py b/imod/tests/fixtures/imod5_cap_data.py index b670abdc1..8e9b29a9b 100644 --- a/imod/tests/fixtures/imod5_cap_data.py +++ b/imod/tests/fixtures/imod5_cap_data.py @@ -143,11 +143,11 @@ def cap_data_sprinkling_points() -> Imod5DataDict: artificial_rch_type = zeros_grid(n) artificial_rch_type[:, 1] = 3000 - artificial_rch_type[:, 2] = 4000 + artificial_rch_type[2, 1] = 4000 data = { - "x": [1.5, 2.5], - "y": [1.5, 2.5], + "x": [2.0, 2.0], + "y": [3.0, 2.0], "layer": [2, 3], "id": [3000, 4000], "capacity": [15.0, 30.0], diff --git a/imod/tests/test_msw/test_sprinkling.py b/imod/tests/test_msw/test_sprinkling.py index 0c43287e9..9c9852add 100644 --- a/imod/tests/test_msw/test_sprinkling.py +++ b/imod/tests/test_msw/test_sprinkling.py @@ -196,12 +196,15 @@ def case_one_point_one_art_cell( """ svat, _ = sprinkling_svat_index case_data = SprinklingPointsCaseData() - case_data.art_grid = xr.full_like(svat.isel(subunit=0, drop=True), 0, dtype=int) + case_data.art_grid = xr.full_like(svat, 0, dtype=int) # fmt: off case_data.art_grid.data = np.array( - [[0, 1, 0], - [0, 2, 0], - [0, 3, 0]] + [[[0, 1, 0], + [0, 2, 0], + [0, 3, 0],], + [[0, 1, 0], + [0, 2, 0], + [0, 3, 0]]] ) # fmt: on case_data.x_p = [2.0, 2.0, 2.0] @@ -228,12 +231,15 @@ def case_multi_point_one_art_cell( """ svat, _ = sprinkling_svat_index case_data = SprinklingPointsCaseData() - case_data.art_grid = xr.full_like(svat.isel(subunit=0, drop=True), 0, dtype=int) + case_data.art_grid = xr.full_like(svat, 0, dtype=int) # fmt: off case_data.art_grid.data = np.array( - [[0, 0, 0], - [0, 1, 0], - [0, 0, 0]] + [[[0, 0, 0], + [0, 1, 0], + [0, 0, 0]], + [[0, 0, 0], + [0, 1, 0], + [0, 0, 0]]] ) # fmt: on case_data.x_p = [2.0, 2.0, 2.0] @@ -257,12 +263,15 @@ def case_one_point_multi_art_cell( """ svat, _ = sprinkling_svat_index case_data = SprinklingPointsCaseData() - case_data.art_grid = xr.full_like(svat.isel(subunit=0, drop=True), 0, dtype=int) + case_data.art_grid = xr.full_like(svat, 0, dtype=int) # fmt: off case_data.art_grid.data = np.array( - [[0, 1, 0], - [0, 1, 0], - [0, 1, 0]] + [[[0, 1, 0], + [0, 1, 0], + [0, 1, 0]], + [[0, 1, 0], + [0, 1, 0], + [0, 1, 0]]] ) # fmt: on case_data.x_p = [2.0] @@ -290,13 +299,17 @@ def case_art_grid_outside( """ svat, _ = sprinkling_svat_index case_data = SprinklingPointsCaseData() - case_data.art_grid = xr.full_like(svat.isel(subunit=0, drop=True), 0, dtype=int) + case_data.art_grid = xr.full_like(svat, 0, dtype=int) # fmt: off case_data.art_grid.data = np.array( - [[0, 0, 4], - [0, 0, 0], - [0, 0, 0]] + [[[0, 0, 4], + [0, 0, 0], + [0, 0, 0]], + [[0, 0, 4], + [0, 0, 0], + [0, 0, 0]]] ) + # fmt: on case_data.x_p = [2.0] case_data.y_p = [2.0] @@ -322,12 +335,15 @@ def case_point_outside( """ svat, _ = sprinkling_svat_index case_data = SprinklingPointsCaseData() - case_data.art_grid = xr.full_like(svat.isel(subunit=0, drop=True), 0, dtype=int) + case_data.art_grid = xr.full_like(svat, 0, dtype=int) # fmt: off case_data.art_grid.data = np.array( - [[0, 0, 0], - [0, 0, 0], - [0, 5, 0]] + [[[0, 0, 0], + [0, 0, 0], + [0, 5, 0]], + [[0, 0, 0], + [0, 0, 0], + [0, 5, 0]]] ) # fmt: on case_data.x_p = [3.0] @@ -570,7 +586,7 @@ def test_sprinklingpoints_from_imod5_data__points(cap_data_sprinkling_points): sprinkling = msw.SprinklingPoints.from_imod5_data(cap_data_sprinkling_points) # Assert - assert sprinkling.dataset.sizes == {"id": 2, "x": 3, "y": 3} + assert sprinkling.dataset.sizes == {"subunit": 2, "id": 2, "x": 3, "y": 3} assert set(sprinkling.dataset.keys()) == expected_vars # No unit conversion is done in SprinklingPoints, as the capacity is already # in m3/d @@ -617,12 +633,8 @@ def test_sprinklingpoints_from_imod5_data_write__points( ) # Assert - np.testing.assert_equal(results["svat"], [1, 2, 3, 4]) - np.testing.assert_equal(results["svat_groundwater"], [1, 2, 4, 4]) - np.testing.assert_equal(results["layer"], [2, 2, 2, 2]) - np.testing.assert_equal( - results["max_abstraction_surfacewater"], [15.0, 15.0, 0.0, 0.0] - ) - np.testing.assert_equal( - results["max_abstraction_groundwater"], [0.0, 0.0, 15.0, 15.0] - ) + np.testing.assert_equal(results["svat"], [1, 2]) + np.testing.assert_equal(results["svat_groundwater"], [1, 2]) + np.testing.assert_equal(results["layer"], [2, 3]) + np.testing.assert_equal(results["max_abstraction_surfacewater"], [0.0, 30.0]) + np.testing.assert_equal(results["max_abstraction_groundwater"], [15.0, 0.0]) From 6d302d5e568725646cd91bc2b63efb54b20af6d3 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 19 Aug 2026 15:53:57 +0200 Subject: [PATCH 34/38] Format --- imod/tests/test_mf6/test_mf6_wel.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/imod/tests/test_mf6/test_mf6_wel.py b/imod/tests/test_mf6/test_mf6_wel.py index 6c975d745..04bbb953d 100644 --- a/imod/tests/test_mf6/test_mf6_wel.py +++ b/imod/tests/test_mf6/test_mf6_wel.py @@ -1166,12 +1166,12 @@ def test_from_imod5_cap_data__big_grid( def test_from_imod5_cap_data__points(cap_data_sprinkling_points, cap_coupled_dis_grid): # Act well = LayeredWell.from_imod5_cap_data( - cap_data_sprinkling_points, cap_coupled_dis_grid - ) + cap_data_sprinkling_points, cap_coupled_dis_grid + ) # Assert ds = well.dataset np.testing.assert_allclose(ds["x"].to_numpy(), np.array([2.0, 2.0])) np.testing.assert_allclose(ds["y"].to_numpy(), np.array([3.0, 2.0])) np.testing.assert_equal(ds["layer"].to_numpy(), np.array([2, 3])) np.testing.assert_allclose(ds["rate"].to_numpy(), np.array([0.0, 0.0])) - np.testing.assert_equal(ds["id"].to_numpy(), np.array(["0", "1"])) \ No newline at end of file + np.testing.assert_equal(ds["id"].to_numpy(), np.array(["0", "1"])) From 2e35dcbf892420f52d272b123efadb75cf1f7486 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 19 Aug 2026 15:58:16 +0200 Subject: [PATCH 35/38] Add test case where only one subunit is active --- imod/tests/test_msw/test_sprinkling.py | 36 ++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/imod/tests/test_msw/test_sprinkling.py b/imod/tests/test_msw/test_sprinkling.py index 9c9852add..06f8a2bcd 100644 --- a/imod/tests/test_msw/test_sprinkling.py +++ b/imod/tests/test_msw/test_sprinkling.py @@ -222,6 +222,42 @@ def case_one_point_one_art_cell( return case_data, expected_data + def case_one_point_one_art_cell__one_subunit( + self, sprinkling_svat_index + ) -> tuple[SprinklingPointsCaseData, ExpectedCaseData]: + """ + Simple test case for sprinkling points. Each point is mapped to one + svat. Only one subunit is used, similar to when imported from iMOD5 + DBASE + """ + svat, _ = sprinkling_svat_index + case_data = SprinklingPointsCaseData() + case_data.art_grid = xr.full_like(svat, 0, dtype=int) + # fmt: off + case_data.art_grid.data = np.array( + [[[0, 1, 0], + [0, 2, 0], + [0, 3, 0],], + [[0, 0, 0], + [0, 0, 0], + [0, 0, 0]]] + ) + # fmt: on + case_data.x_p = [2.0, 2.0, 2.0] + case_data.y_p = [3.0, 2.0, 1.0] + case_data.layer_p = [1, 2, 3] + case_data.id_sprinkling_p = [1, 2, 3] + case_data.capacity_p = [10.0, 20.0, 30.0] + + expected_data = ExpectedCaseData() + expected_data.svat = np.array([1, 2]) + expected_data.svat_gw = np.array([1, 2]) + expected_data.layer = np.array([1, 3]) + expected_data.abs_gw = np.array([10.0, 30.0]) + expected_data.abs_sw = np.array([0.0, 0.0]) + + return case_data, expected_data + def case_multi_point_one_art_cell( self, sprinkling_svat_index ) -> tuple[SprinklingPointsCaseData, ExpectedCaseData]: From baaf40ada2a0841f09e3fbe631d894763fa693a1 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 19 Aug 2026 16:07:40 +0200 Subject: [PATCH 36/38] Update changelog --- docs/api/changelog.rst | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst index cf97902cb..1a8fc6bc0 100644 --- a/docs/api/changelog.rst +++ b/docs/api/changelog.rst @@ -9,6 +9,16 @@ The format is based on `Keep a Changelog`_, and this project adheres to [Unreleased] ------------ +Added +~~~~~ + +- :class:`imod.msw.SprinklingPoints` to specify sprinkling from points for + MetaSWAP models, instead of from grid. You can use this to specify sprinkling + wells from IPF files in an iMOD5 CAP dataset with + :meth:`imod.msw.SprinklingPoints.from_imod5_cap_data`. +- :class:`imod.msw.LayeredWell.from_imod5_cap_data` now also supports loading + wells from IPF files in an iMOD5 CAP dataset. + Fixed ~~~~~ From c8981d20c0037a8f5a06038c6dfc7ce781d4444c Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 19 Aug 2026 17:00:59 +0200 Subject: [PATCH 37/38] Move sprinkling imod5 import utilities into the imod5_converter utilities --- imod/msw/model.py | 8 +- imod/msw/sprinkling.py | 115 +++----------------------- imod/msw/utilities/imod5_converter.py | 111 ++++++++++++++++++++++++- 3 files changed, 127 insertions(+), 107 deletions(-) diff --git a/imod/msw/model.py b/imod/msw/model.py index f295b76d2..c4c349d3c 100644 --- a/imod/msw/model.py +++ b/imod/msw/model.py @@ -46,11 +46,12 @@ from imod.msw.ponding import Ponding from imod.msw.regrid.regrid_schemes import CapDataRegridMethod from imod.msw.scaling_factors import ScalingFactors -from imod.msw.sprinkling import Sprinkling +from imod.msw.sprinkling import Sprinkling, SprinklingPoints from imod.msw.timeutil import to_metaswap_timeformat from imod.msw.utilities.common import find_in_file_list from imod.msw.utilities.imod5_converter import ( has_active_scaling_factor, + is_sprinkling_from_points, ) from imod.msw.utilities.mask import ( MetaSwapActive, @@ -830,7 +831,10 @@ def from_imod5_data( } model["infiltration"] = Infiltration.from_imod5_data(imod5_masked) model["ponding"] = Ponding.from_imod5_data(imod5_masked) - model["sprinkling"] = Sprinkling.from_imod5_data(imod5_masked) + if is_sprinkling_from_points(imod5_masked): + model["sprinkling"] = SprinklingPoints.from_imod5_data(imod5_masked) + else: + model["sprinkling"] = Sprinkling.from_imod5_data(imod5_masked) model["meteo_grid"] = MeteoGridCopy.from_imod5_data(imod5_masked) model["prec_mapping"] = PrecipitationMapping.from_imod5_data(imod5_masked) model["evt_mapping"] = EvapotranspirationMapping.from_imod5_data(imod5_masked) diff --git a/imod/msw/sprinkling.py b/imod/msw/sprinkling.py index 5ba5f384f..a68f7f291 100644 --- a/imod/msw/sprinkling.py +++ b/imod/msw/sprinkling.py @@ -1,5 +1,5 @@ import textwrap -from typing import TextIO, TypedDict, cast +from typing import TextIO, cast import numpy as np import pandas as pd @@ -14,32 +14,13 @@ SprinklingPointsRegridMethod, SprinklingRegridMethod, ) -from imod.msw.utilities.common import concat_imod5 from imod.msw.utilities.imod5_converter import ( - get_cell_area_from_imod5_data, + CapSprinklingDataDict, + is_sprinkling_from_points, + sprinkling_data_from_imod5_grid, + sprinkling_data_from_imod5_ipf, ) -from imod.typing import GridDataArray, GridDataDict, Imod5DataDict, IntArray -from imod.typing.grid import zeros_like - - -# Some additional type aliases for sprinkling data, which is a bit more complex -# than other packages. -class CapSprinklingDataDict(TypedDict, total=False): - artificial_recharge: GridDataArray - artificial_recharge_layer: pd.DataFrame - artificial_recharge_capacity: GridDataArray - - -class SprinklingPointsDataDict(TypedDict, total=False): - x_p: np.ndarray | list[float] - y_p: np.ndarray | list[float] - layer_p: np.ndarray | list[int] - id_sprinkling_p: np.ndarray | list[int] - capacity_p: np.ndarray | list[float] - - -class SprinklingPointsGridDataDict(SprinklingPointsDataDict, total=False): - art_grid: GridDataArray +from imod.typing import Imod5DataDict, IntArray def _ravel_per_subunit(da: xr.DataArray) -> np.ndarray: @@ -49,77 +30,6 @@ def _ravel_per_subunit(da: xr.DataArray) -> np.ndarray: return array_out[np.isfinite(array_out)] -def _sprinkling_data_from_imod5_ipf( - cap_data: CapSprinklingDataDict, -) -> SprinklingPointsGridDataDict: - art_grid = cap_data["artificial_recharge"] - # Set urban landuse irrigation to 0, as sprinkling is not allowed for urban landuse. - subunit_template = xr.DataArray( - np.array([1, 0], dtype=int), dims="subunit", coords={"subunit": [0, 1]} - ) - art_grid = subunit_template * art_grid - - df_points = cap_data["artificial_recharge_layer"] - # Select first 5 columns and enforce column names, iMOD5 expects columns in - # this order. The additional columns are metadata for the user and can be - # ignored. - arl_points = df_points.iloc[:, :5] - arl_points.columns = ["x_p", "y_p", "layer_p", "id_sprinkling_p", "capacity_p"] - # Enforce dtypes - dtype_dict = { - "x_p": float, - "y_p": float, - "layer_p": int, - "id_sprinkling_p": int, - "capacity_p": float, - } - - arl_points = arl_points.astype(dtype_dict) - arl_point_dict = cast( - SprinklingPointsDataDict, - {key: arl_points[key].to_numpy() for key in dtype_dict.keys()}, - ) - - return { - "art_grid": art_grid, - **arl_point_dict, - } - - -def _sprinkling_data_from_imod5_grid(cap_data: GridDataDict) -> GridDataDict: - # Convert units from mm/d to m3/d - msw_area = get_cell_area_from_imod5_data(cap_data) - capacity_mmd = cap_data["artificial_recharge_capacity"] - capacity_m3d = capacity_mmd * 1e-3 * msw_area.sel(subunit=0, drop=True) - - artificial_rch_type = cap_data["artificial_recharge"] - from_groundwater = artificial_rch_type == 1 - from_surfacewater = artificial_rch_type == 2 - is_active = artificial_rch_type != 0 - - zero_where_active = zeros_like(artificial_rch_type).where(is_active) - - # Add zero where active, to have active cells set to 0.0. - max_abstraction_groundwater_rural = zero_where_active.where( - ~from_groundwater, capacity_m3d - ) - max_abstraction_surfacewater_rural = zero_where_active.where( - ~from_surfacewater, capacity_m3d - ) - - # No sprinkling for urban environments - max_abstraction_urban = zero_where_active - - data = {} - data["max_abstraction_groundwater"] = concat_imod5( - max_abstraction_groundwater_rural, max_abstraction_urban - ) - data["max_abstraction_surfacewater"] = concat_imod5( - max_abstraction_surfacewater_rural, max_abstraction_urban - ) - return data - - def _extract_indexer_for_svat(df: pd.DataFrame, columns: list[str]): """ Get the indexer for a dataframe of wells to select the SVAT subunit for each @@ -384,8 +294,7 @@ def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "Sprinkling": ------- Sprinkling package """ - cap_data = imod5_data["cap"] - if isinstance(cap_data["artificial_recharge_layer"], pd.DataFrame): + if is_sprinkling_from_points(imod5_data): msg = textwrap.dedent( """ Unsupported format for artificial_recharge_layer: expected a @@ -394,8 +303,8 @@ def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "Sprinkling": """ ) raise TypeError(msg) - - data = _sprinkling_data_from_imod5_grid(cap_data) + cap_data = imod5_data["cap"] + data = sprinkling_data_from_imod5_grid(cap_data) return cls(**data) @@ -512,9 +421,9 @@ def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "SprinklingPoints": ------- SprinklingPoints package """ - cap_data = cast(CapSprinklingDataDict, imod5_data["cap"]) - if isinstance(cap_data["artificial_recharge_layer"], pd.DataFrame): - data = _sprinkling_data_from_imod5_ipf(cap_data) + if is_sprinkling_from_points(imod5_data): + cap_data = cast(CapSprinklingDataDict, imod5_data["cap"]) + data = sprinkling_data_from_imod5_ipf(cap_data) return cls(**data) else: msg = textwrap.dedent( diff --git a/imod/msw/utilities/imod5_converter.py b/imod/msw/utilities/imod5_converter.py index c6a3cd32d..6656488e8 100644 --- a/imod/msw/utilities/imod5_converter.py +++ b/imod/msw/utilities/imod5_converter.py @@ -1,3 +1,8 @@ +from typing import TypedDict, cast + +import numpy as np +import pandas as pd +import xarray as xr from xarray.core.utils import is_scalar from imod.common.constants import MaskValues @@ -5,11 +10,31 @@ from imod.mf6 import StructuredDiscretization from imod.msw.utilities.common import concat_imod5 from imod.msw.utilities.mask import MetaSwapActive -from imod.typing import GridDataArray, GridDataDict -from imod.typing.grid import ones_like +from imod.typing import GridDataArray, GridDataDict, Imod5DataDict +from imod.typing.grid import ones_like, zeros_like from imod.util.spatial import get_cell_area +# Some additional type aliases for sprinkling data, which is a bit more complex +# than other packages. +class CapSprinklingDataDict(TypedDict, total=False): + artificial_recharge: GridDataArray + artificial_recharge_layer: pd.DataFrame + artificial_recharge_capacity: GridDataArray + + +class SprinklingPointsDataDict(TypedDict, total=False): + x_p: np.ndarray | list[float] + y_p: np.ndarray | list[float] + layer_p: np.ndarray | list[int] + id_sprinkling_p: np.ndarray | list[int] + capacity_p: np.ndarray | list[float] + + +class SprinklingPointsGridDataDict(SprinklingPointsDataDict, total=False): + art_grid: GridDataArray + + def get_cell_area_from_imod5_data( imod5_cap: GridDataDict, ) -> GridDataArray: @@ -116,3 +141,85 @@ def has_active_scaling_factor(imod5_cap: GridDataDict): ) return not scaling_factor_inactive + + +def is_sprinkling_from_points(imod5_data: Imod5DataDict) -> bool: + """ + Check if sprinkling is specified from points, based on the presence of + sprinkling layer and sprinkling points data in the iMOD5 CAP dataset. + """ + cap_data = cast(CapSprinklingDataDict, imod5_data["cap"]) + if isinstance(cap_data.get("artificial_recharge_layer"), pd.DataFrame): + return True + return False + + +def sprinkling_data_from_imod5_ipf( + cap_data: CapSprinklingDataDict, +) -> SprinklingPointsGridDataDict: + art_grid = cap_data["artificial_recharge"] + # Set urban landuse irrigation to 0, as sprinkling is not allowed for urban landuse. + subunit_template = xr.DataArray( + np.array([1, 0], dtype=int), dims="subunit", coords={"subunit": [0, 1]} + ) + art_grid = subunit_template * art_grid + + df_points = cap_data["artificial_recharge_layer"] + # Select first 5 columns and enforce column names, iMOD5 expects columns in + # this order. The additional columns are metadata for the user and can be + # ignored. + arl_points = df_points.iloc[:, :5] + arl_points.columns = ["x_p", "y_p", "layer_p", "id_sprinkling_p", "capacity_p"] + # Enforce dtypes + dtype_dict = { + "x_p": float, + "y_p": float, + "layer_p": int, + "id_sprinkling_p": int, + "capacity_p": float, + } + + arl_points = arl_points.astype(dtype_dict) + arl_point_dict = cast( + SprinklingPointsDataDict, + {key: arl_points[key].to_numpy() for key in dtype_dict.keys()}, + ) + + return { + "art_grid": art_grid, + **arl_point_dict, + } + + +def sprinkling_data_from_imod5_grid(cap_data: GridDataDict) -> GridDataDict: + # Convert units from mm/d to m3/d + msw_area = get_cell_area_from_imod5_data(cap_data) + capacity_mmd = cap_data["artificial_recharge_capacity"] + capacity_m3d = capacity_mmd * 1e-3 * msw_area.sel(subunit=0, drop=True) + + artificial_rch_type = cap_data["artificial_recharge"] + from_groundwater = artificial_rch_type == 1 + from_surfacewater = artificial_rch_type == 2 + is_active = artificial_rch_type != 0 + + zero_where_active = zeros_like(artificial_rch_type).where(is_active) + + # Add zero where active, to have active cells set to 0.0. + max_abstraction_groundwater_rural = zero_where_active.where( + ~from_groundwater, capacity_m3d + ) + max_abstraction_surfacewater_rural = zero_where_active.where( + ~from_surfacewater, capacity_m3d + ) + + # No sprinkling for urban environments + max_abstraction_urban = zero_where_active + + data = {} + data["max_abstraction_groundwater"] = concat_imod5( + max_abstraction_groundwater_rural, max_abstraction_urban + ) + data["max_abstraction_surfacewater"] = concat_imod5( + max_abstraction_surfacewater_rural, max_abstraction_urban + ) + return data From 2c741e0ff195554be915ddc15aee6a85183615b0 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 19 Aug 2026 17:15:41 +0200 Subject: [PATCH 38/38] Reduce script to latest method calls --- imod/tests/_scratch.py | 185 ++--------------------------------------- 1 file changed, 8 insertions(+), 177 deletions(-) diff --git a/imod/tests/_scratch.py b/imod/tests/_scratch.py index 97aafa82a..bcc6bcda7 100644 --- a/imod/tests/_scratch.py +++ b/imod/tests/_scratch.py @@ -1,82 +1,11 @@ # %% -import pandas as pd -import xarray as xr - import imod from imod.util.dims import drop_layer_dim_cap_data - -def get_indexer(df: pd.DataFrame, columns: list[str]): - """ - Get the indexer for a dataframe of wells to select the SVAT subunit for each - well based on its row/col location in the model grid. - - Parameters - ---------- - df : pd.DataFrame - DataFrame containing the wells with columns "subunit", "row", - and "column" (1-based). - - Returns - ------- - np.ndarray - Indexer array for selecting SVAT subunits from svat_da. - """ - if not set("row", "column").issubset(columns): - raise ValueError("columns must contain 'row' and 'column'") - df.loc[:, ["row", "column"]] -= 1 # Convert to 0-based indexing for xarray - - indexer = df.loc[:, columns].to_numpy() - return indexer.T - - -def double_length_df_subunit( - df: pd.DataFrame, subunit_col: str = "subunit" -) -> pd.DataFrame: - """ - Duplicate the rows of a DataFrame for each subunit (0 and 1). - - Parameters - ---------- - df : pd.DataFrame - Input DataFrame to duplicate. - subunit_col : str, optional - Name of the column to assign subunit values, by default "subunit". - - Returns - ------- - pd.DataFrame - DataFrame with duplicated rows for each subunit. - """ - return pd.concat( - [df.assign(**{subunit_col: 0}), df.assign(**{subunit_col: 1})], - ignore_index=True, - ) - - -# %% -df_points = imod.ipf.read( - r"c:\Users\engelen\projects_wdir\imod-python\imod5_converter\NHI_sprint\Peelvenen\BASIS7\METASWAP\grid\Sprinkling\BEREGEN_LOC.IPF" -) -art_grid = ( - imod.idf.open( - r"c:\Users\engelen\projects_wdir\imod-python\imod5_converter\NHI_sprint\Peelvenen\BASIS7\METASWAP\grid\Sprinkling\BEREGENINGS_LOCATIES.IDF" - ) - .compute() - .astype(int) -) - -imod5_data = { - "cap": {"artificial_recharge": art_grid, "artificial_recharge_layer": df_points} -} - -well = imod.mf6.LayeredWell.from_imod5_cap_data( - imod5_data, target_dis=None, regridder_types=None, regrid_cache=None -) # %% # Setup to get example args for _render() of SprinklingPoints prj_data, period_data = imod.formats.prj.open_projectfile_data( - r"c:\Users\engelen\projects_wdir\imod-python\imod5_converter\NHI_sprint\Peelvenen\prjfiles\Peelvenen_relative_paths_dis_npf.PRJ" + r"c:\Users\engelen\projects_wdir\imod-python\imod5_converter\NHI_sprint\Peelvenen\prjfiles\Peelvenen_absolute_paths_sprinkling_ipf.PRJ" ) dis_pkg = imod.mf6.StructuredDiscretization.from_imod5_data(prj_data, validate=False) @@ -87,120 +16,22 @@ def double_length_df_subunit( prj_data = drop_layer_dim_cap_data(prj_data) griddata, msw_active = imod.msw.GridData.from_imod5_data(prj_data, dis_pkg) + +well = imod.mf6.LayeredWell.from_imod5_cap_data( + prj_data, target_dis=None, regridder_types=None, regrid_cache=None +) # Convert to args of Sprinkling._render() mf6_well = well.to_mf6_pkg( dis_pkg["idomain"], dis_pkg["top"], dis_pkg["bottom"], npf_pkg["k"] ) isactive_1d, svat = griddata.generate_isactive_svat_arrays() -# %% -# In from_imod5_cap_data -arl_points = df_points.iloc[:, :5] -arl_points.columns = ["x_p", "y_p", "layer_p", "id_sprinkling_p", "capacity"] -# Enforce dtypes -arl_points = arl_points.astype( - { - "x_p": float, - "y_p": float, - "layer_p": int, - "id_sprinkling_p": int, - "capacity": float, - } -) -arl_points["id"] = arl_points.index.astype(str) -arl_points = arl_points.set_index("id") - -points_ds = arl_points.to_xarray() - -# in def __init__ -art_grid = art_grid.rename("id_sprinkling") -dataset = xr.merge([art_grid, points_ds]) +sprinkling_points = imod.msw.SprinklingPoints.from_imod5_data(prj_data) # %% -# In render() -# Promote id to dim to join datasets -mf6_well_ds = mf6_well.dataset.set_coords("id").swap_dims({"ncellid": "id"}) -# Convert the cellid DataArray to a broad table for easier manipulation. -mf6_cellid_df = mf6_well_ds["cellid"].to_dataset("dim_cellid").to_dataframe() -# Select only the cellid columns we need and reset index to promote id to column -# for merging -dim_cellid = ["layer", "row", "column"] -mf6_cellid_df = mf6_cellid_df.loc[:, dim_cellid].reset_index() -# Merge again to confine to wells actually used in the modflow6 model. -# This drops points that are not in art_grid. -points_mf6_merged_df = arl_points.reset_index().merge( - mf6_cellid_df, on="id", how="right", validate="many_to_one" -) -# %% -# Flatten id_sprinkling grid → (y, x, id_sprinkling) table, drop cells with no well -grids = xr.merge([art_grid, svat]) -art_df = grids.to_dataframe().reset_index().query("(id_sprinkling > 0) & (svat > 0)") -# Drop unnecessary columns. We preserve the x, y coords as they might prove -# useful for debugging. -art_df = art_df.drop(["dx", "dy"], axis=1) +directory = r"c:\Users\engelen\projects_wdir\imod-python\imod5_converter\NHI_sprint\Peelvenen\conversion_output_ipf" +sprinkling_points.write(directory, isactive_1d, svat, dis_pkg, mf6_well) -# Join: each SVAT cell gets the matching well row(s) from arl_points -msw_mf6_merged_df = art_df.merge( - points_mf6_merged_df, # brings id back as a column - left_on="id_sprinkling", - right_on="id_sprinkling_p", - how="inner", - validate="many_to_one", -) - -# %% -# Derive the SVAT subunit for each well from the SVAT grid based on the -# well's row/col location. -indexer = get_indexer(msw_mf6_merged_df, columns=["subunit", "row", "column"]) -# We need to align the SVAT grid with the dis_pkg grid as the SVAT grid might be -# smaller. -idomain_flat = dis_pkg.dataset["idomain"].isel(layer=0, drop=True) -_, svat_aligned = xr.align(idomain_flat, svat, join="left") -svat_groundwater = svat_aligned.data[*indexer] - -msw_mf6_merged_df["svat_groundwater"] = svat_groundwater.astype(int) - -# %% -# Select columns that need to be written to scap_svat.inp -dataframe = msw_mf6_merged_df[["svat", "layer", "svat_groundwater"]] -dataframe["svat"] = dataframe["svat"].astype(int) -capacity = msw_mf6_merged_df["capacity"] -# Set wells with layer > 0 to groundwater abstraction, and wells with layer = 0 -# to surfacewater abstraction. -is_gw_extraction = msw_mf6_merged_df["layer"] > 0 -dataframe["max_abstraction_groundwater"] = capacity.where(is_gw_extraction, 0.0) -dataframe["max_abstraction_surfacewater"] = capacity.where(~is_gw_extraction, 0.0) - -# %% -# -# Deal with edge case: wells that are outside art_grid, but in model domain. -# These will be assigned to surfacewater abstraction. We can identify these by -# checking which wells in mf6_cellid_df are not in msw_mf6_merged_df. -is_outside_art = ~mf6_cellid_df["id"].isin(msw_mf6_merged_df["id"].unique()) -outside_df = points_mf6_merged_df.loc[is_outside_art, ["layer", "capacity"]] -outside_df = double_length_df_subunit(outside_df) -# Select the SVAT subunit for these wells based on their row/col location. -df_cellid_outside = mf6_cellid_df.loc[is_outside_art] -df_cellid_outside = double_length_df_subunit(df_cellid_outside) -indexer_outside = get_indexer(df_cellid_outside, columns=["subunit", "row", "column"]) -svat_outside = svat_aligned.data[ - indexer_outside[0], indexer_outside[1], indexer_outside[2] -] -outside_df["svat_groundwater"] = svat_outside.astype(int) -outside_df["svat"] = svat_outside.astype(int) -# Set capacity to surfacewater abstraction, and set groundwater abstraction to 0. -outside_df = outside_df.rename(columns={"capacity": "max_abstraction_surfacewater"}) -outside_df["max_abstraction_groundwater"] = 0.0 -# drop subunit column as it is no longer needed -outside_df = outside_df.drop(columns=["subunit"]) -# drop wells that are outside the active metaswap model domain (svat = 0) -outside_df = outside_df.query("svat > 0").reset_index(drop=True) - -# %% -# -# Combine -dataframe_out = pd.concat([dataframe, outside_df], axis=0, ignore_index=True) -dataframe_out = dataframe_out.sort_values(by=["svat"]).reset_index(drop=True) # %% # # TODO: Verify if iMOD5 SVAT grid the same as the one derived in this script.