diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst index 7db57326d..d36c88837 100644 --- a/docs/api/changelog.rst +++ b/docs/api/changelog.rst @@ -55,7 +55,10 @@ Changed :class:`imod.mf6.regrid.RiverRegridMethod`, :class:`imod.mf6.regrid.SpecificStorageRegridMethod`, :class:`imod.mf6.regrid.StorageCoefficientRegridMethod`. - +- Renamed ``imod.mf6.LayeredHorizontalFlowBarrier`` classes to + :class:`imod.mf6.SingleLayerHorizontalFlowBarrierResistance`, + :class:`imod.mf6.SingleLayerHorizontalFlowBarrierHydraulicCharacteristic`, + :class:`imod.mf6.SingleLayerHorizontalFlowBarrierMultiplier`, Fixed ~~~~~ diff --git a/imod/mf6/__init__.py b/imod/mf6/__init__.py index d72b1bdb8..d6225b065 100644 --- a/imod/mf6/__init__.py +++ b/imod/mf6/__init__.py @@ -21,9 +21,9 @@ HorizontalFlowBarrierHydraulicCharacteristic, HorizontalFlowBarrierMultiplier, HorizontalFlowBarrierResistance, - LayeredHorizontalFlowBarrierHydraulicCharacteristic, - LayeredHorizontalFlowBarrierMultiplier, - LayeredHorizontalFlowBarrierResistance, + SingleLayerHorizontalFlowBarrierHydraulicCharacteristic, + SingleLayerHorizontalFlowBarrierMultiplier, + SingleLayerHorizontalFlowBarrierResistance, ) from imod.mf6.ic import InitialConditions from imod.mf6.ims import ( diff --git a/imod/mf6/hfb.py b/imod/mf6/hfb.py index 45d4460c6..c2f51773f 100644 --- a/imod/mf6/hfb.py +++ b/imod/mf6/hfb.py @@ -4,11 +4,10 @@ import typing from copy import deepcopy from enum import Enum -from typing import TYPE_CHECKING, List, Optional, Tuple +from typing import TYPE_CHECKING, Dict, Optional, Tuple import cftime import numpy as np -import pandas as pd import xarray as xr import xugrid as xu from fastcore.dispatch import typedispatch @@ -19,7 +18,7 @@ from imod.mf6.mf6_hfb_adapter import Mf6HorizontalFlowBarrier from imod.mf6.package import Package from imod.mf6.utilities.grid import broadcast_to_full_domain -from imod.schemata import EmptyIndexesSchema +from imod.schemata import EmptyIndexesSchema, MaxNUniqueValuesSchema from imod.typing import GridDataArray from imod.util.imports import MissingOptionalModule @@ -343,7 +342,6 @@ def to_mf6_pkg( top: GridDataArray, bottom: GridDataArray, k: GridDataArray, - validate: bool = False, ) -> Mf6HorizontalFlowBarrier: """ Write package to Modflow 6 package. @@ -362,16 +360,11 @@ def to_mf6_pkg( Grid with bottom of model layers. k: GridDataArray Grid with hydraulic conductivities. - validate: bool - Run validation before converting Returns ------- """ - if validate: - self._validate(self._write_schemata) - top, bottom = broadcast_to_full_domain(idomain, top, bottom) k = idomain * k unstructured_grid, top, bottom, k = ( @@ -680,7 +673,7 @@ def __remove_edge_values_connected_to_inactive_cells( class HorizontalFlowBarrierHydraulicCharacteristic(HorizontalFlowBarrierBase): """ - Horizontal Flow Barrier (HFB) package + Horizontal Flow Barrier (HFB) package Input to the Horizontal Flow Barrier (HFB) Package is read from the file that has type "HFB6" in the Name File. Only one HFB Package can be @@ -741,9 +734,11 @@ def _compute_barrier_values( return barrier_values -class LayeredHorizontalFlowBarrierHydraulicCharacteristic(HorizontalFlowBarrierBase): +class SingleLayerHorizontalFlowBarrierHydraulicCharacteristic( + HorizontalFlowBarrierBase +): """ - Horizontal Flow Barrier (HFB) package + Horizontal Flow Barrier (HFB) package Input to the Horizontal Flow Barrier (HFB) Package is read from the file that has type "HFB6" in the Name File. Only one HFB Package can be @@ -755,8 +750,10 @@ class LayeredHorizontalFlowBarrierHydraulicCharacteristic(HorizontalFlowBarrierB geometry: gpd.GeoDataFrame Dataframe that describes: - geometry: the geometries of the barriers, - - hydraulic_characteristic: the hydraulic characteristic of the barriers - - layer: model layer for the barrier + - hydraulic_characteristic: the hydraulic characteristic of the + barriers + - layer: model layer for the barrier, only 1 single layer can be + entered. print_input: bool Examples @@ -775,6 +772,11 @@ class LayeredHorizontalFlowBarrierHydraulicCharacteristic(HorizontalFlowBarrierB """ + _write_schemata = { + "geometry": [EmptyIndexesSchema()], + "layer": [MaxNUniqueValuesSchema(1)], + } + @init_log_decorator() def __init__( self, @@ -806,7 +808,7 @@ def _compute_barrier_values( class HorizontalFlowBarrierMultiplier(HorizontalFlowBarrierBase): """ - Horizontal Flow Barrier (HFB) package + Horizontal Flow Barrier (HFB) package Input to the Horizontal Flow Barrier (HFB) Package is read from the file that has type "HFB6" in the Name File. Only one HFB Package can be @@ -872,24 +874,23 @@ def _compute_barrier_values( return barrier_values -class LayeredHorizontalFlowBarrierMultiplier(HorizontalFlowBarrierBase): +class SingleLayerHorizontalFlowBarrierMultiplier(HorizontalFlowBarrierBase): """ - Horizontal Flow Barrier (HFB) package + Horizontal Flow Barrier (HFB) package Input to the Horizontal Flow Barrier (HFB) Package is read from the file that has type "HFB6" in the Name File. Only one HFB Package can be specified for a GWF model. https://water.usgs.gov/water-resources/software/MODFLOW-6/mf6io_6.2.2.pdf - If parts of the barrier overlap a layer the multiplier is applied to the entire layer. - Parameters ---------- geometry: gpd.GeoDataFrame Dataframe that describes: - geometry: the geometries of the barriers, - multiplier: the multiplier of the barriers - - layer: model layer for the barrier + - layer: model layer for the barrier, only 1 single layer can be + entered. print_input: bool Examples @@ -908,6 +909,11 @@ class LayeredHorizontalFlowBarrierMultiplier(HorizontalFlowBarrierBase): """ + _write_schemata = { + "geometry": [EmptyIndexesSchema()], + "layer": [MaxNUniqueValuesSchema(1)], + } + @init_log_decorator() def __init__( self, @@ -1021,7 +1027,7 @@ def _compute_barrier_values( return barrier_values -class LayeredHorizontalFlowBarrierResistance(HorizontalFlowBarrierBase): +class SingleLayerHorizontalFlowBarrierResistance(HorizontalFlowBarrierBase): """ Horizontal Flow Barrier (HFB) package @@ -1036,7 +1042,8 @@ class LayeredHorizontalFlowBarrierResistance(HorizontalFlowBarrierBase): Dataframe that describes: - geometry: the geometries of the barriers, - resistance: the resistance of the barriers - - layer: model layer for the barrier + - layer: model layer for the barrier, only 1 single layer can be + entered. print_input: bool Examples @@ -1056,6 +1063,11 @@ class LayeredHorizontalFlowBarrierResistance(HorizontalFlowBarrierBase): """ + _write_schemata = { + "geometry": [EmptyIndexesSchema()], + "layer": [MaxNUniqueValuesSchema(1)], + } + @init_log_decorator() def __init__( self, @@ -1084,25 +1096,22 @@ def _compute_barrier_values( return barrier_values @classmethod - def from_imod5_dataset(cls, imod5_data: dict[str, dict[str, GridDataArray]]): + def from_imod5_dataset( + cls, key: str, imod5_data: Dict[str, Dict[str, GridDataArray]] + ): imod5_keys = list(imod5_data.keys()) - hfb_keys = [key for key in imod5_keys if key[0:3] == "hfb"] - if len(hfb_keys) == 0: - raise ValueError("no hfb keys present.") - - dataframe_ls: List[gpd.GeoDataFrame] = [] - for hfb_key in hfb_keys: - hfb_dict = imod5_data[hfb_key] - if not list(hfb_dict.keys()) == ["geodataframe", "layer"]: - raise ValueError("hfb is not a LayeredHorizontalFlowBarrierResistance") - layer = hfb_dict["layer"] - if layer == 0: - raise ValueError( - "assigning to layer 0 is not supported yet for import of HFB's" - ) - geometry_layer = hfb_dict["geodataframe"] - geometry_layer["layer"] = layer - dataframe_ls.append(geometry_layer) - compound_dataframe = pd.concat(dataframe_ls, ignore_index=True) + if key not in imod5_keys: + raise ValueError("hfb key not present.") + + hfb_dict = imod5_data[key] + if not list(hfb_dict.keys()) == ["geodataframe", "layer"]: + raise ValueError("hfb is not a SingleLayerHorizontalFlowBarrierResistance") + layer = hfb_dict["layer"] + if layer == 0: + raise ValueError( + "assigning to layer 0 is not supported yet for import of HFB's" + ) + geometry_layer = hfb_dict["geodataframe"] + geometry_layer["layer"] = layer - return LayeredHorizontalFlowBarrierResistance(compound_dataframe) + return cls(geometry_layer) diff --git a/imod/mf6/model.py b/imod/mf6/model.py index f3a5cc20a..2293425e8 100644 --- a/imod/mf6/model.py +++ b/imod/mf6/model.py @@ -274,7 +274,7 @@ def write( elif isinstance(pkg, imod.mf6.HorizontalFlowBarrierBase): top, bottom, idomain = self.__get_domain_geometry() k = self.__get_k() - mf6_hfb_pkg = pkg.to_mf6_pkg(idomain, top, bottom, k, validate) + mf6_hfb_pkg = pkg.to_mf6_pkg(idomain, top, bottom, k) mf6_hfb_pkg.write( pkgname=pkg_name, globaltimes=globaltimes, diff --git a/imod/mf6/model_gwf.py b/imod/mf6/model_gwf.py index be0b9eb69..6c4f73ae4 100644 --- a/imod/mf6/model_gwf.py +++ b/imod/mf6/model_gwf.py @@ -11,7 +11,7 @@ from imod.mf6.clipped_boundary_condition_creator import create_clipped_boundary from imod.mf6.dis import StructuredDiscretization from imod.mf6.drn import Drainage -from imod.mf6.hfb import LayeredHorizontalFlowBarrierResistance +from imod.mf6.hfb import SingleLayerHorizontalFlowBarrierResistance from imod.mf6.ic import InitialConditions from imod.mf6.model import Modflow6Model from imod.mf6.npf import NodePropertyFlow @@ -266,7 +266,11 @@ def from_imod5_data( # import hfb hfb_keys = [key for key in imod5_keys if key[0:3] == "hfb"] if len(hfb_keys) != 0: - hfb = LayeredHorizontalFlowBarrierResistance.from_imod5_dataset(imod5_data) - result["hfb"] = hfb + for hfb_key in hfb_keys: + result[hfb_key] = ( + SingleLayerHorizontalFlowBarrierResistance.from_imod5_dataset( + hfb_key, imod5_data + ) + ) return result diff --git a/imod/mf6/riv.py b/imod/mf6/riv.py index 50a90c605..6720bf17e 100644 --- a/imod/mf6/riv.py +++ b/imod/mf6/riv.py @@ -203,6 +203,9 @@ def from_imod5_data( Parameters ---------- + key: str + Packagename of the package that needs to be converted to river + package. imod5_data: dict Dictionary with iMOD5 data. This can be constructed from the :func:`imod.formats.prj.open_projectfile_data` method. diff --git a/imod/schemata.py b/imod/schemata.py index a61308d64..95b5d6a3a 100644 --- a/imod/schemata.py +++ b/imod/schemata.py @@ -503,6 +503,21 @@ def validate(self, obj: GridDataArray, **kwargs) -> None: ) +class MaxNUniqueValuesSchema(BaseSchema): + """ + Fails if amount of unique values exceeds a limit. + """ + + def __init__(self, max_unique_values: int): + self.max_unique_values = max_unique_values + + def validate(self, obj: GridDataArray, **kwargs) -> None: + if len(np.unique(obj)) > self.max_unique_values: + raise ValidationError( + f"Amount of unique values exceeds limit of {self.max_unique_values}" + ) + + class UniqueValuesSchema(BaseSchema): def __init__(self, other_value: list) -> None: """ diff --git a/imod/tests/test_mf6/test_mf6_hfb.py b/imod/tests/test_mf6/test_mf6_hfb.py index 81a242d16..eaf822ff0 100644 --- a/imod/tests/test_mf6/test_mf6_hfb.py +++ b/imod/tests/test_mf6/test_mf6_hfb.py @@ -12,9 +12,9 @@ HorizontalFlowBarrierHydraulicCharacteristic, HorizontalFlowBarrierMultiplier, HorizontalFlowBarrierResistance, - LayeredHorizontalFlowBarrierHydraulicCharacteristic, - LayeredHorizontalFlowBarrierMultiplier, - LayeredHorizontalFlowBarrierResistance, + SingleLayerHorizontalFlowBarrierHydraulicCharacteristic, + SingleLayerHorizontalFlowBarrierMultiplier, + SingleLayerHorizontalFlowBarrierResistance, ) from imod.mf6.dis import StructuredDiscretization from imod.mf6.hfb import to_connected_cells_dataset @@ -153,10 +153,10 @@ def test_to_mf6_creates_mf6_adapter( @pytest.mark.parametrize( "barrier_class, barrier_value_name, barrier_value, expected_hydraulic_characteristic", [ - (LayeredHorizontalFlowBarrierResistance, "resistance", 1e3, 1e-3), - (LayeredHorizontalFlowBarrierMultiplier, "multiplier", 1.5, -1.5), + (SingleLayerHorizontalFlowBarrierResistance, "resistance", 1e3, 1e-3), + (SingleLayerHorizontalFlowBarrierMultiplier, "multiplier", 1.5, -1.5), ( - LayeredHorizontalFlowBarrierHydraulicCharacteristic, + SingleLayerHorizontalFlowBarrierHydraulicCharacteristic, "hydraulic_characteristic", 1e-3, 1e-3, @@ -195,7 +195,7 @@ def test_to_mf6_creates_mf6_adapter_layered( hfb = barrier_class(geometry, print_input) # Act. - _ = hfb.to_mf6_pkg(idomain, top, bottom, k, False) + _ = hfb.to_mf6_pkg(idomain, top, bottom, k) # Assert. snapped, _ = xu.snap_to_grid(geometry, grid=idomain, max_snap_distance=0.5) @@ -282,6 +282,69 @@ def test_to_mf6_different_z_boundaries( assert_array_equal(max_values_per_layer, 1.0 / expected_values) +@pytest.mark.parametrize( + "layer, expected_values", + [ + (2, np.array([1e3, 1e3, 1e3, 1e3, 1e3, 1e3, 1e3, 1e3])), # 2nd layer + ], +) +@patch("imod.mf6.mf6_hfb_adapter.Mf6HorizontalFlowBarrier.__new__", autospec=True) +def test_to_mf6_layered_hfb(mf6_flow_barrier_mock, basic_dis, layer, expected_values): + # Arrange. + idomain, top, bottom = basic_dis + k = ones_like(top) + + print_input = False + + barrier_y = [5.5, 5.5, 5.5] + barrier_x = [82.0, 40.0, 0.0] + + geometry = gpd.GeoDataFrame( + geometry=[shapely.linestrings(barrier_x, barrier_y)], + data={ + "resistance": [1e3], + "layer": [layer], + }, + ) + + hfb = SingleLayerHorizontalFlowBarrierResistance(geometry, print_input) + + # Act. + _ = hfb.to_mf6_pkg(idomain, top, bottom, k) + + # Assert. + _, args = mf6_flow_barrier_mock.call_args + barrier_values = args["hydraulic_characteristic"].values + assert_array_equal(barrier_values, 1.0 / expected_values) + expected_layer = np.full((8,), layer) + barrier_layer = args["layer"].values + assert_array_equal(barrier_layer, expected_layer) + + +def test_to_mf6_layered_hfb__error(): + """Throws error because multiple layers attached to one object.""" + # Arrange. + print_input = False + + barrier_y = [5.5, 5.5, 5.5] + barrier_x = [82.0, 40.0, 0.0] + + linestring = shapely.linestrings(barrier_x, barrier_y) + + geometry = gpd.GeoDataFrame( + geometry=[linestring, linestring], + data={ + "resistance": [1e3, 1e3], + "layer": [1, 2], + }, + ) + + hfb = SingleLayerHorizontalFlowBarrierResistance(geometry, print_input) + errors = hfb._validate(hfb._write_schemata) + + assert len(errors) > 0 + + @pytest.mark.parametrize( "barrier_x_loc, expected_number_barriers", [ @@ -451,8 +514,10 @@ def test_hfb_from_imod5(imod5_dataset, tmp_path): imod5_dataset, target_dis.dataset["idomain"] ) - hfb = LayeredHorizontalFlowBarrierResistance.from_imod5_dataset(imod5_dataset) + hfb = SingleLayerHorizontalFlowBarrierResistance.from_imod5_dataset( + "hfb-3", imod5_dataset + ) hfb_package = hfb.to_mf6_pkg( target_dis["idomain"], target_dis["top"], target_dis["bottom"], target_npf["k"] ) - assert list(np.unique(hfb_package["layer"].values)) == [3, 21] + assert list(np.unique(hfb_package["layer"].values)) == [7] diff --git a/imod/tests/test_mf6/test_mf6_simulation.py b/imod/tests/test_mf6/test_mf6_simulation.py index 3e9dda2dd..dccba0b04 100644 --- a/imod/tests/test_mf6/test_mf6_simulation.py +++ b/imod/tests/test_mf6/test_mf6_simulation.py @@ -485,5 +485,10 @@ def test_import_from_imod5(imod5_dataset, tmp_path): simulation.create_time_discretization(["01-01-2003", "02-01-2003"]) + # Remove HFB packages outside domain + # TODO: Build in support for hfb packages outside domain + for hfb_outside in ["hfb-24", "hfb-26"]: + simulation["imported_model"].pop(hfb_outside) + # write and validate the simulation. simulation.write(tmp_path, binary=False, validate=True)