Issue #1787 mask topsystem packages

docs/api/changelog.rst

@@ -18,20 +18,35 @@ Added
   :meth:`imod.msw.SprinklingPoints.from_imod5_data`.
 - :class:`imod.mf6.LayeredWell.from_imod5_cap_data` now also supports loading
   wells from IPF files in an iMOD5 CAP dataset.
+- Added :meth:`imod.mf6.GroundwaterFlowModel.mask_packages` and
+  :meth:`imod.mf6.GroundwaterTransportModel.mask_packages` to mask specific
+  packages of a groundwater flow model and a groundwater transport model
+  respectively.
 
 Fixed
 ~~~~~
 
 - Fixed resampling in :meth:`imod.mf6.Well.from_imod5_data` and
   :meth:`imod.mf6.LayeredWell.from_imod5_data` when simulation timesteps precede
   the first well timestep.
+- :meth:`imod.mf6.GroundwaterFlowModel.from_imod5_data` now masks cells in
+  topsystem packages (:class:`imod.mf6.River`,
+  :class:`imod.mf6.GeneralHeadBoundary`, :class:`imod.mf6.Drainage`,
+  :class:`imod.mf6.Recharge`) where IBOUND is less than 0.
+- :meth:`imod.msw.MetaSwapModel.from_imod5_data` now masks cells where IBOUND is
+  less than 0.
 
 Changed
 ~~~~~~~
 
 - Deprecated :class:`imod.msw.Sprinkling` in favor of
   :class:`imod.msw.SprinklingGrid`. Call :class:`imod.msw.SprinklingGrid` to get
   the same behavior as you were used to.
+- If ``states_for_boundary`` is provided to
+  :meth:`imod.mf6.GroundwaterFlowModel.clip_box`, topsystem packages
+  (:class:`imod.mf6.River`, :class:`imod.mf6.GeneralHeadBoundary`,
+  :class:`imod.mf6.Drainage`, :class:`imod.mf6.Recharge`) will also be masked
+  where constant head cells are placed. 
 
 [1.1.0] - 2026-08-03
 --------------------

docs/api/mf6.rst

@@ -46,6 +46,7 @@ Model objects & methods
     Modflow6Simulation.set_validation_settings
     GroundwaterFlowModel
     GroundwaterFlowModel.mask_all_packages
+    GroundwaterFlowModel.mask_packages
     GroundwaterFlowModel.prepare_wel_for_mf6
     GroundwaterFlowModel.regrid_like
     GroundwaterFlowModel.dump
@@ -60,6 +61,7 @@ Model objects & methods
     GroundwaterFlowModel.get_diskey
     GroundwaterTransportModel
     GroundwaterTransportModel.mask_all_packages
+    GroundwaterTransportModel.mask_packages
     GroundwaterTransportModel.dump
     GroundwaterTransportModel.clip_box
     GroundwaterTransportModel.regrid_like

imod/common/interfaces/iboundarycondition.py

@@ -0,0 +1,12 @@
+from abc import ABC
+from typing import ClassVar
+
+from imod.common.interfaces.ipackage import IPackage
+
+
+class IBoundaryCondition(IPackage, ABC):
+    """
+    Interface for boundary condition packages with stress-period variables.
+    """
+
+    _period_data: ClassVar[tuple[str, ...]]

imod/common/interfaces/imodel.py

@@ -2,6 +2,7 @@
 from typing import Any, Optional, Tuple
 
 from imod.common.interfaces.idict import IDict
+from imod.common.interfaces.ipackage import IPackage
 from imod.common.statusinfo import StatusInfoBase
 from imod.mf6.validation_settings import ValidationSettings
 from imod.typing import GridDataArray
@@ -16,9 +17,21 @@ class IModel(IDict):
     def mask_all_packages(self, mask: GridDataArray):
         raise NotImplementedError
 
+    @abstractmethod
+    def mask_packages(
+        self,
+        package_names: list[str],
+        mask: GridDataArray,
+        ignore_time_purge_empty: bool = False,
+    ):
+        raise NotImplementedError
+
     @abstractmethod
     def purge_empty_packages(
-        self, model_name: Optional[str] = "", ignore_time: bool = False
+        self,
+        model_name: Optional[str] = "",
+        ignore_time: bool = False,
+        package_names: list[str] | None = None,
     ) -> None:
         raise NotImplementedError
 
@@ -56,3 +69,8 @@ def _is_splitting_supported(self) -> Tuple[bool, str]:
     @abstractmethod
     def _is_clipping_supported(self) -> Tuple[bool, str]:
         raise NotImplementedError
+
+    @property
+    @abstractmethod
+    def _boundary_state_pkg_type(self) -> type[IPackage]:
+        raise NotImplementedError

imod/common/interfaces/itopsystembc.py

@@ -0,0 +1,15 @@
+from abc import abstractmethod
+
+from imod.common.interfaces.ipackage import IPackage
+from imod.typing import GridDataDict, GridDataset
+
+
+class ITopSystemBoundaryCondition(IPackage):
+    """
+    Interface for top system boundary condition packages in MODFLOW 6.
+    """
+
+    @classmethod
+    @abstractmethod
+    def aggregate_layers(cls, dataset: GridDataset) -> GridDataDict:
+        raise NotImplementedError

imod/common/utilities/mask.py

@@ -62,15 +62,18 @@ def mask_all_models(
             )
 
 
-def mask_all_packages(
+def mask_packages(
     model: IModel,
+    package_names: list[str],
     mask: GridDataArray,
     ignore_time_purge_empty: bool = False,
-):
+) -> None:
     _validate_coords_mask(mask)
-    for pkgname, pkg in model.items():
-        model[pkgname] = pkg.mask(mask)
-    model.purge_empty_packages(ignore_time=ignore_time_purge_empty)
+    for pkgname in package_names:
+        model[pkgname] = model[pkgname].mask(mask)
+    model.purge_empty_packages(
+        ignore_time=ignore_time_purge_empty, package_names=package_names
+    )
 
 
 def mask_package(package: IPackage, mask: GridDataArray) -> IPackage:

imod/common/utilities/regrid.py

@@ -27,7 +27,7 @@
     is_unstructured,
     ones_like,
 )
-from imod.util.dims import drop_layer_dim_cap_data, enforced_dim_order
+from imod.util.dims import drop_layer_dim_cap_and_bnd_data, enforced_dim_order
 from imod.util.regrid import (
     RegridderType,
     RegridderWeightsCache,
@@ -471,7 +471,7 @@ def _get_regridding_domain(
     return new_idomain
 
 
-def regrid_imod5_cap_data(
+def regrid_imod5_cap_and_bnd_data(
     imod5_data: Imod5DataDict,
     target_dis: IRegridPackage,
     regridder_types: DataclassType,
@@ -486,15 +486,19 @@ def regrid_imod5_cap_data(
     and ``imod.mf6.Recharge.from_imod5_cap_data``.
     """
     # Drop layer coords
-    imod5_cap_no_layer = drop_layer_dim_cap_data(imod5_data)
+    imod5_no_layer = drop_layer_dim_cap_and_bnd_data(imod5_data)
     target_grid = target_dis.dataset["idomain"].isel(layer=0, drop=True)
     # Regrid the input data
     cap_data_regridded = _regrid_package_data(
-        imod5_cap_no_layer["cap"], target_grid, regridder_types, regrid_cache
+        imod5_no_layer["cap"], target_grid, regridder_types, regrid_cache
+    )
+    bnd_data_regridded = _regrid_package_data(
+        imod5_no_layer["bnd"], target_grid, regridder_types, regrid_cache
     )
     extra_paths = imod5_data["extra"]["paths"]
     imod5_regridded: Imod5DataDict = {
         "cap": cap_data_regridded,
+        "bnd": bnd_data_regridded,
         "extra": {"paths": extra_paths},
     }
     return imod5_regridded

imod/mf6/boundary_condition.py

@@ -7,6 +7,7 @@
 import xarray as xr
 import xugrid as xu
 
+from imod.common.interfaces.iboundarycondition import IBoundaryCondition
 from imod.common.utilities.value_filters import enforce_scalar
 from imod.mf6.auxiliary_variables import (
     expand_transient_auxiliary_variables,
@@ -57,7 +58,7 @@ def _disv_recarr(arrdict, layer, notnull):
     return recarr
 
 
-class BoundaryCondition(Package, abc.ABC):
+class BoundaryCondition(Package, IBoundaryCondition, abc.ABC):
     """
     BoundaryCondition is used to share methods for specific stress packages
     with a time component.
@@ -324,7 +325,7 @@ def _get_period_varnames(self) -> list[str]:
         >>> river._get_period_varnames()
         >>> # prints: ['stage', 'conductance', 'bottom_elevation', 'species1', 'species2']
         """
-        result = []
+        result: list[str] = []
         if hasattr(self, "_period_data"):
             result.extend(self._period_data)
         if hasattr(self, "_optional_data"):