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/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/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,47 @@ 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,7 @@
 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.typing import GridDataDict, Imod5DataDict
 from imod.typing.grid import ones_like
 
 
@@ -101,7 +101,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 +120,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) * -9999.0
         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
@@ -223,7 +223,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)
 

imod/msw/model.py

@@ -1,14 +1,17 @@
 import collections
 from copy import copy
+from datetime import datetime
 from pathlib import Path
-from typing import Optional, Tuple, Union
+from typing import Any, Optional, Tuple, Union, cast
 
 import jinja2
 import numpy as np
+import xarray as xr
 
 from imod.mf6.dis import StructuredDiscretization
 from imod.mf6.mf6_wel_adapter import Mf6Wel
 from imod.mf6.utilities.regrid import RegridderWeightsCache
+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
@@ -20,11 +23,20 @@
     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.output_control import TimeOutputControl
+from imod.msw.ponding import Ponding
+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.mask import mask_and_broadcast_cap_data
+from imod.msw.utilities.parse import read_para_sim
+from imod.msw.utilities.select import drop_layer_dim_cap_data
 from imod.msw.vegetation import AnnualCropFactors
+from imod.typing import Imod5DataDict
 from imod.util.regrid_method_type import RegridderType
 
 REQUIRED_PACKAGES = (
@@ -88,6 +100,7 @@ class MetaSwapModel(Model):
     ----------
     unsaturated_database: Path-like or str
         Path to the MetaSWAP soil physical database folder.
+    settings: dict
     """
 
     _pkg_id = "model"
@@ -99,10 +112,18 @@ class MetaSwapModel(Model):
         "{%endfor%}"
     )
 
-    def __init__(self, unsaturated_database):
+    def __init__(
+        self,
+        unsaturated_database: Path | str,
+        settings: Optional[dict[str, Any]] = None,
+    ):
         super().__init__()
 
-        self.simulation_settings = copy(DEFAULT_SETTINGS)
+        if settings is None:
+            self.simulation_settings = copy(DEFAULT_SETTINGS)
+        else:
+            self.simulation_settings = settings
+
         self.simulation_settings["unsa_svat_path"] = (
             self._render_unsaturated_database_path(unsaturated_database)
         )
@@ -210,6 +231,7 @@ def write(
         directory: Union[str, Path],
         mf6_dis: StructuredDiscretization,
         mf6_wel: Mf6Wel,
+        validate: bool = True,
     ):
         """
         Write packages and simulation settings (para_sim.inp).
@@ -221,9 +243,10 @@ def write(
         """
 
         # Model checks
-        self._check_required_packages()
-        self._check_vegetation_indices_in_annual_crop_factors()
-        self._check_landuse_indices_in_lookup_options()
+        if validate:
+            self._check_required_packages()
+            self._check_vegetation_indices_in_annual_crop_factors()
+            self._check_landuse_indices_in_lookup_options()
 
         # Force to Path
         directory = Path(directory)
@@ -258,7 +281,7 @@ def regrid_like(
         regrid_context: Optional[RegridderWeightsCache] = None,
         regridder_types: Optional[dict[str, Tuple[RegridderType, str]]] = None,
     ) -> "MetaSwapModel":
-        unsat_database = self.simulation_settings["unsa_svat_path"]
+        unsat_database = cast(str, self.simulation_settings["unsa_svat_path"])
         regridded_model = MetaSwapModel(unsat_database)
 
         target_grid = mf6_regridded_dis["idomain"]
@@ -282,3 +305,66 @@ def regrid_like(
             regridded_model[mod2svat_name] = CouplerMapping()
 
         return regridded_model
+
+    @classmethod
+    def from_imod5_data(
+        cls,
+        imod5_data: Imod5DataDict,
+        target_dis: StructuredDiscretization,
+        times: list[datetime],
+    ):
+        """
+        Construct a MetaSWAP model from iMOD5 data in the CAP package, loaded
+        with the :func:`imod.formats.prj.open_projectfile_data` function.
+
+        Parameters
+        ----------
+        imod5_data: dict
+            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.
+        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.
+
+        Returns
+        -------
+        MetaSwapModel
+            MetaSWAP model imported from imod5 data.
+        """
+        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)
+        model = cls(unsa_svat_path, parasim_settings)
+        model["grid"], msw_active = GridData.from_imod5_data(
+            imod5_cap_no_layer, target_dis
+        )
+        cap_data_masked = mask_and_broadcast_cap_data(
+            imod5_cap_no_layer["cap"], msw_active
+        )
+        imod5_masked: Imod5DataDict = {
+            "cap": cap_data_masked,
+            "extra": {"paths": extra_paths},
+        }
+        model["infiltration"] = Infiltration.from_imod5_data(imod5_masked)
+        model["ponding"] = Ponding.from_imod5_data(imod5_masked)
+        model["sprinkling"] = Sprinkling.from_imod5_data(imod5_masked)
+        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"]):
+            model["scaling_factor"] = ScalingFactors.from_imod5_data(imod5_masked)
+        area = model["grid"]["area"].isel(subunit=0, drop=True)
+        model["idf_mapping"] = IdfMapping(area, -9999.0)
+        model["coupling"] = CouplerMapping()
+        model["extra_files"] = FileCopier.from_imod5_data(imod5_masked)
+
+        times_da = xr.DataArray(times, coords={"time": times}, dims=("time",))
+        model["time_oc"] = TimeOutputControl(times_da)
+
+        return model

imod/msw/ponding.py

@@ -8,7 +8,7 @@
 from imod.msw.pkgbase import DataDictType, MetaSwapPackage
 from imod.msw.regrid.regrid_schemes import PondingRegridMethod
 from imod.msw.utilities.common import concat_imod5
-from imod.typing import GridDataDict, IntArray
+from imod.typing import Imod5DataDict, IntArray
 
 
 class Ponding(MetaSwapPackage, IRegridPackage):
@@ -73,7 +73,7 @@ def _render(self, file: TextIO, index: IntArray, svat: xr.DataArray, *args: Any)
         return self.write_dataframe_fixed_width(file, dataframe)
 
     @classmethod
-    def from_imod5_data(cls, imod5_data: dict[str, GridDataDict]) -> "Ponding":
+    def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "Ponding":
         """
         Concatenate ponding depths along subunits
         """