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/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

@@ -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_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/model.py

@@ -20,9 +20,8 @@
 from imod.common.interfaces.imodel import IModel
 from imod.common.serializer import EngineType
 from imod.common.statusinfo import NestedStatusInfo, StatusInfo, StatusInfoBase
-from imod.common.utilities.clip import clip_box_dataset
 from imod.common.utilities.dump_model import dump_model
-from imod.common.utilities.mask import mask_all_packages
+from imod.common.utilities.mask import mask_packages
 from imod.common.utilities.regrid import _regrid_like
 from imod.common.utilities.schemata import (
     concatenate_schemata_dicts,
@@ -40,9 +39,9 @@
 from imod.mf6.riv import River
 from imod.mf6.utilities.clipped_bc_creator import (
     StateClassType,
-    StateType,
-    create_clipped_boundary,
+    create_boundary_condition_clipped_boundary,
 )
+from imod.mf6.utilities.mask import mask_topsystem
 from imod.mf6.utilities.mf6hfb import merge_hfb_packages
 from imod.mf6.validation_settings import ValidationSettings
 from imod.mf6.wel import GridAgnosticWell
@@ -62,90 +61,6 @@ def pkg_has_cleanup(pkg: Package):
     return any(isinstance(pkg, pkgtype) for pkgtype in PKGTYPES_WITH_CLEANUP)
 
 
-def _create_boundary_condition_for_unassigned_boundary(
-    model: Modflow6Model,
-    state_for_boundary: Optional[GridDataArray],
-    additional_boundaries: list[Optional[StateType]] = [None],
-) -> Optional[StateType]:
-    if state_for_boundary is None:
-        return None
-
-    pkg_type = model._boundary_state_pkg_type
-    constant_state_packages = [
-        pkg for _, pkg in model.items() if isinstance(pkg, pkg_type)
-    ]
-
-    filtered_boundaries: list[StateType] = [
-        item for item in additional_boundaries or [] if item is not None
-    ]
-
-    constant_state_packages.extend(filtered_boundaries)
-
-    return create_clipped_boundary(
-        model.domain, state_for_boundary, constant_state_packages, pkg_type
-    )
-
-
-def _create_boundary_condition_clipped_boundary(
-    original_model: Modflow6Model,
-    clipped_model: Modflow6Model,
-    state_for_boundary: Optional[GridDataArray],
-    clip_box_args: tuple[Any, ...],
-) -> Optional[StateType]:
-    # Create temporary boundary condition for the original model boundary. This
-    # is used later to see which boundaries can be ignored as they were already
-    # present in the original model. We want to just end up with the boundary
-    # created by the clip.
-    unassigned_boundary_original_domain = (
-        _create_boundary_condition_for_unassigned_boundary(
-            original_model, state_for_boundary
-        )
-    )
-    # Clip the unassigned boundary to the clipped model's domain, required to
-    # avoid topological errors later.
-    if unassigned_boundary_original_domain is not None:
-        unassigned_boundary_clipped = unassigned_boundary_original_domain.clip_box(
-            *clip_box_args
-        )
-    else:
-        unassigned_boundary_clipped = None
-
-    if state_for_boundary is not None:
-        # Clip box as dataset, temporarily add variable name to convert to
-        # dataset, then turn back into DataArray.
-        varname = original_model._boundary_state_pkg_type._period_data[0]
-        state_for_boundary = state_for_boundary.to_dataset(name=varname)
-        state_for_boundary_clipped = clip_box_dataset(
-            state_for_boundary, *clip_box_args
-        )[varname]
-    else:
-        state_for_boundary_clipped = None
-
-    bc_constant_pkg = _create_boundary_condition_for_unassigned_boundary(
-        clipped_model, state_for_boundary_clipped, [unassigned_boundary_clipped]
-    )
-
-    # Remove all indices before first timestep of state_for_clipped_boundary.
-    # This to prevent empty dataarrays unnecessarily being made for these
-    # indices, which can lead to them to be removed when purging empty packages
-    # with ignore_time=True. Unfortunately, this is needs to be handled here and
-    # not in _create_boundary_condition_for_unassigned_boundary, as otherwise
-    # this function is called twice which could result in broadcasting errors in
-    # the second call if the time domain of state_for_boundary and assigned
-    # packages have no overlap.
-    if (
-        (state_for_boundary is not None)
-        and (state_for_boundary.indexes.get("time") is not None)
-        and (bc_constant_pkg is not None)
-    ):
-        start_time = state_for_boundary.indexes["time"][0]
-        bc_constant_pkg.dataset = bc_constant_pkg.dataset.sel(
-            time=slice(start_time, None)
-        )
-
-    return bc_constant_pkg
-
-
 class Modflow6Model(collections.UserDict[str, Package], IModel, abc.ABC):
     _mandatory_packages: tuple[str, ...] = ()
     _init_schemata: SchemataDict = {}
@@ -805,15 +720,29 @@ def clip_box(
             *clip_box_args,
         )
 
-        clipped_boundary_condition = _create_boundary_condition_clipped_boundary(
+        clipped_boundary_condition = create_boundary_condition_clipped_boundary(
             self, clipped, state_for_boundary, clip_box_args
         )
-        state_pkg_id = self._boundary_state_pkg_type._pkg_id
-        pkg_name = f"{state_pkg_id}_clipped"
         if clipped_boundary_condition is not None:
+            # Assign clipped boundary condition package
+            state_pkg_id = self._boundary_state_pkg_type._pkg_id
+            pkg_name = f"{state_pkg_id}_clipped"
+
             clipped[pkg_name] = clipped_boundary_condition
 
-        clipped.purge_empty_packages(ignore_time=ignore_time_purge_empty)
+            # Mask topsystem packages where the state boundary cells have been
+            # added.
+            state_varname = clipped_boundary_condition._period_data[0]
+            # Select the state variable for the first time step as mask.
+            # Its location will be constant through time.
+            state_var = clipped_boundary_condition.dataset[state_varname].isel(
+                time=0, missing_dims="ignore", drop=True
+            )
+            not_added_bc = np.isnan(state_var)
+            # Purge empty packages called by the mask_topsystem function
+            mask_topsystem(clipped, not_added_bc, ignore_time_purge_empty)
+        else:
+            clipped.purge_empty_packages(ignore_time=ignore_time_purge_empty)
 
         return clipped
 
@@ -930,11 +859,41 @@ def mask_all_packages(
             Whether to ignore time dimension when purging empty packages. Can
             improve performance when masking models with many time steps.
         """
+        package_names = list(self.keys())
+        mask_packages(self, package_names, mask, ignore_time_purge_empty)
+
+    def mask_packages(
+        self,
+        package_names: list[str],
+        mask: GridDataArray,
+        ignore_time_purge_empty: bool = False,
+    ) -> None:
+        """
+        This function applies a mask to packages in a model. The mask must
+        be presented as an idomain-like integer array that has 0 (inactive) or
+        <0 (vertical passthrough) values in filtered cells and >0 in active
+        cells.
+        Masking will overwrite idomain with the mask where the mask is <=0.
+        Where the mask is >0, the original value of idomain will be kept. Masking
+        will update the packages accordingly, blanking their input where needed,
+        and is therefore not a reversible operation.
 
-        mask_all_packages(self, mask, ignore_time_purge_empty)
+        Parameters
+        ----------
+        mask: xr.DataArray, xu.UgridDataArray of ints
+            idomain-like integer array. >0 sets cells to active, 0 sets cells to inactive,
+            <0 sets cells to vertical passthrough
+        ignore_time_purge_empty: bool, default False
+            Whether to ignore time dimension when purging empty packages. Can
+            improve performance when masking models with many time steps.
+        """
+        mask_packages(self, package_names, mask, ignore_time_purge_empty)
 
     def purge_empty_packages(
-        self, model_name: Optional[str] = "", ignore_time: bool = False
+        self,
+        model_name: Optional[str] = "",
+        ignore_time: bool = False,
+        package_names: Optional[list[str]] = None,
     ) -> None:
         """
         This method removes empty packages from the model in place.
@@ -948,11 +907,17 @@ def purge_empty_packages(
             timesteps. If True, packages are considered empty if they have no
             data at the first time step. The latter can increase performance
             considerably.
+        package_names: list[str], optional
+            List of package names to check for emptiness. If None, all packages
+            are checked.
         """
+        if package_names is None:
+            package_names = list(self.keys())
+
         empty_packages = [
             package_name
-            for package_name, package in self.items()
-            if package.is_empty(ignore_time=ignore_time)
+            for package_name in package_names
+            if self[package_name].is_empty(ignore_time=ignore_time)
         ]
         logger.info(
             f"packages: {empty_packages} removed in {model_name}, because all empty"

imod/mf6/model_gwf.py

@@ -37,6 +37,7 @@
 from imod.mf6.riv import River
 from imod.mf6.sto import StorageCoefficient
 from imod.mf6.utilities.chd_concat import concat_layered_chd_packages
+from imod.mf6.utilities.imod5_converter import mask_topsystem_packages_with_ibound
 from imod.mf6.validation_settings import ValidationSettings
 from imod.mf6.wel import LayeredWell, Well
 from imod.prepare.topsystem.default_allocation_methods import (
@@ -453,4 +454,12 @@ def from_imod5_data(
             for key, chd_package in chd_packages.items():
                 result[key] = chd_package
 
+        # Mask all topsystem packages where IBOUND < 0
+        mask_topsystem_packages_with_ibound(
+            imod5_data,
+            result,
+            cast(ConstantHeadRegridMethod, regridder_types.get("topsystem_mask")),
+            regrid_cache,
+            ignore_time_purge_empty=True,
+        )
         return result