diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst index 330bf8144..9c067eedd 100644 --- a/docs/api/changelog.rst +++ b/docs/api/changelog.rst @@ -14,6 +14,9 @@ Added - :class:`imod.msw.MeteoGridCopy` to copy existing `mete_grid.inp` files, so ASCII grids in large existing meteo databases do not have to be read. +- :meth:`imod.mf6.LayeredWell.from_imod5_cap_data` to construct a + :class:`imod.mf6.LayeredWell` package from iMOD5 data in the CAP package (for + MetaSWAP). Currently only griddata (IDF) is supported. - :meth:`imod.mf6.Recharge.from_imod5_cap_data` to construct a recharge package for coupling a MODFLOW6 model to MetaSWAP. @@ -35,6 +38,9 @@ Changed ``downward_resistance`` now require a ``subunit`` coordinate. - Variables ``max_abstraction_groundwater`` and ``max_abstraction_surfacewater`` in :class:`imod.msw.Sprinkling` now needs to have a subunit coordinate. +- If ``"cap"`` package present in ``imod5_data``, + :meth:`imod.mf6.GroundwaterFlowModel.from_imod5_data` now automatically adds a + well for metaswap sprinkling named ``"msw-sprinkling"`` [0.18.1] - 2024-11-20 diff --git a/docs/api/mf6.rst b/docs/api/mf6.rst index 4d43ead80..e7683632e 100644 --- a/docs/api/mf6.rst +++ b/docs/api/mf6.rst @@ -97,6 +97,7 @@ Flow Packages HorizontalFlowBarrierResistance.to_mf6_pkg LayeredWell LayeredWell.from_imod5_data + LayeredWell.from_imod5_cap_data LayeredWell.mask LayeredWell.regrid_like LayeredWell.to_mf6_pkg diff --git a/imod/mf6/model_gwf.py b/imod/mf6/model_gwf.py index 9ae5ece06..e83ee5e1c 100644 --- a/imod/mf6/model_gwf.py +++ b/imod/mf6/model_gwf.py @@ -303,9 +303,9 @@ def from_imod5_data( ) # now import the non-singleton packages' + imod5_keys = list(imod5_data.keys()) # import wells - imod5_keys = list(imod5_data.keys()) wel_keys = [key for key in imod5_keys if key[0:3] == "wel"] for wel_key in wel_keys: wel_key_truncated = wel_key[:16] @@ -331,10 +331,10 @@ def from_imod5_data( ) if "cap" in imod5_keys: + result["msw-sprinkling"] = LayeredWell.from_imod5_cap_data(imod5_data) # type: ignore result["msw-rch"] = Recharge.from_imod5_cap_data(imod5_data) # type: ignore # import ghb's - imod5_keys = list(imod5_data.keys()) ghb_keys = [key for key in imod5_keys if key[0:3] == "ghb"] for ghb_key in ghb_keys: ghb_pkg = GeneralHeadBoundary.from_imod5_data( diff --git a/imod/mf6/utilities/imod5_converter.py b/imod/mf6/utilities/imod5_converter.py index e599c8522..c3301430f 100644 --- a/imod/mf6/utilities/imod5_converter.py +++ b/imod/mf6/utilities/imod5_converter.py @@ -1,8 +1,10 @@ -from typing import Union +from typing import Union, cast import numpy as np +import pandas as pd import xarray as xr +from imod.typing import GridDataDict, Imod5DataDict, SelSettingsType from imod.typing.grid import full_like @@ -48,3 +50,72 @@ def fill_missing_layers( """ layer = full.coords["layer"] return source.reindex(layer=layer, fill_value=fillvalue) + + +def _well_from_imod5_cap_point_data(cap_data: GridDataDict) -> dict[str, np.ndarray]: + raise NotImplementedError( + "Assigning sprinkling wells with an IPF file is not supported, please specify them as IDF." + ) + return {} + + +def _well_from_imod5_cap_grid_data(cap_data: GridDataDict) -> dict[str, np.ndarray]: + drop_layer_kwargs: SelSettingsType = { + "layer": 0, + "drop": True, + "missing_dims": "ignore", + } + type = cap_data["artificial_recharge"].isel(**drop_layer_kwargs).compute() + layer = ( + cap_data["artificial_recharge_layer"] + .isel(**drop_layer_kwargs) + .astype(int) + .compute() + ) + + from_groundwater = (type != 0).to_numpy() + coords = type.coords + x_grid, y_grid = np.meshgrid(coords["x"].to_numpy(), coords["y"].to_numpy()) + + data = {} + data["layer"] = layer.data[from_groundwater] + data["y"] = y_grid[from_groundwater] + data["x"] = x_grid[from_groundwater] + data["rate"] = np.zeros_like(data["x"]) + + return data + + +def well_from_imod5_cap_data(imod5_data: Imod5DataDict) -> dict[str, np.ndarray]: + """ + Abstraction data for sprinkling is defined in iMOD5 either with grids (IDF) + or points (IPF) combined with a grid. Depending on the type, the function does + different conversions + + - grids (IDF) + The ``"artifical_recharge_layer"`` variable was defined as grid + (IDF), this grid defines in which layer a groundwater abstraction + well should be placed. The ``"artificial_recharge"`` grid contains + types which point to the type of abstraction: + * 0: no abstraction + * 1: groundwater abstraction + * 2: surfacewater abstraction + The ``"artificial_recharge_capacity"`` grid/constant defines the + capacity of each groundwater or surfacewater abstraction. This is an + ``1:1`` mapping: Each grid cell maps to a separate well. + + - points with grid (IPF & IDF) + The ``"artifical_recharge_layer"`` variable was defined as point + data (IPF), this table contains wellids with an abstraction capacity + and layer. The ``"artificial_recharge"`` grid contains a mapping of + grid cells to wellids in the point data. The + ``"artificial_recharge_capacity"`` is ignored as the abstraction + capacity is already defined in the point data. This is an ``n:1`` + mapping: multiple grid cells can map to one well. + """ + cap_data = cast(GridDataDict, imod5_data["cap"]) + + if isinstance(cap_data["artificial_recharge_layer"], pd.DataFrame): + return _well_from_imod5_cap_point_data(cap_data) + else: + return _well_from_imod5_cap_grid_data(cap_data) diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py index 372cb9f48..d30b58d30 100644 --- a/imod/mf6/wel.py +++ b/imod/mf6/wel.py @@ -30,6 +30,7 @@ from imod.mf6.package import Package from imod.mf6.utilities.dataset import remove_inactive from imod.mf6.utilities.grid import broadcast_to_full_domain +from imod.mf6.utilities.imod5_converter import well_from_imod5_cap_data from imod.mf6.validation import validation_pkg_error_message from imod.mf6.validation_context import ValidationContext from imod.mf6.write_context import WriteContext @@ -44,7 +45,7 @@ ValidationError, ) from imod.select.points import points_indices, points_values -from imod.typing import GridDataArray +from imod.typing import GridDataArray, Imod5DataDict from imod.typing.grid import is_spatial_grid, ones_like from imod.util.expand_repetitions import resample_timeseries from imod.util.structured import values_within_range @@ -1298,6 +1299,46 @@ def _validate_imod5_depth_information( logger.log(loglevel=LogLevel.ERROR, message=log_msg, additional_depth=2) raise ValueError(log_msg) + @classmethod + def from_imod5_cap_data(cls, imod5_data: Imod5DataDict): + """ + Create LayeredWell from imod5_data in "cap" package. Abstraction data + for sprinkling is defined in iMOD5 either with grids (IDF) or points + (IPF) combined with a grid. Depending on the type, the function does + different conversions + + - grids (IDF) + The ``"artifical_recharge_layer"`` variable was defined as grid + (IDF), this grid defines in which layer a groundwater abstraction + well should be placed. The ``"artificial_recharge"`` grid contains + types which point to the type of abstraction: + * 0: no abstraction + * 1: groundwater abstraction + * 2: surfacewater abstraction + The ``"artificial_recharge_capacity"`` grid/constant defines the + capacity of each groundwater or surfacewater abstraction. This is an + ``1:1`` mapping: Each grid cell maps to a separate well. + + - points with grid (IPF & IDF) + The ``"artifical_recharge_layer"`` variable was defined as point + data (IPF), this table contains wellids with an abstraction capacity + and layer. The ``"artificial_recharge"`` grid contains a mapping of + grid cells to wellids in the point data. The + ``"artificial_recharge_capacity"`` is ignored as the abstraction + capacity is already defined in the point data. This is an ``n:1`` + mapping: multiple grid cells can map to one well. + + Parameters + ---------- + imod5_data: dict[str, dict[str, GridDataArray]] + dictionary containing the arrays mentioned in the project file as + xarray datasets, under the key of the package type to which it + belongs, as returned by + :func:`imod.formats.prj.open_projectfile_data`. + """ + data = well_from_imod5_cap_data(imod5_data) + return cls(**data) # type: ignore + class WellDisStructured(DisStructuredBoundaryCondition): """ diff --git a/imod/msw/sprinkling.py b/imod/msw/sprinkling.py index 5d8bd20f1..ea9925112 100644 --- a/imod/msw/sprinkling.py +++ b/imod/msw/sprinkling.py @@ -1,4 +1,4 @@ -from typing import TextIO +from typing import TextIO, cast import numpy as np import pandas as pd @@ -10,7 +10,9 @@ from imod.msw.fixed_format import VariableMetaData from imod.msw.pkgbase import MetaSwapPackage from imod.msw.regrid.regrid_schemes import SprinklingRegridMethod -from imod.typing import IntArray +from imod.msw.utilities.common import concat_imod5 +from imod.typing import GridDataDict, Imod5DataDict, IntArray, SelSettingsType +from imod.typing.grid import zeros_like def _ravel_per_subunit(da: xr.DataArray) -> np.ndarray: @@ -20,6 +22,49 @@ def _ravel_per_subunit(da: xr.DataArray) -> np.ndarray: return array_out[np.isfinite(array_out)] +def _sprinkling_data_from_imod5_ipf(cap_data: GridDataDict) -> GridDataDict: + raise NotImplementedError( + "Assigning sprinkling wells with an IPF file is not supported, please specify them as IDF." + ) + return {} + + +def _sprinkling_data_from_imod5_grid(cap_data: GridDataDict) -> GridDataDict: + drop_layer_kwargs: SelSettingsType = { + "layer": 0, + "drop": True, + "missing_dims": "ignore", + } + type = cap_data["artificial_recharge"].isel(**drop_layer_kwargs) + capacity = cap_data["artificial_recharge_capacity"].isel(**drop_layer_kwargs) + + from_groundwater = type == 1 + from_surfacewater = type == 2 + is_active = type != 0 + + zero_where_active = zeros_like(type).where(is_active) + + # Add zero where active, to have active cells set to 0.0. + max_abstraction_groundwater_rural = zero_where_active.where( + ~from_groundwater, capacity + ) + max_abstraction_surfacewater_rural = zero_where_active.where( + ~from_surfacewater, capacity + ) + + # No sprinkling for urban environments + max_abstraction_urban = zero_where_active + + data = {} + data["max_abstraction_groundwater"] = concat_imod5( + max_abstraction_groundwater_rural, max_abstraction_urban + ) + data["max_abstraction_surfacewater"] = concat_imod5( + max_abstraction_surfacewater_rural, max_abstraction_urban + ) + return data + + class Sprinkling(MetaSwapPackage, IRegridPackage): """ This contains the sprinkling capacities of links between SVAT units and @@ -122,3 +167,51 @@ def _render( self._check_range(dataframe) return self.write_dataframe_fixed_width(file, dataframe) + + @classmethod + def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "Sprinkling": + """ + Import sprinkling data from imod5 data. Abstraction data for sprinkling + is defined in iMOD5 either with grids (IDF) or points (IPF) combined + with a grid. Depending on the type, the method does different conversions: + + - grids (IDF) + The ``"artifical_recharge_layer"`` variable was defined as grid + (IDF), this grid defines in which layer a groundwater abstraction + well should be placed. The ``"artificial_recharge"`` grid contains + types which point to the type of abstraction: + * 0: no abstraction + * 1: groundwater abstraction + * 2: surfacewater abstraction + The ``"artificial_recharge_capacity"`` grid/constant defines the + capacity of each groundwater or surfacewater abstraction. This is an + ``1:1`` mapping: Each grid cell maps to a separate well. + + - points with grid (IPF & IDF) + The ``"artifical_recharge_layer"`` variable was defined as point + data (IPF), this table contains wellids with an abstraction capacity + and layer. The ``"artificial_recharge"`` grid contains a mapping of + grid cells to wellids in the point data. The + ``"artificial_recharge_capacity"`` is ignored as the abstraction + capacity is already defined in the point data. This is an ``n:1`` + mapping: multiple grid cells can map to one well. + + Parameters + ---------- + imod5_data: dict[str, dict[str, GridDataArray]] + dictionary containing the arrays mentioned in the project file as + xarray datasets, under the key of the package type to which it + belongs, as returned by + :func:`imod.formats.prj.open_projectfile_data`. + + Returns + ------- + Sprinkling package + """ + cap_data = cast(GridDataDict, imod5_data["cap"]) + if isinstance(cap_data["artificial_recharge_layer"], pd.DataFrame): + data = _sprinkling_data_from_imod5_ipf(cap_data) + else: + data = _sprinkling_data_from_imod5_grid(cap_data) + + return cls(**data) diff --git a/imod/tests/conftest.py b/imod/tests/conftest.py index 593ffb5bb..0677d4b4f 100644 --- a/imod/tests/conftest.py +++ b/imod/tests/conftest.py @@ -22,6 +22,10 @@ ) from .fixtures.flow_example_fixture import imodflow_model from .fixtures.flow_transport_simulation_fixture import flow_transport_simulation +from .fixtures.imod5_cap_data import ( + cap_data_sprinkling_grid, + cap_data_sprinkling_points, +) from .fixtures.imod5_well_data import ( well_duplication_import_prj, well_mixed_ipfs, diff --git a/imod/tests/fixtures/imod5_cap_data.py b/imod/tests/fixtures/imod5_cap_data.py new file mode 100644 index 000000000..c32749eca --- /dev/null +++ b/imod/tests/fixtures/imod5_cap_data.py @@ -0,0 +1,60 @@ +import numpy as np +import pandas as pd +import pytest +import xarray as xr + +from imod.typing import Imod5DataDict + + +def zeros_grid(): + x = [1.0, 2.0, 3.0] + y = [3.0, 2.0, 1.0] + dx = 1.0 + dy = -1.0 + + coords = {"x": x, "y": y, "dx": dx, "dy": dy} + shape = (len(y), len(x)) + data = np.zeros(shape) + + return xr.DataArray(data, coords=coords, dims=("y", "x")) + + +@pytest.fixture(scope="function") +def cap_data_sprinkling_grid() -> Imod5DataDict: + type = zeros_grid() + type[:, 1] = 1 + type[:, 2] = 2 + layer = xr.ones_like(type) + layer[:, 1] = 2 + + cap_data = { + "artificial_recharge": type, + "artificial_recharge_layer": layer, + "artificial_recharge_capacity": xr.DataArray(25.0), + } + + return {"cap": cap_data} + + +@pytest.fixture(scope="function") +def cap_data_sprinkling_points() -> Imod5DataDict: + type = zeros_grid() + type[:, 1] = 3000 + type[:, 2] = 4000 + + data = { + "id": [3000, 4000], + "layer": [2, 3], + "capacity": [15.0, 30.0], + "y": [1.0, 2.0], + "x": [1.0, 2.0], + } + + layer = pd.DataFrame(data=data) + cap_data = { + "artificial_recharge": type, + "artificial_recharge_layer": layer, + "artificial_recharge_capacity": xr.DataArray(25.0), + } + + return {"cap": cap_data} diff --git a/imod/tests/test_mf6/test_mf6_wel.py b/imod/tests/test_mf6/test_mf6_wel.py index 8fc66d050..91475cef7 100644 --- a/imod/tests/test_mf6/test_mf6_wel.py +++ b/imod/tests/test_mf6/test_mf6_wel.py @@ -1176,3 +1176,25 @@ def test_logmessage_for_missing_filter_settings( in log ) assert message_required == message_present + + +def test_from_imod5_cap_data__grid(cap_data_sprinkling_grid): + # Arrange + expected_layer = np.array([2, 1, 2, 1, 2, 1]) + expected_y = np.array([3.0, 3.0, 2.0, 2.0, 1.0, 1.0]) + expected_x = np.array([2.0, 3.0, 2.0, 3.0, 2.0, 3.0]) + + # Act + well = LayeredWell.from_imod5_cap_data(cap_data_sprinkling_grid) + + # Assert + ds = well.dataset + np.testing.assert_equal(ds["rate"].to_numpy(), 0.0) + np.testing.assert_array_equal(ds["layer"].to_numpy(), expected_layer) + np.testing.assert_array_equal(ds["x"].to_numpy(), expected_x) + np.testing.assert_array_equal(ds["y"].to_numpy(), expected_y) + + +def test_from_imod5_cap_data__points(cap_data_sprinkling_points): + with pytest.raises(NotImplementedError): + LayeredWell.from_imod5_cap_data(cap_data_sprinkling_points) diff --git a/imod/tests/test_msw/test_sprinkling.py b/imod/tests/test_msw/test_sprinkling.py index 3954fcbee..74714d3a4 100644 --- a/imod/tests/test_msw/test_sprinkling.py +++ b/imod/tests/test_msw/test_sprinkling.py @@ -2,6 +2,7 @@ from pathlib import Path import numpy as np +import pytest import xarray as xr from numpy import nan from numpy.testing import assert_almost_equal, assert_equal @@ -269,3 +270,39 @@ def test_simple_model_1_subunit(fixed_format_parser): ) assert_equal(results["layer"], np.array([3, 2])) assert_equal(results["svat_groundwater"], np.array([1, 2])) + + +def test_sprinkling_from_imod5_data__points(cap_data_sprinkling_points): + with pytest.raises(NotImplementedError): + msw.Sprinkling.from_imod5_data(cap_data_sprinkling_points) + + +def test_sprinkling_from_imod5_data__grid(cap_data_sprinkling_grid): + # Arrange + # fmt: off + expected_gw_abstraction = np.array( + [[nan, 25., 0.], + [nan, 25., 0.], + [nan, 25., 0.]] + ) + expected_sw_abstraction = np.array( + [[nan, 0., 25.], + [nan, 0., 25.], + [nan, 0., 25.]] + ) + # fmt: on + + # Act + sprinkling = msw.Sprinkling.from_imod5_data(cap_data_sprinkling_grid) + + # Assert + assert isinstance(sprinkling, msw.Sprinkling) + ds = sprinkling.dataset + assert (ds.sel(subunit=1) == 0).all() + rural_ds = ds.sel(subunit=0) + np.testing.assert_array_equal( + rural_ds["max_abstraction_groundwater"].to_numpy(), expected_gw_abstraction + ) + np.testing.assert_array_equal( + rural_ds["max_abstraction_surfacewater"].to_numpy(), expected_sw_abstraction + ) diff --git a/imod/typing/__init__.py b/imod/typing/__init__.py index 8beb7fd13..0217f9782 100644 --- a/imod/typing/__init__.py +++ b/imod/typing/__init__.py @@ -2,7 +2,7 @@ Module to define type aliases. """ -from typing import TYPE_CHECKING, TypeAlias, TypeVar, Union +from typing import TYPE_CHECKING, Literal, TypeAlias, TypedDict, TypeVar, Union import numpy as np import xarray as xr @@ -20,6 +20,12 @@ IntArray: TypeAlias = NDArray[np.int_] +class SelSettingsType(TypedDict, total=False): + layer: int + drop: bool + missing_dims: Literal["raise", "warn", "ignore"] + + # Types for optional dependencies. if TYPE_CHECKING: import geopandas as gpd