Issue #1285 grid data from imod5

imod/msw/grid_data.py

@@ -1,13 +1,80 @@
+from typing import Optional
+
 import numpy as np
 import xarray as xr
 
+from imod.logging import LogLevel, logger
+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
+from imod.typing import GridDataArray, GridDataDict
+from imod.typing.grid import concat, ones_like
 from imod.util.spatial import get_cell_area, spatial_reference
 
 
+def _concat_subunits(arg1: GridDataArray, arg2: GridDataArray):
+    return concat([arg1, arg2], dim="subunit").assign_coords(subunit=[0, 1])
+
+
+def get_cell_area_from_imod5_data(
+    imod5_cap: GridDataDict,
+) -> GridDataArray:
+    # area's per type of svats
+    mf6_area = get_cell_area(imod5_cap["wetted_area"])
+    wetted_area = imod5_cap["wetted_area"]
+    urban_area = imod5_cap["urban_area"]
+    rural_area = mf6_area - (wetted_area + urban_area)
+    if (wetted_area > mf6_area).any():
+        logger.log(
+            loglevel=LogLevel.WARNING,
+            message=f"wetted area was set to the max cell area of {mf6_area}",
+            additional_depth=0,
+        )
+        wetted_area = wetted_area.where(wetted_area <= mf6_area, other=mf6_area)
+    if (rural_area < 0.0).any():
+        logger.log(
+            loglevel=LogLevel.WARNING,
+            message="found urban area > than (cel-area - wetted area). Urban area was set to 0",
+            additional_depth=0,
+        )
+        urban_area = urban_area.where(rural_area > 0.0, other=0.0)
+    rural_area = mf6_area - (wetted_area + urban_area)
+    return _concat_subunits(rural_area, urban_area)
+
+
+def get_landuse_from_imod5_data(
+    imod5_cap: GridDataDict,
+) -> GridDataArray:
+    """
+    Get landuse from imod5 capillary zone data. This adds two subunits, one
+    based on the landuse grid, which specifies rural landuse. The other
+    specifies urban landuse, which is coded to value 18.
+    """
+    rural_landuse = imod5_cap["landuse"]
+    # Urban landuse = 18
+    urban_landuse = ones_like(rural_landuse) * 18
+    return _concat_subunits(rural_landuse, urban_landuse)
+
+
+def get_rootzone_depth_from_imod5_data(
+    imod5_cap: GridDataDict,
+) -> GridDataArray:
+    """
+    Get rootzone depth from imod5 capillary zone data. Also does a unit
+    conversion: iMOD5 specifies rootzone thickness in centimeters, whereas
+    MetaSWAP requires rootzone depth in meters.
+    """
+    rootzone_thickness = imod5_cap["rootzone_thickness"] * 0.01
+    # rootzone depth is equal for both svats.
+    return _concat_subunits(rootzone_thickness, rootzone_thickness)
+
+
 class GridData(MetaSwapPackage, IRegridPackage):
     """
     This contains the grid data of MetaSWAP.
@@ -111,3 +178,37 @@ def _pkgcheck(self):
             raise ValueError(
                 "Provided area grid with total areas larger than cell area"
             )
+
+    @classmethod
+    def from_imod5_data(
+        cls,
+        imod5_data: dict[str, GridDataDict],
+        target_dis: StructuredDiscretization,
+        regridder_types: Optional[RegridMethodType] = None,
+        regrid_cache: RegridderWeightsCache = RegridderWeightsCache(),
+    ) -> "GridData":
+        # Get iMOD5 capillary zone data
+        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"]
+
+        mf6_top_active = target_dis["idomain"].isel(layer=0, drop=True)
+        subunit_active = (
+            (imod5_cap["boundary"] == 1) & (data["area"] > 0) & (mf6_top_active >= 1)
+        )
+        active = subunit_active.all(dim="subunit")
+        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)

imod/msw/pkgbase.py

@@ -240,3 +240,6 @@ def regrid_like(
 
     def get_regrid_methods(self) -> RegridMethodType:
         return deepcopy(self._regrid_method)
+
+    def from_imod5_data(self, *args, **kwargs):
+        raise NotImplementedError("Method not implemented for this package.")