From 4d207e062b5363ce742ba3cb35bb5145c5505749 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Fri, 29 Nov 2024 17:12:09 +0100 Subject: [PATCH 01/10] Start working on import --- imod/mf6/rch.py | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/imod/mf6/rch.py b/imod/mf6/rch.py index ba917cbd8..b202c57ac 100644 --- a/imod/mf6/rch.py +++ b/imod/mf6/rch.py @@ -218,7 +218,7 @@ def from_imod5_data( # create an array indicating in which cells rch is active is_rch_cell = allocate_rch_cells( ALLOCATION_OPTION.at_first_active, - new_idomain == 1, + new_idomain > 0, planar_rate_regridded, ) @@ -229,6 +229,23 @@ def from_imod5_data( return Recharge(**new_package_data, validate=True, fixed_cell=False) + @classmethod + def from_imod5_cap_data( + cls, + imod5_data: dict[str, dict[str, GridDataArray]], + target_dis: StructuredDiscretization, + ) -> "Recharge": + """ + Construct an rch-package from iMOD5 data in the CAP package, loaded with + the :func:`imod.formats.prj.open_projectfile_data` function. + """ + cap_data = imod5_data["cap"] + new_idomain = target_dis.dataset["idomain"] + + # Find which cells are active, basically same logic as in + # msw.GridData.from_imod5_data. Maybe cache that somewhere? + + @classmethod def get_regrid_methods(cls) -> RechargeRegridMethod: return deepcopy(cls._regrid_method) From a3e5d98eb6bce4fbf14c904d44241ba17aca7fd5 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Mon, 2 Dec 2024 10:27:17 +0100 Subject: [PATCH 02/10] Move util functions to imod5_converter utilities --- imod/msw/grid_data.py | 71 +++------------------- imod/msw/utilities/imod5_converter.py | 87 +++++++++++++++++++++++++++ 2 files changed, 96 insertions(+), 62 deletions(-) create mode 100644 imod/msw/utilities/imod5_converter.py diff --git a/imod/msw/grid_data.py b/imod/msw/grid_data.py index 3e08a69b0..ad7299206 100644 --- a/imod/msw/grid_data.py +++ b/imod/msw/grid_data.py @@ -3,7 +3,6 @@ import numpy as np import xarray as xr -from imod.logging import LogLevel, logger from imod.mf6 import StructuredDiscretization from imod.mf6.interfaces.iregridpackage import IRegridPackage from imod.mf6.regrid.regrid_schemes import ( @@ -13,65 +12,16 @@ from imod.msw.fixed_format import VariableMetaData from imod.msw.pkgbase import MetaSwapPackage from imod.msw.regrid.regrid_schemes import GridDataRegridMethod -from imod.msw.utilities.common import concat_imod5 -from imod.typing import GridDataArray, GridDataDict -from imod.typing.grid import ones_like +from imod.msw.utilities.imod5_converter import ( + get_cell_area_from_imod5_data, + get_landuse_from_imod5_data, + get_rootzone_depth_from_imod5_data, + is_msw_active_cell, +) +from imod.typing import GridDataDict from imod.util.spatial import get_cell_area, spatial_reference -def get_cell_area_from_imod5_data( - imod5_cap: GridDataDict, -) -> GridDataArray: - # area's per type of svats - mf6_area = get_cell_area(imod5_cap["wetted_area"]) - wetted_area = imod5_cap["wetted_area"] - urban_area = imod5_cap["urban_area"] - rural_area = mf6_area - (wetted_area + urban_area) - if (wetted_area > mf6_area).any(): - logger.log( - loglevel=LogLevel.WARNING, - message=f"wetted area was set to the max cell area of {mf6_area}", - additional_depth=0, - ) - wetted_area = wetted_area.where(wetted_area <= mf6_area, other=mf6_area) - if (rural_area < 0.0).any(): - logger.log( - loglevel=LogLevel.WARNING, - message="found urban area > than (cel-area - wetted area). Urban area was set to 0", - additional_depth=0, - ) - urban_area = urban_area.where(rural_area > 0.0, other=0.0) - rural_area = mf6_area - (wetted_area + urban_area) - return concat_imod5(rural_area, urban_area) - - -def get_landuse_from_imod5_data( - imod5_cap: GridDataDict, -) -> GridDataArray: - """ - Get landuse from imod5 capillary zone data. This adds two subunits, one - based on the landuse grid, which specifies rural landuse. The other - specifies urban landuse, which is coded to value 18. - """ - rural_landuse = imod5_cap["landuse"] - # Urban landuse = 18 - urban_landuse = ones_like(rural_landuse) * 18 - return concat_imod5(rural_landuse, urban_landuse) - - -def get_rootzone_depth_from_imod5_data( - imod5_cap: GridDataDict, -) -> GridDataArray: - """ - Get rootzone depth from imod5 capillary zone data. Also does a unit - conversion: iMOD5 specifies rootzone thickness in centimeters, whereas - MetaSWAP requires rootzone depth in meters. - """ - rootzone_thickness = imod5_cap["rootzone_thickness"] * 0.01 - # rootzone depth is equal for both svats. - return concat_imod5(rootzone_thickness, rootzone_thickness) - - class GridData(MetaSwapPackage, IRegridPackage): """ This contains the grid data of MetaSWAP. @@ -194,11 +144,8 @@ def from_imod5_data( data["surface_elevation"] = imod5_cap["surface_elevation"] data["soil_physical_unit"] = imod5_cap["soil_physical_unit"] - mf6_top_active = target_dis["idomain"].isel(layer=0, drop=True) - subunit_active = ( - (imod5_cap["boundary"] == 1) & (data["area"] > 0) & (mf6_top_active >= 1) - ) - active = subunit_active.all(dim="subunit") + active, subunit_active = is_msw_active_cell(target_dis, imod5_cap, data["area"]) + data_active = { key: ( griddata.where(subunit_active) diff --git a/imod/msw/utilities/imod5_converter.py b/imod/msw/utilities/imod5_converter.py new file mode 100644 index 000000000..70f88c982 --- /dev/null +++ b/imod/msw/utilities/imod5_converter.py @@ -0,0 +1,87 @@ +from imod.logging import LogLevel, logger +from imod.mf6 import StructuredDiscretization +from imod.msw.utilities.common import concat_imod5 +from imod.typing import GridDataArray, GridDataDict +from imod.typing.grid import ones_like +from imod.util.spatial import get_cell_area + + +def get_cell_area_from_imod5_data( + imod5_cap: GridDataDict, +) -> GridDataArray: + # area's per type of svats + mf6_area = get_cell_area(imod5_cap["wetted_area"]) + wetted_area = imod5_cap["wetted_area"] + urban_area = imod5_cap["urban_area"] + rural_area = mf6_area - (wetted_area + urban_area) + if (wetted_area > mf6_area).any(): + logger.log( + loglevel=LogLevel.WARNING, + message=f"wetted area was set to the max cell area of {mf6_area}", + additional_depth=0, + ) + wetted_area = wetted_area.where(wetted_area <= mf6_area, other=mf6_area) + if (rural_area < 0.0).any(): + logger.log( + loglevel=LogLevel.WARNING, + message="found urban area > than (cel-area - wetted area). Urban area was set to 0", + additional_depth=0, + ) + urban_area = urban_area.where(rural_area > 0.0, other=0.0) + rural_area = mf6_area - (wetted_area + urban_area) + return concat_imod5(rural_area, urban_area) + + +def get_landuse_from_imod5_data( + imod5_cap: GridDataDict, +) -> GridDataArray: + """ + Get landuse from imod5 capillary zone data. This adds two subunits, one + based on the landuse grid, which specifies rural landuse. The other + specifies urban landuse, which is coded to value 18. + """ + rural_landuse = imod5_cap["landuse"] + # Urban landuse = 18 + urban_landuse = ones_like(rural_landuse) * 18 + return concat_imod5(rural_landuse, urban_landuse) + + +def get_rootzone_depth_from_imod5_data( + imod5_cap: GridDataDict, +) -> GridDataArray: + """ + Get rootzone depth from imod5 capillary zone data. Also does a unit + conversion: iMOD5 specifies rootzone thickness in centimeters, whereas + MetaSWAP requires rootzone depth in meters. + """ + rootzone_thickness = imod5_cap["rootzone_thickness"] * 0.01 + # rootzone depth is equal for both svats. + return concat_imod5(rootzone_thickness, rootzone_thickness) + + +def is_msw_active_cell( + target_dis: StructuredDiscretization, + imod5_cap: GridDataDict, + msw_area: GridDataArray, +) -> tuple[GridDataArray, GridDataArray]: + """ + Return grid of cells that are active in the coupled computation, based on + following criteria: + + - Active in top layer MODFLOW6 + - Active in boundary array in CAP package + - MetaSWAP area > 0 + + Returns + ------- + active: xr.DataArray + Active cells in any of the subunits + subunit_active: xr.DataArray + Cells active per subunit + """ + mf6_top_active = target_dis["idomain"].isel(layer=0, drop=True) + subunit_active = ( + (imod5_cap["boundary"] == 1) & (msw_area > 0) & (mf6_top_active >= 1) + ) + active = subunit_active.any(dim="subunit") + return active, subunit_active From c52059998912c82b337dad1e853bcda8b4ce0557 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Mon, 2 Dec 2024 11:19:52 +0100 Subject: [PATCH 03/10] Cache cell area utility function for reuse --- imod/msw/utilities/imod5_converter.py | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/imod/msw/utilities/imod5_converter.py b/imod/msw/utilities/imod5_converter.py index 70f88c982..891fc7c00 100644 --- a/imod/msw/utilities/imod5_converter.py +++ b/imod/msw/utilities/imod5_converter.py @@ -1,3 +1,5 @@ +from functools import lru_cache + from imod.logging import LogLevel, logger from imod.mf6 import StructuredDiscretization from imod.msw.utilities.common import concat_imod5 @@ -9,10 +11,20 @@ def get_cell_area_from_imod5_data( imod5_cap: GridDataDict, ) -> GridDataArray: - # area's per type of svats - mf6_area = get_cell_area(imod5_cap["wetted_area"]) + # Unpack grids and call into private function, so that only 2 grids have to + # be cached. wetted_area = imod5_cap["wetted_area"] urban_area = imod5_cap["urban_area"] + return _get_cell_area_from_imod5_data(wetted_area, urban_area) + + +@lru_cache(maxsize=2) +def _get_cell_area_from_imod5_data( + wetted_area: GridDataArray, urban_area: GridDataArray +) -> GridDataArray: + # area's per type of svats + mf6_area = get_cell_area(wetted_area) + rural_area = mf6_area - (wetted_area + urban_area) if (wetted_area > mf6_area).any(): logger.log( From 942503fa8a3fb4d653bf16c857218c5616767ff1 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Mon, 2 Dec 2024 11:20:48 +0100 Subject: [PATCH 04/10] Add Recharge.from_imod5_cap_data classmethod --- imod/mf6/rch.py | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/imod/mf6/rch.py b/imod/mf6/rch.py index b202c57ac..deb11e565 100644 --- a/imod/mf6/rch.py +++ b/imod/mf6/rch.py @@ -1,7 +1,8 @@ from copy import deepcopy -from typing import Optional +from typing import Optional, cast import numpy as np +import xarray as xr from imod.logging import init_log_decorator from imod.mf6.boundary_condition import BoundaryCondition @@ -11,6 +12,10 @@ from imod.mf6.utilities.imod5_converter import convert_unit_rch_rate from imod.mf6.utilities.regrid import RegridderWeightsCache, _regrid_package_data from imod.mf6.validation import BOUNDARY_DIMS_SCHEMA, CONC_DIMS_SCHEMA +from imod.msw.utilities.imod5_converter import ( + get_cell_area_from_imod5_data, + is_msw_active_cell, +) from imod.prepare.topsystem.allocation import ALLOCATION_OPTION, allocate_rch_cells from imod.schemata import ( AllInsideNoDataSchema, @@ -23,7 +28,7 @@ IndexesSchema, OtherCoordsSchema, ) -from imod.typing import GridDataArray +from imod.typing import Imod5DataDict, GridDataDict, GridDataArray from imod.typing.grid import ( enforce_dim_order, is_planar_grid, @@ -227,24 +232,27 @@ def from_imod5_data( rch_rate = enforce_dim_order(rch_rate) new_package_data["rate"] = rch_rate - return Recharge(**new_package_data, validate=True, fixed_cell=False) + return cls(**new_package_data, validate=True, fixed_cell=False) @classmethod def from_imod5_cap_data( cls, - imod5_data: dict[str, dict[str, GridDataArray]], + imod5_data: Imod5DataDict, target_dis: StructuredDiscretization, ) -> "Recharge": """ Construct an rch-package from iMOD5 data in the CAP package, loaded with the :func:`imod.formats.prj.open_projectfile_data` function. """ - cap_data = imod5_data["cap"] - new_idomain = target_dis.dataset["idomain"] + cap_data = cast(GridDataDict, imod5_data["cap"]) + + msw_area = get_cell_area_from_imod5_data(cap_data) + active, _ = is_msw_active_cell(target_dis, cap_data, msw_area) - # Find which cells are active, basically same logic as in - # msw.GridData.from_imod5_data. Maybe cache that somewhere? + data = {} + data["rate"] = xr.DataArray(0.0).where(active) + return cls(**data, validate=True, fixed_cell=False) @classmethod def get_regrid_methods(cls) -> RechargeRegridMethod: From dddc75e5ae50de72a52ed61fa6a55aed257b73f7 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Mon, 2 Dec 2024 11:45:22 +0100 Subject: [PATCH 05/10] Revert "Cache cell area utility function for reuse" This reverts commit c52059998912c82b337dad1e853bcda8b4ce0557. --- imod/msw/utilities/imod5_converter.py | 16 ++-------------- 1 file changed, 2 insertions(+), 14 deletions(-) diff --git a/imod/msw/utilities/imod5_converter.py b/imod/msw/utilities/imod5_converter.py index 891fc7c00..70f88c982 100644 --- a/imod/msw/utilities/imod5_converter.py +++ b/imod/msw/utilities/imod5_converter.py @@ -1,5 +1,3 @@ -from functools import lru_cache - from imod.logging import LogLevel, logger from imod.mf6 import StructuredDiscretization from imod.msw.utilities.common import concat_imod5 @@ -11,20 +9,10 @@ def get_cell_area_from_imod5_data( imod5_cap: GridDataDict, ) -> GridDataArray: - # Unpack grids and call into private function, so that only 2 grids have to - # be cached. + # area's per type of svats + mf6_area = get_cell_area(imod5_cap["wetted_area"]) wetted_area = imod5_cap["wetted_area"] urban_area = imod5_cap["urban_area"] - return _get_cell_area_from_imod5_data(wetted_area, urban_area) - - -@lru_cache(maxsize=2) -def _get_cell_area_from_imod5_data( - wetted_area: GridDataArray, urban_area: GridDataArray -) -> GridDataArray: - # area's per type of svats - mf6_area = get_cell_area(wetted_area) - rural_area = mf6_area - (wetted_area + urban_area) if (wetted_area > mf6_area).any(): logger.log( From b5e1c895e9462f189b66f542fe33f4107ccf55f9 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Mon, 2 Dec 2024 13:09:00 +0100 Subject: [PATCH 06/10] Recharge first layer to zero --- imod/mf6/rch.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/imod/mf6/rch.py b/imod/mf6/rch.py index deb11e565..1a070c098 100644 --- a/imod/mf6/rch.py +++ b/imod/mf6/rch.py @@ -28,7 +28,7 @@ IndexesSchema, OtherCoordsSchema, ) -from imod.typing import Imod5DataDict, GridDataDict, GridDataArray +from imod.typing import GridDataArray, GridDataDict, Imod5DataDict from imod.typing.grid import ( enforce_dim_order, is_planar_grid, @@ -250,7 +250,11 @@ def from_imod5_cap_data( active, _ = is_msw_active_cell(target_dis, cap_data, msw_area) data = {} - data["rate"] = xr.DataArray(0.0).where(active) + layer_da = xr.full_like( + target_dis.dataset.coords["layer"], np.nan, dtype=np.float64 + ) + layer_da.loc[{"layer": 1}] = 0.0 + data["rate"] = layer_da.where(active) return cls(**data, validate=True, fixed_cell=False) From fed4b063526640b5de3ba753793f71002cd52831 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Mon, 2 Dec 2024 13:09:27 +0100 Subject: [PATCH 07/10] Add test for from_imod5_cap_data --- imod/tests/test_mf6/test_mf6_rch.py | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/imod/tests/test_mf6/test_mf6_rch.py b/imod/tests/test_mf6/test_mf6_rch.py index d8da6a66b..fac9f2a95 100644 --- a/imod/tests/test_mf6/test_mf6_rch.py +++ b/imod/tests/test_mf6/test_mf6_rch.py @@ -455,3 +455,30 @@ def test_non_planar_rch_from_imod5_transient(imod5_dataset, tmp_path): ) assert rendered_rch.count("begin period") == 3 assert "maxbound 33856" in rendered_rch + + +@pytest.mark.usefixtures("imod5_dataset") +def test_from_imod5_cap_data(imod5_dataset): + # Arrange + data = deepcopy(imod5_dataset[0]) + target_discretization = StructuredDiscretization.from_imod5_data(data) + data["cap"] = {} + msw_bound = data["bnd"]["ibound"].isel(layer=0, drop=True) + data["cap"]["boundary"] = msw_bound + data["cap"]["wetted_area"] = xr.ones_like(msw_bound) * 100 + data["cap"]["urban_area"] = xr.ones_like(msw_bound) * 200 + # Set to total cellsize, cell needs to be deactivated. + data["cap"]["wetted_area"][100, 100] = 625.0 + # Set to zero, cell needs to be deactivated. + data["cap"]["urban_area"][100, 100] = 0.0 + # Act + rch = imod.mf6.Recharge.from_imod5_cap_data(data, target_discretization) + rate = rch.dataset["rate"] + # Assert + np.testing.assert_array_equal(np.unique(rate), np.array([0.0, np.nan])) + # Boundaries inactive in MetaSWAP + assert np.isnan(rate[0, :, 0]).all() + assert np.isnan(rate[0, :, -1]).all() + assert np.isnan(rate[0, 0, :]).all() + assert np.isnan(rate[0, -1, :]).all() + assert np.isnan(rate[:, 100, 100]).all() From a13bc316ec5063d2585feb26084525c1a4dadb35 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Mon, 2 Dec 2024 13:23:56 +0100 Subject: [PATCH 08/10] Update changelog --- docs/api/changelog.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst index dbba23f5e..330bf8144 100644 --- a/docs/api/changelog.rst +++ b/docs/api/changelog.rst @@ -14,6 +14,8 @@ 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.Recharge.from_imod5_cap_data` to construct a recharge package + for coupling a MODFLOW6 model to MetaSWAP. Fixed ~~~~~ From ef7aee460dd3551653f22bcf96190bdffd8c1a87 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Mon, 2 Dec 2024 13:32:31 +0100 Subject: [PATCH 09/10] Extend docstring --- imod/mf6/rch.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/imod/mf6/rch.py b/imod/mf6/rch.py index 1a070c098..392cc9d56 100644 --- a/imod/mf6/rch.py +++ b/imod/mf6/rch.py @@ -242,7 +242,9 @@ def from_imod5_cap_data( ) -> "Recharge": """ Construct an rch-package from iMOD5 data in the CAP package, loaded with - the :func:`imod.formats.prj.open_projectfile_data` function. + the :func:`imod.formats.prj.open_projectfile_data` function. Package is + used to couple MODFLOW6 to MetaSWAP models. Active cells will have a + recharge rate of 0.0. """ cap_data = cast(GridDataDict, imod5_data["cap"]) From d58e47e0babf56460bc2f99b480145bdce41a72b Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Mon, 2 Dec 2024 13:44:55 +0100 Subject: [PATCH 10/10] Add Recharge.from_imod5_cap_data if 'cap' in imod5_data --- imod/mf6/model_gwf.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/imod/mf6/model_gwf.py b/imod/mf6/model_gwf.py index 3197cfa91..9ae5ece06 100644 --- a/imod/mf6/model_gwf.py +++ b/imod/mf6/model_gwf.py @@ -330,6 +330,9 @@ def from_imod5_data( wel_key, imod5_data, times ) + if "cap" in imod5_keys: + 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"]