Issue #1311 sprinkling from imod5 data (#1312)

Fixes #1311

# Description

Adds the following:

- ``imod.mf6.LayeredWell.from_imod5_cap_data``, to construct
LayeredWells for MODLOW6 from the CAP package in iMOD5 Data. Throws an
error if data is provided as point data (IPF, which is read as pandas
DataFrame). This has rate 0.0, actual rates will be inserted by in the
coupling scheme of iMOD Coupler.
- ``imod.mf6.GroundwaterFlowModel.from_imod5_data`` adds Sprinkling well
if CAP package present. This does not affect MODFLOW6 model if not
coupled to MetaSWAP: The wells have rates 0.0, so will not affect the
simulation results.
- ``Sprinkling.from_imod5_data`` to construct a Sprinkling package from
iMOD5 data.

# Checklist
<!---
Before requesting review, please go through this checklist:
-->

- [x] Links to correct issue
- [x] Update changelog, if changes affect users
- [x] PR title starts with ``Issue #nr``, e.g. ``Issue #737``
- [x] Unit tests were added
- [ ] **If feature added**: Added/extended example

docs/api/changelog.rst

@@ -14,6 +14,9 @@ Added
 
 - :class:`imod.msw.MeteoGridCopy` to copy existing `mete_grid.inp` files, so
   ASCII grids in large existing meteo databases do not have to be read.
+- :meth:`imod.mf6.LayeredWell.from_imod5_cap_data` to construct a
+  :class:`imod.mf6.LayeredWell` package from iMOD5 data in the CAP package (for
+  MetaSWAP). Currently only griddata (IDF) is supported.
 - :meth:`imod.mf6.Recharge.from_imod5_cap_data` to construct a recharge package
   for coupling a MODFLOW6 model to MetaSWAP.
 
@@ -35,6 +38,9 @@ Changed
   ``downward_resistance`` now require a ``subunit`` coordinate.
 - Variables ``max_abstraction_groundwater`` and ``max_abstraction_surfacewater``
   in :class:`imod.msw.Sprinkling` now needs to have a subunit coordinate.
+- If ``"cap"`` package present in ``imod5_data``,
+  :meth:`imod.mf6.GroundwaterFlowModel.from_imod5_data` now automatically adds a
+  well for metaswap sprinkling named ``"msw-sprinkling"``
 
 
 [0.18.1] - 2024-11-20

docs/api/mf6.rst

@@ -97,6 +97,7 @@ Flow Packages
     HorizontalFlowBarrierResistance.to_mf6_pkg
     LayeredWell
     LayeredWell.from_imod5_data
+    LayeredWell.from_imod5_cap_data
     LayeredWell.mask
     LayeredWell.regrid_like
     LayeredWell.to_mf6_pkg

imod/mf6/model_gwf.py

@@ -303,9 +303,9 @@ def from_imod5_data(
             )
 
         # now import the non-singleton packages'
+        imod5_keys = list(imod5_data.keys())
 
         # import wells
-        imod5_keys = list(imod5_data.keys())
         wel_keys = [key for key in imod5_keys if key[0:3] == "wel"]
         for wel_key in wel_keys:
             wel_key_truncated = wel_key[:16]
@@ -331,10 +331,10 @@ 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)  # type: ignore
 
         # import ghb's
-        imod5_keys = list(imod5_data.keys())
         ghb_keys = [key for key in imod5_keys if key[0:3] == "ghb"]
         for ghb_key in ghb_keys:
             ghb_pkg = GeneralHeadBoundary.from_imod5_data(

imod/mf6/utilities/imod5_converter.py

@@ -1,8 +1,10 @@
-from typing import Union
+from typing import Union, cast
 
 import numpy as np
+import pandas as pd
 import xarray as xr
 
+from imod.typing import GridDataDict, Imod5DataDict, SelSettingsType
 from imod.typing.grid import full_like
 
 
@@ -48,3 +50,72 @@ def fill_missing_layers(
     """
     layer = full.coords["layer"]
     return source.reindex(layer=layer, fill_value=fillvalue)
+
+
+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."
+    )
+    return {}
+
+
+def _well_from_imod5_cap_grid_data(cap_data: GridDataDict) -> dict[str, np.ndarray]:
+    drop_layer_kwargs: SelSettingsType = {
+        "layer": 0,
+        "drop": True,
+        "missing_dims": "ignore",
+    }
+    type = cap_data["artificial_recharge"].isel(**drop_layer_kwargs).compute()
+    layer = (
+        cap_data["artificial_recharge_layer"]
+        .isel(**drop_layer_kwargs)
+        .astype(int)
+        .compute()
+    )
+
+    from_groundwater = (type != 0).to_numpy()
+    coords = type.coords
+    x_grid, y_grid = np.meshgrid(coords["x"].to_numpy(), coords["y"].to_numpy())
+
+    data = {}
+    data["layer"] = layer.data[from_groundwater]
+    data["y"] = y_grid[from_groundwater]
+    data["x"] = x_grid[from_groundwater]
+    data["rate"] = np.zeros_like(data["x"])
+
+    return data
+
+
+def well_from_imod5_cap_data(imod5_data: Imod5DataDict) -> 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
+    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.
+    """
+    cap_data = cast(GridDataDict, imod5_data["cap"])
+
+    if isinstance(cap_data["artificial_recharge_layer"], pd.DataFrame):
+        return _well_from_imod5_cap_point_data(cap_data)
+    else:
+        return _well_from_imod5_cap_grid_data(cap_data)

imod/mf6/wel.py

@@ -30,6 +30,7 @@
 from imod.mf6.package import Package
 from imod.mf6.utilities.dataset import remove_inactive
 from imod.mf6.utilities.grid import broadcast_to_full_domain
+from imod.mf6.utilities.imod5_converter import well_from_imod5_cap_data
 from imod.mf6.validation import validation_pkg_error_message
 from imod.mf6.validation_context import ValidationContext
 from imod.mf6.write_context import WriteContext
@@ -44,7 +45,7 @@
     ValidationError,
 )
 from imod.select.points import points_indices, points_values
-from imod.typing import GridDataArray
+from imod.typing import GridDataArray, Imod5DataDict
 from imod.typing.grid import is_spatial_grid, ones_like
 from imod.util.expand_repetitions import resample_timeseries
 from imod.util.structured import values_within_range
@@ -1298,6 +1299,46 @@ def _validate_imod5_depth_information(
             logger.log(loglevel=LogLevel.ERROR, message=log_msg, additional_depth=2)
             raise ValueError(log_msg)
 
+    @classmethod
+    def from_imod5_cap_data(cls, imod5_data: Imod5DataDict):
+        """
+        Create LayeredWell from imod5_data in "cap" package. 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
+        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.
+
+        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`.
+        """
+        data = well_from_imod5_cap_data(imod5_data)
+        return cls(**data)  # type: ignore
+
 
 class WellDisStructured(DisStructuredBoundaryCondition):
     """

imod/msw/sprinkling.py

@@ -1,4 +1,4 @@
-from typing import TextIO
+from typing import TextIO, cast
 
 import numpy as np
 import pandas as pd
@@ -10,7 +10,9 @@
 from imod.msw.fixed_format import VariableMetaData
 from imod.msw.pkgbase import MetaSwapPackage
 from imod.msw.regrid.regrid_schemes import SprinklingRegridMethod
-from imod.typing import IntArray
+from imod.msw.utilities.common import concat_imod5
+from imod.typing import GridDataDict, Imod5DataDict, IntArray, SelSettingsType
+from imod.typing.grid import zeros_like
 
 
 def _ravel_per_subunit(da: xr.DataArray) -> np.ndarray:
@@ -20,6 +22,49 @@ 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:
+    raise NotImplementedError(
+        "Assigning sprinkling wells with an IPF file is not supported, please specify them as IDF."
+    )
+    return {}
+
+
+def _sprinkling_data_from_imod5_grid(cap_data: GridDataDict) -> GridDataDict:
+    drop_layer_kwargs: SelSettingsType = {
+        "layer": 0,
+        "drop": True,
+        "missing_dims": "ignore",
+    }
+    type = cap_data["artificial_recharge"].isel(**drop_layer_kwargs)
+    capacity = cap_data["artificial_recharge_capacity"].isel(**drop_layer_kwargs)
+
+    from_groundwater = type == 1
+    from_surfacewater = type == 2
+    is_active = type != 0
+
+    zero_where_active = zeros_like(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
+    )
+    max_abstraction_surfacewater_rural = zero_where_active.where(
+        ~from_surfacewater, capacity
+    )
+
+    # 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
+
+
 class Sprinkling(MetaSwapPackage, IRegridPackage):
     """
     This contains the sprinkling capacities of links between SVAT units and
@@ -122,3 +167,51 @@ def _render(
         self._check_range(dataframe)
 
         return self.write_dataframe_fixed_width(file, dataframe)
+
+    @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.
+
+        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`.
+
+        Returns
+        -------
+        Sprinkling package
+        """
+        cap_data = cast(GridDataDict, 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)
+
+        return cls(**data)

imod/tests/conftest.py

@@ -22,6 +22,10 @@
 )
 from .fixtures.flow_example_fixture import imodflow_model
 from .fixtures.flow_transport_simulation_fixture import flow_transport_simulation
+from .fixtures.imod5_cap_data import (
+    cap_data_sprinkling_grid,
+    cap_data_sprinkling_points,
+)
 from .fixtures.imod5_well_data import (
     well_duplication_import_prj,
     well_mixed_ipfs,

imod/tests/fixtures/imod5_cap_data.py

@@ -0,0 +1,60 @@
+import numpy as np
+import pandas as pd
+import pytest
+import xarray as xr
+
+from imod.typing import Imod5DataDict
+
+
+def zeros_grid():
+    x = [1.0, 2.0, 3.0]
+    y = [3.0, 2.0, 1.0]
+    dx = 1.0
+    dy = -1.0
+
+    coords = {"x": x, "y": y, "dx": dx, "dy": dy}
+    shape = (len(y), len(x))
+    data = np.zeros(shape)
+
+    return xr.DataArray(data, coords=coords, dims=("y", "x"))
+
+
+@pytest.fixture(scope="function")
+def cap_data_sprinkling_grid() -> Imod5DataDict:
+    type = zeros_grid()
+    type[:, 1] = 1
+    type[:, 2] = 2
+    layer = xr.ones_like(type)
+    layer[:, 1] = 2
+
+    cap_data = {
+        "artificial_recharge": type,
+        "artificial_recharge_layer": layer,
+        "artificial_recharge_capacity": xr.DataArray(25.0),
+    }
+
+    return {"cap": cap_data}
+
+
+@pytest.fixture(scope="function")
+def cap_data_sprinkling_points() -> Imod5DataDict:
+    type = zeros_grid()
+    type[:, 1] = 3000
+    type[:, 2] = 4000
+
+    data = {
+        "id": [3000, 4000],
+        "layer": [2, 3],
+        "capacity": [15.0, 30.0],
+        "y": [1.0, 2.0],
+        "x": [1.0, 2.0],
+    }
+
+    layer = pd.DataFrame(data=data)
+    cap_data = {
+        "artificial_recharge": type,
+        "artificial_recharge_layer": layer,
+        "artificial_recharge_capacity": xr.DataArray(25.0),
+    }
+
+    return {"cap": cap_data}