From 4930661553bf086c854f0a3ed042a2d0edbf046c Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Tue, 21 Apr 2026 16:59:16 +0200 Subject: [PATCH 1/6] Allow user-supplied target grid. When not supplied, default to iMOD5's default setting instead of computing a target grid based on the smallest possible DIS grid. --- docs/faq/imod5_backwards_compatibility.rst | 6 ---- imod/common/utilities/grid.py | 19 ------------ imod/mf6/dis.py | 8 +++-- imod/mf6/model_gwf.py | 7 ++++- imod/mf6/simulation.py | 5 ++++ .../test_mf6/test_utilities/test_grid.py | 29 ------------------- 6 files changed, 17 insertions(+), 57 deletions(-) delete mode 100644 imod/tests/test_mf6/test_utilities/test_grid.py diff --git a/docs/faq/imod5_backwards_compatibility.rst b/docs/faq/imod5_backwards_compatibility.rst index 9698a6ee8..1467dc144 100644 --- a/docs/faq/imod5_backwards_compatibility.rst +++ b/docs/faq/imod5_backwards_compatibility.rst @@ -51,12 +51,6 @@ Notes preset from MODFLOW 6 ("Moderate") is set. This is because the solvers between iMODLFOW and MODFLOW 6 are different. You are advised to test settings yourself. -- The imported iMOD5 discretization for the model is created by taking the - smallest grid and finest resolution amongst the TOP, BOT, and BND grids. This - differs from iMOD5, where the first BND grid is used as target grid. All input - grids are regridded towards this target grid. Therefore, be careful when you - have a very fine resolution in one of these packages. - Files ----- diff --git a/imod/common/utilities/grid.py b/imod/common/utilities/grid.py index a6ca79e22..bcf0e50ec 100644 --- a/imod/common/utilities/grid.py +++ b/imod/common/utilities/grid.py @@ -81,22 +81,3 @@ def create_geometric_grid_info(active: xr.DataArray) -> pd.DataFrame: "dy": dy.flatten(), } ) - - -def create_smallest_target_grid(*grids: xr.DataArray) -> xr.DataArray: - """ - Create smallest target grid from multiple structured grids. This is the grid - with smallest extent and finest resolution amongst all provided grids. - """ - dx_ls, xmin_ls, xmax_ls, dy_ls, ymin_ls, ymax_ls = zip( - *[imod.util.spatial.spatial_reference(grid) for grid in grids] - ) - - dx = min(dx_ls) - xmin = max(xmin_ls) - xmax = min(xmax_ls) - dy = max(dy_ls) - ymax = min(ymax_ls) - ymin = max(ymin_ls) - - return imod.util.spatial.empty_2d(dx, xmin, xmax, dy, ymin, ymax) diff --git a/imod/mf6/dis.py b/imod/mf6/dis.py index e2878044b..ac97b56ff 100644 --- a/imod/mf6/dis.py +++ b/imod/mf6/dis.py @@ -7,7 +7,6 @@ from imod.common.interfaces.imaskingsettings import IMaskingSettings from imod.common.interfaces.iregridpackage import IRegridPackage from imod.common.utilities.dataclass_type import DataclassType -from imod.common.utilities.grid import create_smallest_target_grid from imod.common.utilities.regrid import ( _regrid_like, _regrid_package_data, @@ -178,6 +177,7 @@ def from_imod5_data( regridder_types: Optional[DataclassType] = None, regrid_cache: RegridderWeightsCache = RegridderWeightsCache(), validate: bool = True, + target_grid: Optional[GridDataArray] = None, ) -> "StructuredDiscretization": """ Construct package from iMOD5 data, loaded with the @@ -202,6 +202,9 @@ def from_imod5_data( regrid_cache: RegridderWeightsCache, optional stores regridder weights for different regridders. Can be used to speed up regridding, if the same regridders are used several times for regridding different arrays. + target_grid: GridDataArray, optional + the target grid to which the data should be regridded. If not + provided, the first grid of the BND package is used as target grid. Returns ------- @@ -214,7 +217,8 @@ def from_imod5_data( "bottom": imod5_data["bot"]["bottom"], } - target_grid = create_smallest_target_grid(*data.values()) + if target_grid is None: + target_grid = data["idomain"].isel(layer=0, drop=True, missing_dims="ignore") if regridder_types is None: regridder_types = StructuredDiscretization.get_regrid_methods() diff --git a/imod/mf6/model_gwf.py b/imod/mf6/model_gwf.py index 201996272..e14997d21 100644 --- a/imod/mf6/model_gwf.py +++ b/imod/mf6/model_gwf.py @@ -207,6 +207,7 @@ def from_imod5_data( allocation_options: Optional[SimulationAllocationOptions] = None, distributing_options: Optional[SimulationDistributingOptions] = None, regridder_types: Optional[dict[str, DataclassType]] = None, + target_grid: Optional[GridDataArray] = None, ) -> "GroundwaterFlowModel": """ Imports a GroundwaterFlowModel (GWF) from the data in an iMOD5 project @@ -244,7 +245,10 @@ def from_imod5_data( then it should be imported separately regridder_types: dict[str, RegridMethodType] the key is the package name. The value is a subclass of RegridMethodType. - + target_grid: GridDataArray, optional + the target grid to which the data should be regridded. If not + provided, the first grid of the BND package is used as target grid. + Returns ------- A GWF model containing the packages that could be imported form IMOD5. Users must still @@ -268,6 +272,7 @@ def from_imod5_data( cast(DiscretizationRegridMethod, regridder_types.get("dis")), regrid_cache, False, + target_grid=target_grid, ) grid = dis_pkg.dataset["idomain"] result["dis"] = dis_pkg diff --git a/imod/mf6/simulation.py b/imod/mf6/simulation.py index 1c9b32801..e2ff57b05 100644 --- a/imod/mf6/simulation.py +++ b/imod/mf6/simulation.py @@ -1751,6 +1751,7 @@ def from_imod5_data( allocation_options: Optional[SimulationAllocationOptions] = None, distributing_options: Optional[SimulationDistributingOptions] = None, regridder_types: Optional[dict[str, DataclassType]] = None, + target_grid: Optional[GridDataArray] = None, ) -> "Modflow6Simulation": """ Imports a GroundwaterFlowModel (GWF) from the data in an iMOD5 project @@ -1795,6 +1796,9 @@ def from_imod5_data( the key is the package name. The value is the RegridMethodType object containing the settings for regridding the package with the specified key. + target_grid: GridDataArray, optional + the target grid to which the data should be regridded. If not + provided, the first grid of the BND package is used as target grid. Returns ------- @@ -1848,6 +1852,7 @@ def from_imod5_data( allocation_options, distributing_options, regridder_types, + target_grid, ) simulation["imported_model"] = gwf_model diff --git a/imod/tests/test_mf6/test_utilities/test_grid.py b/imod/tests/test_mf6/test_utilities/test_grid.py deleted file mode 100644 index cc6d7b0a2..000000000 --- a/imod/tests/test_mf6/test_utilities/test_grid.py +++ /dev/null @@ -1,29 +0,0 @@ -import numpy as np - -from imod.common.utilities.grid import create_smallest_target_grid -from imod.util.spatial import empty_2d - - -def test_create_smallest_target_grid(): - # Three grids with aligned cell edges at each 100m - grid1 = empty_2d(dx=25.0, xmin=100.0, xmax=300.0, dy=-25.0, ymin=100.0, ymax=300.0) - grid2 = empty_2d(dx=50.0, xmin=0.0, xmax=200.0, dy=-50.0, ymin=0.0, ymax=200.0) - grid3 = empty_2d(dx=20.0, xmin=0.0, xmax=400.0, dy=-20.0, ymin=0.0, ymax=400.0) - - actual = create_smallest_target_grid(grid1, grid2, grid3) - - assert actual.coords["dx"] == 20.0 - assert actual.coords["dy"] == -20.0 - assert np.all(actual.coords["x"].values == [110.0, 130.0, 150.0, 170.0, 190.0]) - assert np.all(actual.coords["y"].values == [190.0, 170.0, 150.0, 130.0, 110.0]) - - # Two grids with barely aligned cell edges. - grid1 = empty_2d(dx=50.0, xmin=110.0, xmax=220.0, dy=-50.0, ymin=110.0, ymax=220.0) - grid2 = empty_2d(dx=20.0, xmin=0.0, xmax=400.0, dy=-20.0, ymin=0.0, ymax=400.0) - - actual = create_smallest_target_grid(grid1, grid2) - - assert actual.coords["dx"] == 20.0 - assert actual.coords["dy"] == -20.0 - assert np.all(actual.coords["x"].values == [120.0, 140.0, 160.0, 180.0, 200.0]) - assert np.all(actual.coords["y"].values == [210.0, 190.0, 170.0, 150.0, 130.0]) From cbda0c0c8b910a7d24f71e95d3b7778c5b52a447 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Tue, 21 Apr 2026 17:03:36 +0200 Subject: [PATCH 2/6] Format --- imod/mf6/dis.py | 4 +++- imod/mf6/model_gwf.py | 2 +- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/imod/mf6/dis.py b/imod/mf6/dis.py index ac97b56ff..01fec46a2 100644 --- a/imod/mf6/dis.py +++ b/imod/mf6/dis.py @@ -218,7 +218,9 @@ def from_imod5_data( } if target_grid is None: - target_grid = data["idomain"].isel(layer=0, drop=True, missing_dims="ignore") + target_grid = data["idomain"].isel( + layer=0, drop=True, missing_dims="ignore" + ) if regridder_types is None: regridder_types = StructuredDiscretization.get_regrid_methods() diff --git a/imod/mf6/model_gwf.py b/imod/mf6/model_gwf.py index e14997d21..5e7b544af 100644 --- a/imod/mf6/model_gwf.py +++ b/imod/mf6/model_gwf.py @@ -248,7 +248,7 @@ def from_imod5_data( target_grid: GridDataArray, optional the target grid to which the data should be regridded. If not provided, the first grid of the BND package is used as target grid. - + Returns ------- A GWF model containing the packages that could be imported form IMOD5. Users must still From e4fa7b8e3d792ba356409e43dbbe30404ead760b Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 22 Apr 2026 11:16:49 +0200 Subject: [PATCH 3/6] Update docstring --- imod/mf6/dis.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/imod/mf6/dis.py b/imod/mf6/dis.py index 01fec46a2..03b0de26c 100644 --- a/imod/mf6/dis.py +++ b/imod/mf6/dis.py @@ -183,9 +183,8 @@ def from_imod5_data( Construct package from iMOD5 data, loaded with the :func:`imod.formats.prj.open_projectfile_data` function. - Method regrids all variables to a target grid with the smallest extent - and smallest cellsize available in all the grids. Consequently it - converts iMODFLOW data to MODFLOW 6 data. + Method regrids all variables to a target grid. If not provided, the + first grid of the BND package is used as target grid. .. note:: From 59a2254cc4f003a4c8c03a8304f0301442d28f9e Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 22 Apr 2026 11:20:57 +0200 Subject: [PATCH 4/6] Add unittests to verify target_grid arg does what it should do --- imod/tests/test_mf6/test_mf6_dis.py | 36 ++++++++++++++++++++-- imod/tests/test_mf6/test_mf6_simulation.py | 28 +++++++++++++++-- 2 files changed, 59 insertions(+), 5 deletions(-) diff --git a/imod/tests/test_mf6/test_mf6_dis.py b/imod/tests/test_mf6/test_mf6_dis.py index 1ad1fa2ed..2822baf99 100644 --- a/imod/tests/test_mf6/test_mf6_dis.py +++ b/imod/tests/test_mf6/test_mf6_dis.py @@ -336,12 +336,13 @@ def test_from_imod5_data__idomain_values(imod5_dataset): assert (dis["idomain"] == 2).sum() == 688329 -def test_from_imod5_data__grid_extent(imod5_dataset): +def test_from_imod5_data_regridding__default(imod5_dataset): imod5_data = imod5_dataset[0] dis = imod.mf6.StructuredDiscretization.from_imod5_data(imod5_data) - # Test if regridded to smallest grid resolution + # Test if regridded to BND resolution, which is 25 by 25 m. TOP and BOT were + # 100 by 100 m. assert dis["top"].dx == 25.0 assert dis["top"].dy == -25.0 assert (dis.dataset.coords["x"][1] - dis.dataset.coords["x"][0]) == 25.0 @@ -354,6 +355,37 @@ def test_from_imod5_data__grid_extent(imod5_dataset): assert dis.dataset.coords["x"].max() == 199287.5 +def test_from_imod5_data_regridding__target_grid(imod5_dataset): + imod5_data = imod5_dataset[0] + + xmin = 195000.0 + xmax = 199000.0 + ymin = 361000.0 + ymax = 365000.0 + dx = 100.0 + dy = 100.0 + + target_grid = imod.util.empty_2d( + dx=dx, xmin=xmin, xmax=xmax, dy=dy, ymin=ymin, ymax=ymax + ) + dis = imod.mf6.StructuredDiscretization.from_imod5_data( + imod5_data, target_grid=target_grid, validate=False + ) + + # Test if regridded to BND resolution, which is 25 by 25 m. TOP and BOT were + # 100 by 100 m. + assert dis["top"].dx == dx + assert dis["top"].dy == -dy + assert (dis.dataset.coords["x"][1] - dis.dataset.coords["x"][0]) == dx + assert (dis.dataset.coords["y"][1] - dis.dataset.coords["y"][0]) == -dy + + # Test extent + assert dis.dataset.coords["y"].min() == ymin + dy / 2 + assert dis.dataset.coords["y"].max() == ymax - dy / 2 + assert dis.dataset.coords["x"].min() == xmin + dx / 2 + assert dis.dataset.coords["x"].max() == xmax - dx / 2 + + def test_from_imod5_data__write(imod5_dataset, tmp_path): directory = tmp_path / "dis_griddata" directory.mkdir() diff --git a/imod/tests/test_mf6/test_mf6_simulation.py b/imod/tests/test_mf6/test_mf6_simulation.py index 888daaa0e..ef7386a01 100644 --- a/imod/tests/test_mf6/test_mf6_simulation.py +++ b/imod/tests/test_mf6/test_mf6_simulation.py @@ -721,6 +721,7 @@ def test_import_from_imod5__well_steady_state(imod5_dataset): @pytest.mark.unittest_jit def test_import_from_imod5__nonstandard_regridding(imod5_dataset, tmp_path): + # Arrange imod5_data = imod5_dataset[0] period_data = imod5_dataset[1] @@ -730,13 +731,24 @@ def test_import_from_imod5__nonstandard_regridding(imod5_dataset, tmp_path): regridding_option["sto"] = StorageCoefficientRegridMethod() times = pd.date_range(start="1/1/2018", end="12/1/2018", freq="ME") + xmin = 195000.0 + xmax = 199000.0 + ymin = 361000.0 + ymax = 365000.0 + dx = 200.0 + dy = 200.0 + + target_grid = imod.util.empty_2d( + dx=dx, xmin=xmin, xmax=xmax, dy=dy, ymin=ymin, ymax=ymax + ) + + # Act simulation = Modflow6Simulation.from_imod5_data( imod5_data, period_data, times, - SimulationAllocationOptions, - SimulationDistributingOptions, - regridding_option, + regridder_types=regridding_option, + target_grid=target_grid, ) simulation["imported_model"]["oc"] = OutputControl( save_head="last", save_budget="last" @@ -750,8 +762,18 @@ def test_import_from_imod5__nonstandard_regridding(imod5_dataset, tmp_path): # Align NoData to domain idomain = simulation["imported_model"].domain simulation.mask_all_models(idomain) + + # Assert # write and validate the simulation. simulation.write(tmp_path, binary=False, validate=True) + # Check that storage package regridded to target_grid + coords = simulation["imported_model"]["sto"].dataset.coords + assert coords["x"][1] - coords["x"][0] == dx + assert coords["y"][1] - coords["y"][0] == -dy + assert coords["x"].min() == xmin + dx / 2 + assert coords["x"].max() == xmax - dx / 2 + assert coords["y"].min() == ymin + dy / 2 + assert coords["y"].max() == ymax - dy / 2 @pytest.mark.unittest_jit From f608df712c10c894a80c86231b4b1d6db4d484ee Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 22 Apr 2026 11:59:02 +0200 Subject: [PATCH 5/6] Expand docstring with examples --- imod/mf6/simulation.py | 31 +++++++++++++++++++++++++++---- 1 file changed, 27 insertions(+), 4 deletions(-) diff --git a/imod/mf6/simulation.py b/imod/mf6/simulation.py index e2ff57b05..5749bb223 100644 --- a/imod/mf6/simulation.py +++ b/imod/mf6/simulation.py @@ -1817,16 +1817,39 @@ def from_imod5_data( >>> times = [np.datetime("2001-01-01"), np.datetime("2002-01-01"), np.datetime("2003-01-01")] >>> mf6_sim = imod.mf6.Modflow6Simulation.from_imod5_data(imod5_data, period_data, times) - Allocate rivers differently: + You can override solver settings if needed after importing: + + >>> mf6_sim["imported_model"]["ims"] = SolutionPresetSimple() + + To allocate rivers differently: >>> from imod.prepare.topsystem import SimulationAllocationOptions, ALLOCATION_OPTION >>> allocation_options = SimulationAllocationOptions() >>> allocation_options.riv = ALLOCATION_OPTION.at_elevation - >>> mf6_sim = imod.mf6.Modflow6Simulation.from_imod5_data(imod5_data, period_data, times, allocation_options) + >>> mf6_sim = imod.mf6.Modflow6Simulation.from_imod5_data( + >>> imod5_data, period_data, times, allocation_options=allocation_options + >>> ) - You can override solver settings if needed after importing: + To regrid the model to a specific grid upon import: + + >>> target_grid = imod.util.empty_2d( + >>> dx=100.0, xmin=195000.0, xmax=199000.0, dy=100.0, ymin=361000.0, ymax=365000.0 + >>> ) + >>> mf6_sim = imod.mf6.Modflow6Simulation.from_imod5_data( + >>> imod5_data, period_data, times, target_grid=target_grid + >>> ) + + To set regridding methods for specific packages upon import: + + >>> from imod.util.regrid import RegridderType + >>> from imod.mf6.regrid import NodePropertyFlowRegridMethod + >>> regridder_types = { + >>> "npf": NodePropertyFlowRegridMethod(k=(RegridderType.OVERLAP, "mean")), + >>> } + >>> mf6_sim = imod.mf6.Modflow6Simulation.from_imod5_data( + >>> imod5_data, period_data, times, regridder_types=regridder_types, target_grid=target_grid + >>> ) - >>> mf6_sim["imported_model"]["ims"] = SolutionPresetComplex() """ if allocation_options is None: From 7e82e8d7e27b3ae417ea53b163d62f32cb7a0022 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 22 Apr 2026 14:45:32 +0200 Subject: [PATCH 6/6] Elaborate further on allocation --- imod/mf6/simulation.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/imod/mf6/simulation.py b/imod/mf6/simulation.py index 5749bb223..2589fad14 100644 --- a/imod/mf6/simulation.py +++ b/imod/mf6/simulation.py @@ -1821,7 +1821,10 @@ def from_imod5_data( >>> mf6_sim["imported_model"]["ims"] = SolutionPresetSimple() - To allocate rivers differently: + To allocate rivers to model layers differently than the default, you can + set the allocation options for the river package before importing. + :doc:`For more information on topsystem allocation see the user guide. + ` >>> from imod.prepare.topsystem import SimulationAllocationOptions, ALLOCATION_OPTION >>> allocation_options = SimulationAllocationOptions()