Issue #365 add regrid to msw packages (#982)

Fixes #365

# Description
Adds regrid functionality to a number of metaswap packages
Adds unittests for regridding these packages

# Checklist

- [X] Links to correct issue
- [ ] Update changelog, if changes affect users
- [X] PR title starts with ``Issue #nr``, e.g. ``Issue #737``
- [X] Unit tests were added
- [ ] **If feature added**: Added/extended example

imod/msw/infiltration.py

@@ -1,10 +1,14 @@
 import xarray as xr
 
+from imod.mf6.interfaces.iregridpackage import IRegridPackage
+from imod.mf6.utilities.regrid import (
+    RegridderType,
+)
 from imod.msw.fixed_format import VariableMetaData
 from imod.msw.pkgbase import MetaSwapPackage
 
 
-class Infiltration(MetaSwapPackage):
+class Infiltration(MetaSwapPackage, IRegridPackage):
     """
     This contains the infiltration data.
 
@@ -54,6 +58,21 @@ class Infiltration(MetaSwapPackage):
     )
     _to_fill = ()
 
+    _regrid_method = {
+        "infiltration_capacity": (RegridderType.OVERLAP, "mean"),
+        "downward_resistance": (RegridderType.OVERLAP, "mean"),
+        "upward_resistance": (
+            RegridderType.OVERLAP,
+            "mean",
+        ),
+        "longitudinal_vertical": (
+            RegridderType.OVERLAP,
+            "mean",
+        ),
+        "bottom_resistance": (RegridderType.OVERLAP, "mean"),
+        "extra_storage_coefficient": (RegridderType.OVERLAP, "mean"),
+    }
+
     def __init__(
         self,
         infiltration_capacity: xr.DataArray,

imod/msw/initial_conditions.py

@@ -1,17 +1,20 @@
 import pathlib
 import shutil
 
+from imod.mf6.interfaces.iregridpackage import IRegridPackage
 from imod.msw.fixed_format import VariableMetaData
 from imod.msw.pkgbase import MetaSwapPackage
 
 
-class InitialConditionsEquilibrium(MetaSwapPackage):
+class InitialConditionsEquilibrium(MetaSwapPackage, IRegridPackage):
     """
     Use an equilibrium profile to initialize the model.
 
     This class is responsible for the file `init_svat.inp`
     """
 
+    _regrid_method = {}
+
     _file_name = "init_svat.inp"
     _option = "Equilibrium"
     _metadata_dict = {}
@@ -23,7 +26,7 @@ def _render(self, file, *args):
         file.write(self._option + "\n")
 
 
-class InitialConditionsRootzonePressureHead(MetaSwapPackage):
+class InitialConditionsRootzonePressureHead(MetaSwapPackage, IRegridPackage):
     """
     Use the pF-value of the root zone pressure head as initial condition.
 
@@ -40,6 +43,7 @@ class InitialConditionsRootzonePressureHead(MetaSwapPackage):
     _metadata_dict = {
         "initial_pF": VariableMetaData(6, 0.0, 6.0, float),
     }
+    _regrid_method = {}
 
     def __init__(self, initial_pF=2.2):
         super().__init__()
@@ -53,7 +57,7 @@ def _render(self, file, *args):
         self.write_dataframe_fixed_width(file, dataframe)
 
 
-class InitialConditionsPercolation(MetaSwapPackage):
+class InitialConditionsPercolation(MetaSwapPackage, IRegridPackage):
     """
     The precipitation intensity at the starting time (iybg, tdbg in
     PARA_SIM.INP) is used for initializing the percolation flux in the profiles.
@@ -68,6 +72,7 @@ class InitialConditionsPercolation(MetaSwapPackage):
     _file_name = "init_svat.inp"
     _option = "MeteoInputP"
     _metadata_dict = {}
+    _regrid_method = {}
 
     def __init__(self):
         super().__init__()

imod/msw/landuse.py

@@ -1,8 +1,9 @@
+from imod.mf6.interfaces.iregridpackage import IRegridPackage
 from imod.msw.fixed_format import VariableMetaData
 from imod.msw.pkgbase import MetaSwapPackage
 
 
-class LanduseOptions(MetaSwapPackage):
+class LanduseOptions(MetaSwapPackage, IRegridPackage):
     """
     Land use options. This object is responsible for luse_svat.inp
 
@@ -53,6 +54,12 @@ class LanduseOptions(MetaSwapPackage):
         Hoyningen.
     interception_intercept: array of floats (xr.DataArray)
         Intercept of the interception evaporation curve. Pun unintended.
+    interception_capacity_per_LAI_Rutter: this initializer argument is for
+        internal use only. Users can ignore it and should only provide
+        interception_capacity_per_LAI
+    interception_capacity_per_LAI_VonHoyningen: this initializer argument is for
+        internal use only. Users can ignore it and should only provide
+        interception_capacity_per_LAI
 
     Notes
     -----
@@ -116,6 +123,8 @@ class LanduseOptions(MetaSwapPackage):
 
     _file_name = "luse_svat.inp"
 
+    _regrid_method = {}
+
     def __init__(
         self,
         landuse_name,
@@ -137,8 +146,10 @@ def __init__(
         start_sprinkling_season,
         end_sprinkling_season,
         interception_option,
-        interception_capacity_per_LAI,
-        interception_intercept,
+        interception_capacity_per_LAI=None,
+        interception_capacity_per_LAI_Rutter=None,
+        interception_capacity_per_LAI_VonHoyningen=None,
+        interception_intercept=None,
     ):
         super().__init__()
         self.dataset["landuse_name"] = landuse_name
@@ -160,12 +171,22 @@ def __init__(
         self.dataset["start_sprinkling_season"] = start_sprinkling_season
         self.dataset["end_sprinkling_season"] = end_sprinkling_season
         self.dataset["interception_option"] = interception_option
-        self.dataset["interception_capacity_per_LAI_Rutter"] = (
-            interception_capacity_per_LAI
-        )
-        self.dataset["interception_capacity_per_LAI_VonHoyningen"] = (
-            interception_capacity_per_LAI
-        )
+        if interception_capacity_per_LAI is not None:
+            self.dataset["interception_capacity_per_LAI_Rutter"] = (
+                interception_capacity_per_LAI
+            )
+        else:
+            self.dataset["interception_capacity_per_LAI_Rutter"] = (
+                interception_capacity_per_LAI_Rutter
+            )
+        if interception_capacity_per_LAI is not None:
+            self.dataset["interception_capacity_per_LAI_VonHoyningen"] = (
+                interception_capacity_per_LAI
+            )
+        else:
+            self.dataset["interception_capacity_per_LAI_VonHoyningen"] = (
+                interception_capacity_per_LAI_VonHoyningen
+            )
         self.dataset["interception_intercept"] = interception_intercept
 
         self._pkgcheck()

imod/msw/pkgbase.py

@@ -1,12 +1,18 @@
 import abc
 from pathlib import Path
-from typing import Union
+from typing import Any, Optional, Tuple, Union
 
 import numpy as np
 import pandas as pd
 import xarray as xr
 
+from imod.mf6.utilities.regrid import (
+    RegridderType,
+    RegridderWeightsCache,
+    _regrid_like,
+)
 from imod.msw.fixed_format import format_fixed_width
+from imod.typing.grid import GridDataArray, GridDataset
 
 
 class MetaSwapPackage(abc.ABC):
@@ -28,6 +34,14 @@ def __getitem__(self, key):
     def __setitem__(self, key, value):
         self.dataset.__setitem__(key, value)
 
+    @property
+    def dataset(self) -> GridDataset:
+        return self._dataset
+
+    @dataset.setter
+    def dataset(self, value):
+        self._dataset = value
+
     def isel(self):
         raise NotImplementedError(
             f"Selection on packages not yet supported. "
@@ -131,3 +145,49 @@ def _pkgcheck(self):
                     f"Variable '{var}' in {self.__class__} should not "
                     "contain 'subunit' coordinate"
                 )
+
+    def _valid(self, value):
+        return True
+
+    def get_non_grid_data(self, grid_names: list[str]) -> dict[str, Any]:
+        """
+        This function copies the attributes of a dataset that are scalars, such as options.
+
+        parameters
+        ----------
+        grid_names: list of str
+            the names of the attribbutes of a dataset that are grids.
+        """
+        result = {}
+        all_non_grid_data = list(self.dataset.keys())
+        for name in (
+            gridname for gridname in grid_names if gridname in all_non_grid_data
+        ):
+            all_non_grid_data.remove(name)
+        for name in all_non_grid_data:
+            result[name] = self.dataset[name]
+        return result
+
+    @property
+    def auxiliary_data_fields(self) -> dict[str, str]:
+        return {}
+
+    def is_regridding_supported(self) -> bool:
+        return
+
+    def regrid_like(
+        self,
+        target_grid: GridDataArray,
+        regrid_context: RegridderWeightsCache,
+        regridder_types: Optional[dict[str, Tuple[RegridderType, str]]] = None,
+    ) -> "MetaSwapPackage":
+        try:
+            result = _regrid_like(self, target_grid, regrid_context, regridder_types)
+        except ValueError as e:
+            raise e
+        except Exception as e:
+            raise ValueError(f"package could not be regridded:{e}")
+        return result
+
+    def get_regrid_methods(self) -> Optional[dict[str, Tuple[RegridderType, str]]]:
+        return self._regrid_method

imod/msw/ponding.py

@@ -1,10 +1,14 @@
 import pandas as pd
 
+from imod.mf6.interfaces.iregridpackage import IRegridPackage
+from imod.mf6.utilities.regrid import (
+    RegridderType,
+)
 from imod.msw.fixed_format import VariableMetaData
 from imod.msw.pkgbase import MetaSwapPackage
 
 
-class Ponding(MetaSwapPackage):
+class Ponding(MetaSwapPackage, IRegridPackage):
     """
     Set ponding related parameters for MetaSWAP. This class is responsible for
     the svat2swnr_roff.inp file. Currently, we do not support ponds coupled to
@@ -39,6 +43,15 @@ class Ponding(MetaSwapPackage):
     _without_subunit = ()
     _to_fill = ()
 
+    _regrid_method = {
+        "ponding_depth": (RegridderType.OVERLAP, "mean"),
+        "runon_resistance": (RegridderType.OVERLAP, "mean"),
+        "runoff_resistance": (
+            RegridderType.OVERLAP,
+            "mean",
+        ),
+    }
+
     def __init__(self, ponding_depth, runon_resistance, runoff_resistance) -> None:
         super().__init__()
         self.dataset["ponding_depth"] = ponding_depth

imod/msw/scaling_factors.py

@@ -1,8 +1,12 @@
+from imod.mf6.interfaces.iregridpackage import IRegridPackage
+from imod.mf6.utilities.regrid import (
+    RegridderType,
+)
 from imod.msw.fixed_format import VariableMetaData
 from imod.msw.pkgbase import MetaSwapPackage
 
 
-class ScalingFactors(MetaSwapPackage):
+class ScalingFactors(MetaSwapPackage, IRegridPackage):
     """
     This package allows you to do three things:
         1. Set scaling factors for some inputs in the soil physical database,
@@ -58,6 +62,19 @@ class ScalingFactors(MetaSwapPackage):
     _without_subunit = ("depth_perched_water_table",)
     _to_fill = ()
 
+    _regrid_method = {
+        "scale_soil_moisture": (RegridderType.OVERLAP, "mean"),
+        "scale_hydraulic_conductivity": (RegridderType.OVERLAP, "mean"),
+        "scale_pressure_head": (
+            RegridderType.OVERLAP,
+            "mean",
+        ),
+        "depth_perched_water_table": (
+            RegridderType.OVERLAP,
+            "mean",
+        ),
+    }
+
     def __init__(
         self,
         scale_soil_moisture,

imod/msw/vegetation.py

@@ -1,11 +1,12 @@
 import numpy as np
 import xarray as xr
 
+from imod.mf6.interfaces.iregridpackage import IRegridPackage
 from imod.msw.fixed_format import VariableMetaData
 from imod.msw.pkgbase import MetaSwapPackage
 
 
-class AnnualCropFactors(MetaSwapPackage):
+class AnnualCropFactors(MetaSwapPackage, IRegridPackage):
     """
     For each vegetation type specify a yearly trend in vegetation factors and
     interception characteristics. These are used if WOFOST is not used.
@@ -52,6 +53,8 @@ class AnnualCropFactors(MetaSwapPackage):
         "ponding_factor": VariableMetaData(8, 0.01, 10.0, float),
     }
 
+    _regrid_method = {}
+
     def __init__(
         self,
         soil_cover: xr.DataArray,