Issue #1290 infiltration from imod5

docs/api/changelog.rst

@@ -7,6 +7,16 @@ The format is based on `Keep a Changelog`_, and this project adheres to
 `Semantic Versioning`_.
 
 
+[Unreleased]
+------------
+
+Changed
+~~~~~~~
+
+- :class:`imod.msw.Infiltration`'s variables ``upward_resistance`` and
+  ``downward_resistance`` now require a ``subunit`` coordinate.
+
+
 [0.18.0]
 --------
 

imod/msw/infiltration.py

@@ -1,9 +1,35 @@
+from textwrap import dedent
+
 import xarray as xr
 
+from imod.logging import LogLevel, logger
 from imod.mf6.interfaces.iregridpackage import IRegridPackage
 from imod.msw.fixed_format import VariableMetaData
 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.grid import ones_like
+
+
+def deactivate_small_resistances_in_data(data: GridDataDict):
+    """
+    Deactivate cells where resistance smaller than 5 days are set to
+    -9999.0.
+    """
+    message = dedent("""Detected cells with resistances smaller than 5.0 in {var}, set
+    to inactive""")
+
+    for var in ["downward_resistance", "upward_resistance"]:
+        to_deactivate = data[var] < 5.0
+        if to_deactivate.any():
+            logger.log(
+                loglevel=LogLevel.WARNING,
+                message=message.format(var=var),
+                additional_depth=1,
+            )
+            data[var] = data[var].where(~to_deactivate, -9999.0)
+    return data
 
 
 class Infiltration(MetaSwapPackage, IRegridPackage):
@@ -19,11 +45,11 @@ class Infiltration(MetaSwapPackage, IRegridPackage):
         a subunit coordinate to describe different land uses.
     downward_resistance: array of floats (xr.DataArray)
         Describes the downward resisitance of SVAT units. Set to -9999.0 to make
-        MetaSWAP ignore this resistance. This array must not have a subunit
+        MetaSWAP ignore this resistance. This array must have a subunit
         coordinate.
     upward_resistance: array of floats (xr.DataArray)
         Describes the upward resistance of SVAT units. Set to -9999.0 to make
-        MetaSWAP ignore this resistance. This array must not have a subunit
+        MetaSWAP ignore this resistance. This array must have a subunit
         coordinate.
     bottom_resistance: array of floats (xr.DataArray)
         Describes the infiltration capacity of SVAT units. Set to -9999.0 to
@@ -32,9 +58,6 @@ class Infiltration(MetaSwapPackage, IRegridPackage):
     extra_storage_coefficient: array of floats (xr.DataArray)
         Extra storage coefficient of phreatic layer. This array must not have a
         subunit coordinate.
-    active: array of bools (xr.DataArray)
-        Describes whether SVAT units are active or not. This array must not have
-        a subunit coordinate.
     """
 
     _file_name = "infi_svat.inp"
@@ -47,10 +70,12 @@ class Infiltration(MetaSwapPackage, IRegridPackage):
         "extra_storage_coefficient": VariableMetaData(8, 0.01, 1.0, float),
     }
 
-    _with_subunit = ("infiltration_capacity",)
-    _without_subunit = (
+    _with_subunit = (
+        "infiltration_capacity",
         "downward_resistance",
         "upward_resistance",
+    )
+    _without_subunit = (
         "bottom_resistance",
         "extra_storage_coefficient",
     )
@@ -74,3 +99,28 @@ def __init__(
         self.dataset["extra_storage_coefficient"] = extra_storage_coefficient
 
         self._pkgcheck()
+
+    @classmethod
+    def from_imod5_data(cls, imod5_data: dict[str, GridDataDict]) -> "Infiltration":
+        cap_data = imod5_data["cap"]
+        data = {}
+        # Use runon resistance as downward resistance, and runoff for downward
+        # resistance
+        key_mapping = {
+            "infiltration_capacity": "infiltration_capacity",
+            "downward_resistance": "runon_resistance",
+            "upward_resistance": "runoff_resistance",
+        }
+        for var_rename, var_key in key_mapping.items():
+            data_ls = [
+                cap_data[f"{landuse}_{var_key}"] for landuse in ["rural", "urban"]
+            ]
+            data[var_rename] = concat_imod5(*data_ls)
+
+        data = deactivate_small_resistances_in_data(data)
+
+        like = data["downward_resistance"].isel(subunit=0, drop=True)
+        data["bottom_resistance"] = ones_like(like)
+        data["extra_storage_coefficient"] = ones_like(like)
+
+        return cls(**data)

imod/tests/fixtures/msw_model_fixture.py

@@ -226,8 +226,8 @@ def make_msw_model():
     # %% Infiltration
     msw_model["infiltration"] = msw.Infiltration(
         infiltration_capacity=xr.full_like(area, 1.0),
-        downward_resistance=xr.full_like(msw_grid, -9999.0),
-        upward_resistance=xr.full_like(msw_grid, -9999.0),
+        downward_resistance=xr.full_like(area, -9999.0),
+        upward_resistance=xr.full_like(area, -9999.0),
         bottom_resistance=xr.full_like(msw_grid, -9999.0),
         extra_storage_coefficient=xr.full_like(msw_grid, 0.1),
     )

imod/tests/test_msw/test_infiltration.py

@@ -1,78 +1,114 @@
 import tempfile
+from copy import deepcopy
 from pathlib import Path
 
 import numpy as np
+import pytest
 import xarray as xr
 from hypothesis import given, settings
 from hypothesis.strategies import floats
 from numpy import nan
 from numpy.testing import assert_almost_equal, assert_equal
+from pytest_cases import case, parametrize_with_cases
 
 from imod.mf6.utilities.regrid import (
     RegridderWeightsCache,
 )
 from imod.msw import Infiltration
 from imod.msw.fixed_format import format_fixed_width
+from imod.typing import GridDataDict
 
 
-def setup_infiltration_package(subunit, y, x, dy, dx):
-    infiltration_capacity = xr.DataArray(
+@pytest.fixture(scope="function")
+def coords_planar() -> dict:
+    x = [1.0, 2.0, 3.0]
+    y = [3.0, 2.0, 1.0]
+    dx = 1.0
+    dy = 1.0
+    return {"y": y, "x": x, "dx": dx, "dy": dy}
+
+
+@pytest.fixture(scope="function")
+def coords_subunit(coords_planar: dict) -> dict:
+    coords_subunit = deepcopy(coords_planar)
+    coords_subunit["subunit"] = [0, 1]
+    return coords_subunit
+
+
+@pytest.fixture(scope="function")
+def svat_index(coords_subunit: dict) -> tuple[xr.DataArray, np.ndarray]:
+    svat = xr.DataArray(
+        np.array(
+            [
+                [[0, 1, 0], [0, 0, 0], [0, 2, 0]],
+                [[0, 3, 0], [0, 4, 0], [0, 0, 0]],
+            ]
+        ),
+        dims=("subunit", "y", "x"),
+        coords=coords_subunit,
+    )
+    index = (svat != 0).values.ravel()
+    return svat, index
+
+
+@pytest.fixture(scope="function")
+def setup_infiltration_data(coords_planar, coords_subunit) -> GridDataDict:
+    data = {}
+    data["infiltration_capacity"] = xr.DataArray(
         np.array(
             [
                 [[0.5, 0.5, 0.5], [nan, nan, nan], [1.0, 1.0, 1.0]],
                 [[0.5, 0.5, 0.5], [1.0, 1.0, 1.0], [nan, nan, nan]],
             ]
         ),
         dims=("subunit", "y", "x"),
-        coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy},
+        coords=coords_subunit,
     )
-
-    downward_resistance = xr.DataArray(
-        np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]),
-        dims=("y", "x"),
-        coords={"y": y, "x": x, "dx": dx, "dy": dy},
+    data["downward_resistance"] = xr.DataArray(
+        np.array(
+            [
+                [[1.0, 2.0, 3.0], [nan, nan, nan], [7.0, 8.0, 9.0]],
+                [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [nan, nan, nan]],
+            ]
+        ),
+        dims=("subunit", "y", "x"),
+        coords=coords_subunit,
     )
-
-    upward_resistance = xr.DataArray(
-        np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]),
-        dims=("y", "x"),
-        coords={"y": y, "x": x, "dx": dx, "dy": dy},
+    data["upward_resistance"] = xr.DataArray(
+        np.array(
+            [
+                [[1.0, 2.0, 3.0], [nan, nan, nan], [7.0, 8.0, 9.0]],
+                [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [nan, nan, nan]],
+            ]
+        ),
+        dims=("subunit", "y", "x"),
+        coords=coords_subunit,
     )
-
-    bottom_resistance = xr.DataArray(
+    data["bottom_resistance"] = xr.DataArray(
         np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]),
         dims=("y", "x"),
-        coords={"y": y, "x": x, "dx": dx, "dy": dy},
+        coords=coords_planar,
     )
-
-    extra_storage_coefficient = xr.DataArray(
+    data["extra_storage_coefficient"] = xr.DataArray(
         np.array([[0.1, 0.2, 0.3], [0.4, 0.5, 0.6], [0.7, 0.8, 0.9]]),
         dims=("y", "x"),
-        coords={"y": y, "x": x, "dx": dx, "dy": dy},
+        coords=coords_planar,
     )
 
-    svat = xr.DataArray(
-        np.array(
-            [
-                [[0, 1, 0], [0, 0, 0], [0, 2, 0]],
-                [[0, 3, 0], [0, 4, 0], [0, 0, 0]],
-            ]
-        ),
-        dims=("subunit", "y", "x"),
-        coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy},
-    )
-    # fmt: on
-    index = (svat != 0).values.ravel()
+    return data
 
-    infiltration = Infiltration(
-        infiltration_capacity,
-        downward_resistance,
-        upward_resistance,
-        bottom_resistance,
-        extra_storage_coefficient,
-    )
 
-    return infiltration, svat, index
+@case(tags="r_low")
+def case_low_resistance(setup_infiltration_data: GridDataDict) -> GridDataDict:
+    return setup_infiltration_data
+
+
+@case(tags="r_high")
+def case_high_resistance(setup_infiltration_data: GridDataDict) -> GridDataDict:
+    data = setup_infiltration_data
+    data["downward_resistance"] += 10.0
+    data["upward_resistance"] += 10.0
+    return data
 
 
 @given(
@@ -108,8 +144,8 @@ def test_write(
 ):
     infiltration = Infiltration(
         xr.DataArray(infiltration_capacity).expand_dims(subunit=[0]),
-        xr.DataArray(downward_resistance),
-        xr.DataArray(upward_resistance),
+        xr.DataArray(downward_resistance).expand_dims(subunit=[0]),
+        xr.DataArray(upward_resistance).expand_dims(subunit=[0]),
         xr.DataArray(bottom_resistance),
         xr.DataArray(extra_storage_coefficient),
     )
@@ -175,14 +211,10 @@ def test_write(
     )
 
 
-def test_simple_model(fixed_format_parser):
-    x = [1.0, 2.0, 3.0]
-    y = [1.0, 2.0, 3.0]
-    subunit = [0, 1]
-    dx = 1.0
-    dy = 1.0
-    # fmt: off
-    infiltration, svat, index = setup_infiltration_package(subunit, y, x, dy, dx)
+@parametrize_with_cases("infiltration_data", cases=".", has_tag="r_low")
+def test_simple_model(fixed_format_parser, svat_index, infiltration_data):
+    svat, index = svat_index
+    infiltration = Infiltration(**infiltration_data)
 
     with tempfile.TemporaryDirectory() as output_dir:
         output_dir = Path(output_dir)
@@ -204,28 +236,58 @@ def test_simple_model(fixed_format_parser):
     )
 
 
-def test_regrid():
-    x = [1.0, 2.0, 3.0]
-    y = [3.0, 2.0, 1.0]
-    subunit = [0, 1]
-    dx = 1.0
-    dy = 1.0
-
-    infiltration, _, _ = setup_infiltration_package(subunit, y, x, dy, dx)
+@parametrize_with_cases("infiltration_data", cases=".", has_tag="r_low")
+def test_regrid(infiltration_data):
+    infiltration = Infiltration(**infiltration_data)
 
     x = [1.0, 1.5, 2.0, 2.5, 3.0]
     y = [3.0, 2.5, 2.0, 1.5, 1.0]
     subunit = [0, 1]
     dx = 0.5
     dy = 0.5
-    # fmt: off
     new_grid = xr.DataArray(
         dims=("subunit", "y", "x"),
-        coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy}
+        coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy},
     )
-    new_grid.values[:,:,:] = 1
+    new_grid.values[:, :, :] = 1
 
     regrid_context = RegridderWeightsCache()
-    regridded = infiltration.regrid_like(new_grid, regrid_context )
+    regridded = infiltration.regrid_like(new_grid, regrid_context)
     assert_almost_equal(regridded.dataset.coords["x"].values, x)
     assert_almost_equal(regridded.dataset.coords["y"].values, y)
+
+
+@parametrize_with_cases("data_infiltration", cases=".")
+def test_from_imod5_data(data_infiltration):
+    expected_pkg = Infiltration(**data_infiltration)
+    # Deactivate cells which have a resistance lower than 5.0
+    for var in ["upward_resistance", "downward_resistance"]:
+        da = expected_pkg.dataset[var]
+        to_deactivate = da < 5.0
+        expected_pkg.dataset[var] = da.where(~to_deactivate, -9999.0)
+
+    cap_data = {}
+    mapping_ls = [
+        ("rural_infiltration_capacity", "infiltration_capacity", 0),
+        ("urban_infiltration_capacity", "infiltration_capacity", 1),
+        ("rural_runoff_resistance", "upward_resistance", 0),
+        ("urban_runoff_resistance", "upward_resistance", 1),
+        ("rural_runon_resistance", "downward_resistance", 0),
+        ("urban_runon_resistance", "downward_resistance", 1),
+    ]
+    for cap_key, pkg_key, subunit_nr in mapping_ls:
+        cap_data[cap_key] = data_infiltration[pkg_key].sel(
+            subunit=subunit_nr, drop=True
+        )
+
+    imod5_data = {"cap": cap_data}
+    actual_pkg = Infiltration.from_imod5_data(imod5_data)
+
+    ones_vars = ["bottom_resistance", "extra_storage_coefficient"]
+    expected = expected_pkg.dataset.drop_vars(ones_vars)
+    actual = actual_pkg.dataset.drop_vars(ones_vars)
+
+    xr.testing.assert_equal(actual, expected)
+
+    for var in ones_vars:
+        assert (actual_pkg.dataset[var] == 1.0).all()

imod/tests/test_msw/test_ponding.py

@@ -116,7 +116,6 @@ def test_from_imod5_data():
     # Create cap data
     cap_data = {}
     mapping_ls = [
-        ("rural_runoff_resistance", "runoff_resistance", 0),
         ("rural_runoff_resistance", "runoff_resistance", 0),
         ("urban_runoff_resistance", "runoff_resistance", 1),
         ("rural_runon_resistance", "runon_resistance", 0),