Issue #1316 msw model from imod5 data

docs/api/changelog.rst

@@ -21,6 +21,12 @@ Added
   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.
 
 Fixed
 ~~~~~

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,

imod/mf6/rch.py

@@ -249,7 +249,8 @@ def from_imod5_cap_data(
         cap_data = cast(GridDataDict, imod5_data["cap"])
 
         msw_area = get_cell_area_from_imod5_data(cap_data)
-        active, _ = is_msw_active_cell(target_dis, cap_data, msw_area)
+        msw_active = is_msw_active_cell(target_dis, cap_data, msw_area)
+        active = msw_active.all
 
         data = {}
         layer_da = xr.full_like(

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/wel.py

@@ -1312,9 +1312,11 @@ def from_imod5_cap_data(cls, imod5_data: Imod5DataDict):
             (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.

imod/msw/grid_data.py

@@ -1,14 +1,8 @@
-from typing import Optional
-
 import numpy as np
 import xarray as xr
 
 from imod.mf6 import StructuredDiscretization
 from imod.mf6.interfaces.iregridpackage import IRegridPackage
-from imod.mf6.regrid.regrid_schemes import (
-    RegridMethodType,
-)
-from imod.mf6.utilities.regrid import RegridderWeightsCache
 from imod.msw.fixed_format import VariableMetaData
 from imod.msw.pkgbase import MetaSwapPackage
 from imod.msw.regrid.regrid_schemes import GridDataRegridMethod
@@ -18,7 +12,8 @@
     get_rootzone_depth_from_imod5_data,
     is_msw_active_cell,
 )
-from imod.typing import GridDataDict
+from imod.msw.utilities.mask import MetaSwapActive, mask_and_broadcast_pkg_data
+from imod.typing import Imod5DataDict
 from imod.util.spatial import get_cell_area, spatial_reference
 
 
@@ -129,30 +124,48 @@ def _pkgcheck(self):
     @classmethod
     def from_imod5_data(
         cls,
-        imod5_data: dict[str, GridDataDict],
+        imod5_data: Imod5DataDict,
         target_dis: StructuredDiscretization,
-        regridder_types: Optional[RegridMethodType] = None,
-        regrid_cache: RegridderWeightsCache = RegridderWeightsCache(),
-    ) -> "GridData":
-        # Get iMOD5 capillary zone data
+    ) -> tuple["GridData", MetaSwapActive]:
+        """
+        Construct a MetaSWAP GridData package from iMOD5 data in the CAP
+        package, loaded with the :func:`imod.formats.prj.open_projectfile_data`
+        function.
+
+        Method does the following things:
+
+        - Computes rural area from the wetted area and urban area grids.
+        - Sets a second subunit for urban landuse (== 18).
+        - Rootzone depth is converted from cm's to m's.
+        - Determines where MetaSWAP should be active based on area grids,
+          boundary grid, and MODFLOW6 idomain.
+
+        Parameters
+        ----------
+        imod5_data: Imod5DataDict
+            iMOD5 data as returned by
+            :func:`imod.formats.prj.open_projectfile_data`
+        target_dis: imod.mf6.StructuredDiscretization
+            MODFLOW6 discretization to couple the MetaSWAP model to.
+
+        Returns
+        -------
+        imod.msw.GridData
+            GridData package
+        MetaSwapActive
+            Dataclass containing where MetaSwAP is active, per subunit, as well
+            as aggregated over subunits.
+        """
         imod5_cap = imod5_data["cap"]
 
         data = {}
         data["area"] = get_cell_area_from_imod5_data(imod5_cap)
         data["landuse"] = get_landuse_from_imod5_data(imod5_cap)
         data["rootzone_depth"] = get_rootzone_depth_from_imod5_data(imod5_cap)
         data["surface_elevation"] = imod5_cap["surface_elevation"]
-        data["soil_physical_unit"] = imod5_cap["soil_physical_unit"]
-
-        active, subunit_active = is_msw_active_cell(target_dis, imod5_cap, data["area"])
+        data["soil_physical_unit"] = imod5_cap["soil_physical_unit"].astype(int)
 
-        data_active = {
-            key: (
-                griddata.where(subunit_active)
-                if key in cls._with_subunit
-                else griddata.where(active)
-            )
-            for key, griddata in data.items()
-        }
-        data_active["active"] = active
-        return cls(**data_active)
+        msw_active = is_msw_active_cell(target_dis, imod5_cap, data["area"])
+        data_active = mask_and_broadcast_pkg_data(cls, data, msw_active)
+        data_active["active"] = msw_active.all
+        return cls(**data_active), msw_active

imod/msw/infiltration.py

@@ -8,7 +8,8 @@
 from imod.msw.pkgbase import MetaSwapPackage
 from imod.msw.regrid.regrid_schemes import InfiltrationRegridMethod
 from imod.msw.utilities.common import concat_imod5
-from imod.typing import GridDataDict
+from imod.msw.utilities.mask import MaskValues
+from imod.typing import GridDataDict, Imod5DataDict
 from imod.typing.grid import ones_like
 
 
@@ -28,7 +29,7 @@ def deactivate_small_resistances_in_data(data: GridDataDict):
                 message=message.format(var=var),
                 additional_depth=1,
             )
-            data[var] = data[var].where(~to_deactivate, -9999.0)
+            data[var] = data[var].where(~to_deactivate, MaskValues.default)
     return data
 
 
@@ -101,7 +102,7 @@ def __init__(
         self._pkgcheck()
 
     @classmethod
-    def from_imod5_data(cls, imod5_data: dict[str, GridDataDict]) -> "Infiltration":
+    def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "Infiltration":
         cap_data = imod5_data["cap"]
         data = {}
         # Use runon resistance as downward resistance, and runoff for downward
@@ -120,7 +121,7 @@ def from_imod5_data(cls, imod5_data: dict[str, GridDataDict]) -> "Infiltration":
         data = deactivate_small_resistances_in_data(data)
 
         like = data["downward_resistance"].isel(subunit=0, drop=True)
-        data["bottom_resistance"] = ones_like(like)
+        data["bottom_resistance"] = ones_like(like) * MaskValues.default
         data["extra_storage_coefficient"] = ones_like(like)
 
         return cls(**data)

imod/msw/meteo_grid.py

@@ -1,7 +1,7 @@
 import csv
 from pathlib import Path
 from shutil import copyfile
-from typing import Optional, Union, cast
+from typing import Optional, Union
 
 import numpy as np
 import pandas as pd
@@ -13,6 +13,7 @@
 from imod.msw.regrid.regrid_schemes import MeteoGridRegridMethod
 from imod.msw.timeutil import to_metaswap_timeformat
 from imod.msw.utilities.common import find_in_file_list
+from imod.msw.utilities.mask import MaskValues
 from imod.typing import Imod5DataDict
 from imod.util.regrid_method_type import EmptyRegridMethod, RegridMethodType
 
@@ -176,7 +177,9 @@ def write(self, directory: Union[str, Path], *args):
                 path = (directory / self._meteo_dirname / str(varname)).with_suffix(
                     ".asc"
                 )
-                imod.rasterio.save(path, self.dataset[str(varname)], nodata=-9999.0)
+                imod.rasterio.save(
+                    path, self.dataset[str(varname)], nodata=MaskValues.default
+                )
 
     def _pkgcheck(self):
         for varname in self.dataset.data_vars:
@@ -223,7 +226,7 @@ def write(self, directory: Path | str, *args):
 
     @classmethod
     def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "MeteoGridCopy":
-        paths = cast(list[str], imod5_data["extra"]["paths"])
+        paths = imod5_data["extra"]["paths"]
         filepath = find_in_file_list(cls._file_name, paths)
 
         return cls(filepath)

imod/msw/meteo_mapping.py

@@ -1,7 +1,7 @@
 from copy import deepcopy
 from pathlib import Path
 from textwrap import dedent
-from typing import Any, Optional, TextIO, cast
+from typing import Any, Optional, TextIO
 
 import numpy as np
 import pandas as pd
@@ -58,7 +58,7 @@ def open_first_meteo_grid(mete_grid_path: str | Path, column_nr: int) -> xr.Data
 
 
 def open_first_meteo_grid_from_imod5_data(imod5_data: Imod5DataDict, column_nr: int):
-    paths = cast(list[str], imod5_data["extra"]["paths"])
+    paths = imod5_data["extra"]["paths"]
     metegrid_path = find_in_file_list("mete_grid.inp", paths)
     return open_first_meteo_grid(metegrid_path, column_nr=column_nr)