Issue #1260 from imod5 data metaswap (#1335)

Fixes #1260

# Description
Adds ``MetaSwapModel.from_imod5_data`` 

# 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

Author: JoerivanEngelen

docs/api/changelog.rst

@@ -12,6 +12,21 @@ The format is based on `Keep a Changelog`_, and this project adheres to
 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.
+- :class:`imod.msw.CopyFiles` to copy settings and lookup tables in existing
+  ``.inp`` files.
+- :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.
+- :meth:`imod.msw.MetaSwapModel.from_imod5_data` to construct a MetaSWAP model
+  from data in an iMOD5 projectfile.
+- :meth:`imod.msw.MetaSwapModel.write` has a ``validate`` argument, which can be
+  used to turn off validation upon writing, use at your own risk!
+- :class:`imod.msw.MetaSwapModel` got ``settings`` argument to set simulation
+  settings.
 - :func:`imod.data.tutorial_03` to load data for the iMOD Documentation
   tutorial.
 
@@ -35,8 +50,13 @@ Fixed
 Changed
 ~~~~~~~
 
+- :class:`imod.msw.Infiltration`'s variables ``upward_resistance`` and
+  ``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

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

docs/api/msw.rst

@@ -7,10 +7,15 @@ MetaSWAP
     :toctree: generated/msw
 
     GridData
+    GridData.from_imod5_data
     Infiltration
+    Infiltration.from_imod5_data
     Ponding
+    Ponding.from_imod5_data
     ScalingFactors
+    ScalingFactors.from_imod5_data
     Sprinkling
+    Sprinkling.from_imod5_data
 
     IdfMapping
     TimeOutputControl
@@ -25,9 +30,15 @@ MetaSWAP
     AnnualCropFactors
 
     MeteoGrid
+    MeteoGridCopy
+    MeteoGridCopy.from_imod5_data
     EvapotranspirationMapping
+    EvapotranspirationMapping.from_imod5_data
     PrecipitationMapping
+    PrecipitationMapping.from_imod5_data
 
     CouplerMapping
 
     MetaSwapModel
+    MetaSwapModel.write
+    MetaSwapModel.from_imod5_data

imod/mf6/model_gwf.py

@@ -218,9 +218,11 @@ def from_imod5_data(
         times: list[datetime]
             Time discretization of the simulation. These times are used for the
             following:
-               * Times of wells with associated timeseries are resampled to these times
-               * Start- and end times in the list are used to repeat the stresses
-                 of periodic data (e.g. river stages in iMOD5 for "summer", "winter")
+
+                * Times of wells with associated timeseries are resampled to these times
+                * Start- and end times in the list are used to repeat the stresses
+                  of periodic data (e.g. river stages in iMOD5 for "summer", "winter")
+
         allocation_options: SimulationAllocationOptions
             object containing the allocation options per package type.
             If you want a package to have a different allocation option,
@@ -303,9 +305,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]
@@ -330,8 +332,11 @@ def from_imod5_data(
                     wel_key, imod5_data, times
                 )
 
+        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/rch.py

@@ -1,7 +1,8 @@
 from copy import deepcopy
-from typing import Optional
+from typing import Optional, cast
 
 import numpy as np
+import xarray as xr
 
 from imod.logging import init_log_decorator
 from imod.mf6.boundary_condition import BoundaryCondition
@@ -11,6 +12,10 @@
 from imod.mf6.utilities.imod5_converter import convert_unit_rch_rate
 from imod.mf6.utilities.regrid import RegridderWeightsCache, _regrid_package_data
 from imod.mf6.validation import BOUNDARY_DIMS_SCHEMA, CONC_DIMS_SCHEMA
+from imod.msw.utilities.imod5_converter import (
+    get_cell_area_from_imod5_data,
+    is_msw_active_cell,
+)
 from imod.prepare.topsystem.allocation import ALLOCATION_OPTION, allocate_rch_cells
 from imod.schemata import (
     AllInsideNoDataSchema,
@@ -23,7 +28,7 @@
     IndexesSchema,
     OtherCoordsSchema,
 )
-from imod.typing import GridDataArray
+from imod.typing import GridDataArray, GridDataDict, Imod5DataDict
 from imod.typing.grid import (
     enforce_dim_order,
     is_planar_grid,
@@ -218,7 +223,7 @@ def from_imod5_data(
             # create an array indicating in which cells rch is active
             is_rch_cell = allocate_rch_cells(
                 ALLOCATION_OPTION.at_first_active,
-                new_idomain == 1,
+                new_idomain > 0,
                 planar_rate_regridded,
             )
 
@@ -227,7 +232,34 @@ def from_imod5_data(
             rch_rate = enforce_dim_order(rch_rate)
             new_package_data["rate"] = rch_rate
 
-        return Recharge(**new_package_data, validate=True, fixed_cell=False)
+        return cls(**new_package_data, validate=True, fixed_cell=False)
+
+    @classmethod
+    def from_imod5_cap_data(
+        cls,
+        imod5_data: Imod5DataDict,
+        target_dis: StructuredDiscretization,
+    ) -> "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 = cast(GridDataDict, imod5_data["cap"])
+
+        msw_area = get_cell_area_from_imod5_data(cap_data)
+        msw_active = is_msw_active_cell(target_dis, cap_data, msw_area)
+        active = msw_active.all
+
+        data = {}
+        layer_da = xr.full_like(
+            target_dis.dataset.coords["layer"], np.nan, dtype=np.float64
+        )
+        layer_da.loc[{"layer": 1}] = 0.0
+        data["rate"] = layer_da.where(active)
+
+        return cls(**data, validate=True, fixed_cell=False)
 
     @classmethod
     def get_regrid_methods(cls) -> RechargeRegridMethod:

imod/mf6/simulation.py

@@ -1361,11 +1361,13 @@ def from_imod5_data(
         times:  list[datetime]
             time discretization of the model to be imported. This is used for
             the following:
+
                 * Times of wells with associated timeseries are resampled to these times
                 * Start and end times in the list are used to repeat the stresses
                   of periodic data (e.g. river stages in iMOD5 for "summer", "winter")
                 * The simulation is discretized to these times, using
                   :meth:`imod.mf6.Modflow6Simulation.create_time_discretization`
+
         allocation_options: SimulationAllocationOptions, optional
             object containing the allocation options per package type. If you
             want a specific package to have a different allocation option, then

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,71 @@ 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."
+    )
+
+
+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,48 @@ 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/__init__.py

@@ -1,3 +1,4 @@
+from imod.msw.copy_files import FileCopier
 from imod.msw.coupler_mapping import CouplerMapping
 from imod.msw.grid_data import GridData
 from imod.msw.idf_mapping import IdfMapping
@@ -9,7 +10,7 @@
     InitialConditionsSavedState,
 )
 from imod.msw.landuse import LanduseOptions
-from imod.msw.meteo_grid import MeteoGrid
+from imod.msw.meteo_grid import MeteoGrid, MeteoGridCopy
 from imod.msw.meteo_mapping import EvapotranspirationMapping, PrecipitationMapping
 from imod.msw.model import MetaSwapModel
 from imod.msw.output_control import TimeOutputControl, VariableOutputControl