Issue #1823 regrid imod5 cap data

docs/api/changelog.rst

@@ -64,6 +64,10 @@ Fixed
 - :meth:`imod.mf6.NodePropertyFlow.regrid_like` now regrids ``k33`` using the
   correct method, namely ``mean`` instead of ``harmonic_mean``. As this is the
   appropriate method for horizontal regridding of ``k33``.
+- :meth:`imod.msw.MetaSwapModel.from_imod5_data`,
+  :meth:`imod.mf6.Recharge.from_imod5_cap_data`,
+  :meth:`imod.mf6.LayeredWell.from_imod5_cap_data` now regrids the iMOD5 CAP
+  data to the MODFLOW6 target discretization.
 
 Changed
 ~~~~~~~

imod/common/utilities/regrid.py

@@ -20,13 +20,14 @@
 from imod.common.utilities.dataclass_type import DataclassType, EmptyRegridMethod
 from imod.common.utilities.dtype import is_integer
 from imod.common.utilities.value_filters import is_valid
+from imod.typing import Imod5DataDict
 from imod.typing.grid import (
     GridDataArray,
     GridDataset,
     is_unstructured,
     ones_like,
 )
-from imod.util.dims import enforced_dim_order
+from imod.util.dims import drop_layer_dim_cap_data, enforced_dim_order
 from imod.util.regrid import (
     RegridderType,
     RegridderWeightsCache,
@@ -77,7 +78,7 @@ def _regrid_array(
 
     # the dataarray might be a scalar. If it is, then it does not need regridding.
     if scalar_da:
-        return da.values[()]
+        return da.squeeze(drop=True)
 
     if isinstance(da, xr.DataArray):
         coords = da.coords
@@ -459,3 +460,32 @@ def _get_regridding_domain(
     ).astype(int)
 
     return new_idomain
+
+
+def regrid_imod5_cap_data(
+    imod5_data: Imod5DataDict,
+    target_dis: IRegridPackage,
+    regridder_types: DataclassType,
+    regrid_cache: RegridderWeightsCache,
+) -> Imod5DataDict:
+    """
+    Regrid iMOD5 CAP data to consistent grid. This is necessary to be able to
+    use iMOD5 data in MetaSWAP, as the grid of the iMOD5 CAP data is not
+    necessarily the same as the grid of the target MODFLOW6 discretization. The
+    regridding process ensures consistency between the iMOD5 CAP data and the
+    target MODFLOW6 grid. Used in both ``imod.msw.MetaSwapModel.from_imod5_data``
+    and ``imod.mf6.Recharge.from_imod5_cap_data``.
+    """
+    # Drop layer coords
+    imod5_cap_no_layer = drop_layer_dim_cap_data(imod5_data)
+    target_grid = target_dis.dataset["idomain"].isel(layer=0, drop=True)
+    # Regrid the input data
+    cap_data_regridded = _regrid_package_data(
+        imod5_cap_no_layer["cap"], target_grid, regridder_types, regrid_cache
+    )
+    extra_paths = imod5_data["extra"]["paths"]
+    imod5_regridded: Imod5DataDict = {
+        "cap": cap_data_regridded,
+        "extra": {"paths": extra_paths},
+    }
+    return imod5_regridded

imod/formats/prj/prj.py

@@ -98,7 +98,7 @@
     "rural_infiltration_capacity",
     "perched_water_table_level",
     "soil_moisture_fraction",
-    "conductivitiy_factor",
+    "conductivity_factor",
     "plot_number",
     "steering_location",
     "plot_drainage_level",

imod/mf6/model_gwf.py

@@ -45,7 +45,7 @@
     SimulationDistributingOptions,
 )
 from imod.schemata import TypeSchema
-from imod.typing import GridDataArray, StressPeriodTimesType
+from imod.typing import GridDataArray, Imod5DataDict, StressPeriodTimesType
 from imod.typing.grid import zeros_like
 from imod.util.regrid import RegridderWeightsCache
 
@@ -354,8 +354,13 @@ def from_imod5_data(
             )
 
         if "cap" in imod5_keys:
-            result["msw-sprinkling"] = LayeredWell.from_imod5_cap_data(imod5_data)  # type: ignore
-            result["msw-rch"] = Recharge.from_imod5_cap_data(imod5_data, dis_pkg)  # type: ignore
+            imod5_cap_data_dict = cast(Imod5DataDict, imod5_data)
+            result["msw-sprinkling"] = LayeredWell.from_imod5_cap_data(
+                imod5_cap_data_dict, dis_pkg, regrid_cache=regrid_cache
+            )
+            result["msw-rch"] = Recharge.from_imod5_cap_data(
+                imod5_cap_data_dict, dis_pkg, regrid_cache=regrid_cache
+            )
 
         # import ghb's
         ghb_keys = [key for key in imod5_keys if key[0:3] == "ghb"]

imod/mf6/rch.py

@@ -6,10 +6,14 @@
 
 from imod.common.interfaces.iregridpackage import IRegridPackage
 from imod.common.utilities.dataclass_type import DataclassType
+from imod.common.utilities.regrid import regrid_imod5_cap_data
 from imod.logging import init_log_decorator
 from imod.mf6.aggregate.aggregate_schemes import RechargeAggregationMethod
 from imod.mf6.dis import StructuredDiscretization, VerticesDiscretization
-from imod.mf6.regrid.regrid_schemes import RechargeRegridMethod
+from imod.mf6.regrid.regrid_schemes import (
+    CapDataRechargeRegridMethod,
+    RechargeRegridMethod,
+)
 from imod.mf6.topsystem import TopSystemBoundaryCondition
 from imod.mf6.utilities.imod5_converter import (
     convert_unit_rch_rate,
@@ -39,7 +43,6 @@
     enforce_dim_order,
     is_planar_grid,
 )
-from imod.util.dims import drop_layer_dim_cap_data
 from imod.util.regrid import RegridderWeightsCache
 
 
@@ -286,14 +289,18 @@ def from_imod5_cap_data(
         cls,
         imod5_data: Imod5DataDict,
         target_dis: StructuredDiscretization,
+        regridder_types: CapDataRechargeRegridMethod = CapDataRechargeRegridMethod(),
+        regrid_cache: RegridderWeightsCache = RegridderWeightsCache(),
     ) -> "Recharge":
         """
         Construct an rch-package from iMOD5 data in the CAP package, loaded with
         the :func:`imod.formats.prj.open_projectfile_data` function. Package is
         used to couple MODFLOW6 to MetaSWAP models. Active cells will have a
         recharge rate of 0.0.
         """
-        cap_data = drop_layer_dim_cap_data(imod5_data)["cap"]
+        cap_data = regrid_imod5_cap_data(
+            imod5_data, target_dis, regridder_types, regrid_cache
+        )["cap"]
 
         msw_area = get_cell_area_from_imod5_data(cap_data)
         msw_active = is_msw_active_cell(target_dis, cap_data, msw_area)

imod/mf6/regrid/regrid_schemes.py

@@ -450,3 +450,28 @@ class StorageCoefficientRegridMethod(DataclassType):
     convertible: RegridVarType = (RegridderType.OVERLAP, "mode")
     storage_coefficient: RegridVarType = (RegridderType.OVERLAP, "mean")
     specific_yield: RegridVarType = (RegridderType.OVERLAP, "mean")
+
+
+@dataclass(config=_CONFIG)
+class CapDataRechargeRegridMethod(DataclassType):
+    """
+    Object containing regridder methods for CAP data for the
+    :class:`imod.mf6.Recharge.from_imod5_cap_data` method. This contains regridder
+    methods for only the relevant CAP variables for the recharge package.
+    """
+
+    boundary: RegridVarType = (RegridderType.OVERLAP, "mode")
+    wetted_area: RegridVarType = (RegridderType.RELATIVEOVERLAP, "conductance")
+    urban_area: RegridVarType = (RegridderType.RELATIVEOVERLAP, "conductance")
+
+
+@dataclass(config=_CONFIG)
+class CapDataWellRegridMethod(DataclassType):
+    """
+    Object containing regridder methods for CAP data for the
+    :class:`imod.mf6.LayeredWell.from_imod5_cap_data` method. This contains
+    regridder methods for only the relevant CAP variables for the well package.
+    """
+
+    artificial_recharge: RegridVarType = (RegridderType.OVERLAP, "mean")
+    artificial_recharge_layer: RegridVarType = (RegridderType.OVERLAP, "mode")

imod/mf6/riv.py

@@ -141,7 +141,7 @@ def _create_drain_from_leftover_riv_imod5_data(
     else:
         drain_leftover_data = infiltration_drn_data
 
-    return Drainage(**drain_leftover_data)  # type: ignore[arg-type]
+    return Drainage(**drain_leftover_data)  # type: ignore
 
 
 class River(TopSystemBoundaryCondition, IRegridPackage):

imod/mf6/utilities/imod5_converter.py

@@ -4,12 +4,12 @@
 import pandas as pd
 import xarray as xr
 
+from imod.common.interfaces.iregridpackage import IRegridPackage
 from imod.common.utilities.dataclass_type import DataclassType
-from imod.common.utilities.regrid import _regrid_package_data
+from imod.common.utilities.regrid import _regrid_package_data, regrid_imod5_cap_data
 from imod.mf6.package import Package
 from imod.typing import GridDataArray, GridDataDict, Imod5DataDict
 from imod.typing.grid import full_like
-from imod.util.dims import drop_layer_dim_cap_data
 from imod.util.regrid import RegridderWeightsCache
 
 
@@ -83,7 +83,12 @@ def _well_from_imod5_cap_grid_data(cap_data: GridDataDict) -> dict[str, np.ndarr
     return data
 
 
-def well_from_imod5_cap_data(imod5_data: Imod5DataDict) -> dict[str, np.ndarray]:
+def well_from_imod5_cap_data(
+    imod5_data: Imod5DataDict,
+    target_dis: IRegridPackage,
+    regridder_types: DataclassType,
+    regrid_cache: RegridderWeightsCache,
+) -> dict[str, np.ndarray]:
     """
     Abstraction data for sprinkling is defined in iMOD5 either with grids (IDF)
     or points (IPF) combined with a grid. Depending on the type, the function does
@@ -110,12 +115,16 @@ def well_from_imod5_cap_data(imod5_data: Imod5DataDict) -> dict[str, np.ndarray]
         capacity is already defined in the point data. This is an ``n:1``
         mapping: multiple grid cells can map to one well.
     """
-    cap_data = drop_layer_dim_cap_data(imod5_data)["cap"]
+    cap_data = imod5_data["cap"]
+    has_ipf_well = isinstance(cap_data["artificial_recharge_layer"], pd.DataFrame)
 
-    if isinstance(cap_data["artificial_recharge_layer"], pd.DataFrame):
+    if has_ipf_well:
         return _well_from_imod5_cap_point_data(cap_data)
     else:
-        return _well_from_imod5_cap_grid_data(cap_data)
+        cap_data_regridded = regrid_imod5_cap_data(
+            imod5_data, target_dis, regridder_types, regrid_cache
+        )["cap"]
+        return _well_from_imod5_cap_grid_data(cap_data_regridded)
 
 
 def regrid_imod5_pkg_data(

imod/mf6/wel.py

@@ -28,6 +28,7 @@
 )
 from imod.mf6.mf6_wel_adapter import Mf6Wel, concat_indices_to_cellid
 from imod.mf6.package import Package
+from imod.mf6.regrid.regrid_schemes import CapDataWellRegridMethod
 from imod.mf6.utilities.dataset import remove_inactive
 from imod.mf6.utilities.imod5_converter import well_from_imod5_cap_data
 from imod.mf6.validation_settings import ValidationSettings
@@ -45,6 +46,7 @@
 from imod.typing import GridDataArray, Imod5DataDict, StressPeriodTimesType
 from imod.typing.grid import is_spatial_grid, ones_like
 from imod.util.expand_repetitions import average_timeseries, resample_timeseries
+from imod.util.regrid import RegridderWeightsCache
 from imod.util.structured import values_within_range
 
 ABSTRACT_METH_ERROR_MSG = "Method in abstract base class called"
@@ -1443,7 +1445,13 @@ def _validate_imod5_depth_information(
             raise ValueError(log_msg)
 
     @classmethod
-    def from_imod5_cap_data(cls, imod5_data: Imod5DataDict):
+    def from_imod5_cap_data(
+        cls,
+        imod5_data: Imod5DataDict,
+        target_dis: StructuredDiscretization,
+        regridder_types: CapDataWellRegridMethod = CapDataWellRegridMethod(),
+        regrid_cache: RegridderWeightsCache = RegridderWeightsCache(),
+    ):
         """
         Create LayeredWell from imod5_data in "cap" package. Abstraction data
         for sprinkling is defined in iMOD5 either with grids (IDF) or points
@@ -1481,7 +1489,9 @@ def from_imod5_cap_data(cls, imod5_data: Imod5DataDict):
             belongs, as returned by
             :func:`imod.formats.prj.open_projectfile_data`.
         """
-        data = well_from_imod5_cap_data(imod5_data)
+        data = well_from_imod5_cap_data(
+            imod5_data, target_dis, regridder_types, regrid_cache
+        )
         return cls(**data)  # type: ignore
 
     @standard_log_decorator()

imod/msw/model.py

@@ -13,6 +13,7 @@
 from imod.common.constants import MaskValues
 from imod.common.utilities.clip import clip_by_grid
 from imod.common.utilities.partitioninfo import create_partition_info
+from imod.common.utilities.regrid import regrid_imod5_cap_data
 from imod.common.utilities.value_filters import enforce_scalar
 from imod.common.utilities.version import prepend_content_with_version_info
 from imod.mf6.dis import StructuredDiscretization
@@ -38,16 +39,18 @@
 from imod.msw.output_control import TimeOutputControl
 from imod.msw.pkgbase import MetaSwapPackage
 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.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
+from imod.msw.utilities.imod5_converter import (
+    has_active_scaling_factor,
+)
 from imod.msw.utilities.mask import mask_and_broadcast_cap_data
 from imod.msw.utilities.parse import read_para_sim
 from imod.msw.vegetation import AnnualCropFactors
 from imod.typing import GridDataArray, Imod5DataDict
-from imod.util.dims import drop_layer_dim_cap_data
 from imod.util.regrid import RegridderWeightsCache
 from imod.util.time import to_datetime_internal
 
@@ -605,7 +608,9 @@ def from_imod5_data(
         imod5_data: Imod5DataDict,
         target_dis: StructuredDiscretization,
         times: list[datetime],
-    ):
+        regridder_types: CapDataRegridMethod = CapDataRegridMethod(),
+        regrid_cache: RegridderWeightsCache = RegridderWeightsCache(),
+    ) -> "MetaSwapModel":
         """
         Construct a MetaSWAP model from iMOD5 data in the CAP package, loaded
         with the :func:`imod.formats.prj.open_projectfile_data` function.
@@ -616,29 +621,43 @@ def from_imod5_data(
             Dictionary with iMOD5 data. This can be constructed from the
             :func:`imod.formats.prj.open_projectfile_data` method.
         target_dis: imod.mf6.StructuredDiscretization
-            Target discretization, cells where MODLOW6 is inactive will be
-            inactive in MetaSWAP as well.
+            Target discretization, iMOD5 CAP data will be regridded to this
+            discretization. Cells where MODLOW6 is inactive will be inactive in
+            MetaSWAP as well.
         times: list[datetime]
             List of datetimes, will be used to set the output control times.
             Is also used to infer the starttime of the simulation.
+        regridder_types: CapDataRegridMethod, default CapDataRegridMethod()
+            Custom regrid method for CAP data.
+        regrid_cache: RegridderWeightsCache, default RegridderWeightsCache()
+            Cache for regridder weights, can be used to speed up regridding if the
+            same regridders are used multiple times.
 
         Returns
         -------
         MetaSwapModel
             MetaSWAP model imported from imod5 data.
         """
+        # Path and settings management
         extra_paths = imod5_data["extra"]["paths"]
         path_to_parasim = find_in_file_list("para_sim.inp", extra_paths)
         parasim_settings = read_para_sim(path_to_parasim)
         unsa_svat_path = cast(str, parasim_settings["unsa_svat_path"])
-        # Drop layer coord
-        imod5_cap_no_layer = drop_layer_dim_cap_data(imod5_data)
+        # Regrid iMOD5 CAP data to target discretization.
+        imod5_regridded = regrid_imod5_cap_data(
+            imod5_data, target_dis, regridder_types, regrid_cache
+        )
+        # Test with regridded data instead of masked, as masking broadcasts
+        # scalars to grids, which causes the is_scalar check in
+        # has_active_scaling_factor to always return False.
+        active_scaling_factor = has_active_scaling_factor(imod5_regridded["cap"])
+        # Setup model
         model = cls(unsa_svat_path, parasim_settings)
         model["grid"], msw_active = GridData.from_imod5_data(
-            imod5_cap_no_layer, target_dis
+            imod5_regridded, target_dis
         )
         cap_data_masked = mask_and_broadcast_cap_data(
-            imod5_cap_no_layer["cap"], msw_active
+            imod5_regridded["cap"], msw_active
         )
         imod5_masked: Imod5DataDict = {
             "cap": cap_data_masked,
@@ -650,7 +669,7 @@ def from_imod5_data(
         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)
-        if has_active_scaling_factor(imod5_cap_no_layer["cap"]):
+        if active_scaling_factor:
             model["scaling_factor"] = ScalingFactors.from_imod5_data(imod5_masked)
         area = model["grid"]["area"].isel(subunit=0, drop=True)
         model["idf_mapping"] = IdfMapping(area, MaskValues.msw_default)