diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst index 59642ac66..68ebc2591 100644 --- a/docs/api/changelog.rst +++ b/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] -------- diff --git a/imod/msw/infiltration.py b/imod/msw/infiltration.py index bf300761f..830365681 100644 --- a/imod/msw/infiltration.py +++ b/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) diff --git a/imod/tests/fixtures/msw_model_fixture.py b/imod/tests/fixtures/msw_model_fixture.py index 3b7c73c80..1fcaf9015 100644 --- a/imod/tests/fixtures/msw_model_fixture.py +++ b/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), ) diff --git a/imod/tests/test_msw/test_infiltration.py b/imod/tests/test_msw/test_infiltration.py index d3cb0a003..b87940a88 100644 --- a/imod/tests/test_msw/test_infiltration.py +++ b/imod/tests/test_msw/test_infiltration.py @@ -1,22 +1,60 @@ 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]], @@ -24,55 +62,53 @@ def setup_infiltration_package(subunit, y, x, dy, dx): ] ), 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() diff --git a/imod/tests/test_msw/test_ponding.py b/imod/tests/test_msw/test_ponding.py index e7bb47e06..94f0d64d0 100644 --- a/imod/tests/test_msw/test_ponding.py +++ b/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),