Issue #1775 improve target grid

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
 -----

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)

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,14 +177,14 @@ 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
         :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::
 
@@ -202,6 +201,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 +216,10 @@ 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()

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,6 +245,9 @@ 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
         -------
@@ -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

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
         -------
@@ -1813,16 +1817,42 @@ 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 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.
+        </user-guide/09-topsystem>`
 
         >>> 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:
@@ -1848,6 +1878,7 @@ def from_imod5_data(
             allocation_options,
             distributing_options,
             regridder_types,
+            target_grid,
         )
         simulation["imported_model"] = gwf_model
 

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()

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