From 69187862bbe22e8e219296b85ac5f1751de869f2 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Thu, 14 Nov 2024 17:18:48 +0100 Subject: [PATCH 01/13] Issue #1285 grid data from imod5 (#1286) Fixes #1285 # Description - This PR implements ``GridData.from_imod5_data``, which has the following requirements: - Area: Split area in rural & urban area - Active: Split where active/inactive - Landuse: Set urban_landuse to 18 - Rootzone thickness: convert from centimeter to meter - Refactors the ``test_grid_data.py`` to separate cases, reduces code duplication - Add unittest for ``GridData.from_imod5_data`` It has a arguments for regridding, but these don't do anything yet. I wasn't sure if I should keep this to already have a future-proof signature, or not add these, as the regridder arguments only confuse people. I don't think regridding the MetaSWAP iMOD5 grids (like we do with MODFLOW6 data) is a part of iMOD Python now, as we are first focusing on the LHM, for which this is not a hard requirement. # Checklist - [x] Links to correct issue - [ ] Update changelog, if changes affect users - [x] PR title starts with ``Issue #nr``, e.g. ``Issue #737`` - [x] Unit tests were added - [ ] **If feature added**: Added/extended example --- imod/msw/grid_data.py | 101 ++++++ imod/msw/pkgbase.py | 3 + imod/tests/test_msw/test_grid_data.py | 504 +++++++++----------------- imod/typing/__init__.py | 1 + 4 files changed, 267 insertions(+), 342 deletions(-) diff --git a/imod/msw/grid_data.py b/imod/msw/grid_data.py index 25003716b..27c08829c 100644 --- a/imod/msw/grid_data.py +++ b/imod/msw/grid_data.py @@ -1,13 +1,80 @@ +from typing import Optional + 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 ( + RegridMethodType, +) +from imod.mf6.utilities.regrid import RegridderWeightsCache from imod.msw.fixed_format import VariableMetaData from imod.msw.pkgbase import MetaSwapPackage from imod.msw.regrid.regrid_schemes import GridDataRegridMethod +from imod.typing import GridDataArray, GridDataDict +from imod.typing.grid import concat, ones_like from imod.util.spatial import get_cell_area, spatial_reference +def _concat_subunits(arg1: GridDataArray, arg2: GridDataArray): + return concat([arg1, arg2], dim="subunit").assign_coords(subunit=[0, 1]) + + +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_subunits(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_subunits(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_subunits(rootzone_thickness, rootzone_thickness) + + class GridData(MetaSwapPackage, IRegridPackage): """ This contains the grid data of MetaSWAP. @@ -111,3 +178,37 @@ def _pkgcheck(self): raise ValueError( "Provided area grid with total areas larger than cell area" ) + + @classmethod + def from_imod5_data( + cls, + imod5_data: dict[str, GridDataDict], + target_dis: StructuredDiscretization, + regridder_types: Optional[RegridMethodType] = None, + regrid_cache: RegridderWeightsCache = RegridderWeightsCache(), + ) -> "GridData": + # Get iMOD5 capillary zone data + imod5_cap = imod5_data["cap"] + + data = {} + data["area"] = get_cell_area_from_imod5_data(imod5_cap) + data["landuse"] = get_landuse_from_imod5_data(imod5_cap) + data["rootzone_depth"] = get_rootzone_depth_from_imod5_data(imod5_cap) + 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") + data_active = { + key: ( + griddata.where(subunit_active) + if key in cls._with_subunit + else griddata.where(active) + ) + for key, griddata in data.items() + } + data_active["active"] = active + return cls(**data_active) diff --git a/imod/msw/pkgbase.py b/imod/msw/pkgbase.py index 5c99c1053..08b306e6d 100644 --- a/imod/msw/pkgbase.py +++ b/imod/msw/pkgbase.py @@ -240,3 +240,6 @@ def regrid_like( def get_regrid_methods(self) -> RegridMethodType: return deepcopy(self._regrid_method) + + def from_imod5_data(self, *args, **kwargs): + raise NotImplementedError("Method not implemented for this package.") diff --git a/imod/tests/test_msw/test_grid_data.py b/imod/tests/test_msw/test_grid_data.py index a3b62a64b..0d0fbf0c3 100644 --- a/imod/tests/test_msw/test_grid_data.py +++ b/imod/tests/test_msw/test_grid_data.py @@ -1,4 +1,5 @@ import tempfile +from copy import deepcopy from pathlib import Path import numpy as np @@ -8,7 +9,9 @@ from hypothesis.strategies import floats from numpy import nan from numpy.testing import assert_almost_equal, assert_equal +from pytest_cases import case, parametrize_with_cases +from imod.mf6.dis import StructuredDiscretization from imod.mf6.utilities.regrid import RegridderWeightsCache from imod.msw import GridData from imod.msw.fixed_format import format_fixed_width @@ -103,147 +106,107 @@ def test_write( assert_equal(results["soil_physical_unit"][0], int(soil_physical_unit)) -def test_generate_index_array(): +@pytest.fixture(scope="function") +def coords_planar() -> dict: x = [1.0, 2.0, 3.0] - y = [1.0, 2.0, 3.0] - subunit = [0, 1] + y = [3.0, 2.0, 1.0] dx = 1.0 dy = 1.0 + return {"y": y, "x": x, "dx": dx, "dy": dy} + + +@pytest.fixture(scope="function") +def coords_one_subunit(coords_planar: dict) -> dict: + coords_subunit = deepcopy(coords_planar) + coords_subunit["subunit"] = [0] + return coords_subunit + + +@pytest.fixture(scope="function") +def coords_two_subunit(coords_planar: dict) -> dict: + coords_subunit = deepcopy(coords_planar) + coords_subunit["subunit"] = [0, 1] + return coords_subunit + + +@case(tags="one_subunit") +def case_grid_data_one_subunits( + coords_one_subunit: dict, coords_planar: dict +) -> dict[str, xr.DataArray]: + data = {} # fmt: off - area = xr.DataArray( + data["area"] = xr.DataArray( np.array( [ [[0.5, 0.5, 0.5], - [nan, nan, nan], + [0.7, 0.7, 0.7], [1.0, 1.0, 1.0]], - [[0.5, 0.5, 0.5], - [1.0, 1.0, 1.0], - [nan, nan, nan]], ] ), dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} + coords=coords_one_subunit ) - landuse = xr.DataArray( + data["landuse"] = xr.DataArray( np.array( [ [[1.0, 1.0, 1.0], - [nan, nan, nan], + [1.0, 1.0, 1.0], [1.0, 1.0, 1.0]], - - [[2.0, 2.0, 2.0], - [2.0, 2.0, 2.0], - [nan, nan, nan]], ] ), dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} + coords=coords_one_subunit ) - rootzone_depth = xr.DataArray( + data["rootzone_depth"] = xr.DataArray( np.array( [ - [[1.0, 1.0, 1.0], - [nan, nan, nan], - [1.0, 1.0, 1.0]], - [[1.0, 1.0, 1.0], [1.0, 1.0, 1.0], - [nan, nan, nan]], + [1.0, 1.0, 1.0]], ] ), dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} + coords=coords_one_subunit ) - surface_elevation = xr.DataArray( + data["surface_elevation"] = xr.DataArray( np.array( [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]), dims=("y", "x"), - coords={"y": y, "x": x, "dx": dx, "dy": dy} + coords=coords_planar ) - soil_physical_unit = xr.DataArray( + data["soil_physical_unit"] = xr.DataArray( np.array( [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]), dims=("y", "x"), - coords={"y": y, "x": x, "dx": dx, "dy": dy} + coords=coords_planar ) - active = xr.DataArray( + data["active"] = xr.DataArray( np.array( [[False, True, False], [False, True, False], [False, True, False]]), dims=("y", "x"), - coords={"y": y, "x": x} - ) - # fmt: on - grid_data = GridData( - area, - landuse, - rootzone_depth, - surface_elevation, - soil_physical_unit, - active, - ) - - index, svat = grid_data.generate_index_array() - - index_expected = [ - False, - True, - False, - False, - False, - False, - False, - True, - False, - False, - True, - False, - False, - True, - False, - False, - False, - False, - ] - - # fmt: off - svat_expected = xr.DataArray( - np.array( - [ - [[0, 1, 0], - [0, 0, 0], - [0, 2, 0]], - - [[0, 3, 0], - [0, 4, 0], - [0, 0, 0]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} + coords=coords_planar ) # fmt: on - assert_equal(index, np.array(index_expected)) - assert_equal(svat.values, svat_expected.values) + return data -def simple_model(): - x = [1.0, 2.0, 3.0] - y = [3.0, 2.0, 1.0] - subunit = [0, 1] - dx = 1.0 - dy = 1.0 +@case(tags="two_subunit") +def case_grid_data_two_subunits( + coords_two_subunit: dict, coords_planar: dict +) -> dict[str, xr.DataArray]: + data = {} # fmt: off - area = xr.DataArray( + data["area"] = xr.DataArray( np.array( [ [[0.5, 0.5, 0.5], @@ -256,9 +219,9 @@ def simple_model(): ] ), dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} + coords=coords_two_subunit ) - landuse = xr.DataArray( + data["landuse"] = xr.DataArray( np.array( [ [[1.0, 1.0, 1.0], @@ -271,9 +234,9 @@ def simple_model(): ] ), dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} + coords=coords_two_subunit ) - rootzone_depth = xr.DataArray( + data["rootzone_depth"] = xr.DataArray( np.array( [ [[1.0, 1.0, 1.0], @@ -286,51 +249,92 @@ def simple_model(): ] ), dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} + coords=coords_two_subunit ) - surface_elevation = xr.DataArray( + data["surface_elevation"] = xr.DataArray( np.array( [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]), dims=("y", "x"), - coords={"y": y, "x": x, "dx": dx, "dy": dy} + coords=coords_planar ) - soil_physical_unit = xr.DataArray( + data["soil_physical_unit"] = xr.DataArray( np.array( [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]), dims=("y", "x"), - coords={"y": y, "x": x, "dx": dx, "dy": dy} + coords=coords_planar ) - active = xr.DataArray( + data["active"] = xr.DataArray( np.array( [[False, True, False], [False, True, False], [False, True, False]]), dims=("y", "x"), - coords={"y": y, "x": x} + coords=coords_planar ) # fmt: on + return data - grid_data = GridData( - area, - landuse, - rootzone_depth, - surface_elevation, - soil_physical_unit, - active, - ) - return grid_data +@parametrize_with_cases("grid_data_dict", cases=".", has_tag="two_subunit") +def test_generate_index_array( + grid_data_dict: dict[str, xr.DataArray], coords_two_subunit: dict +): + grid_data = GridData(**grid_data_dict) + + index, svat = grid_data.generate_index_array() + index_expected = [ + False, + True, + False, + False, + False, + False, + False, + True, + False, + False, + True, + False, + False, + True, + False, + False, + False, + False, + ] -def test_simple_model(fixed_format_parser): - grid_data = simple_model() + # fmt: off + svat_expected = xr.DataArray( + np.array( + [ + [[0, 1, 0], + [0, 0, 0], + [0, 2, 0]], + + [[0, 3, 0], + [0, 4, 0], + [0, 0, 0]], + ] + ), + dims=("subunit", "y", "x"), + coords=coords_two_subunit, + ) + # fmt: on + assert_equal(index, np.array(index_expected)) + assert_equal(svat.values, svat_expected.values) + + +@parametrize_with_cases("grid_data_dict", cases=".", has_tag="two_subunit") +def test_simple_model(fixed_format_parser, grid_data_dict: dict[str, xr.DataArray]): + grid_data = GridData(**grid_data_dict) index, svat = grid_data.generate_index_array() @@ -350,8 +354,13 @@ def test_simple_model(fixed_format_parser): assert_almost_equal(results["rootzone_depth"], np.array([1.0, 1.0, 1.0, 1.0])) -def test_simple_model_regrid(simple_2d_grid_with_subunits): - grid_data = simple_model() +# Only use two subunit case, as this one sums to total area whereas the one +# subunit case doesn't sum to 1 across subunit. +@parametrize_with_cases("grid_data_dict", cases=".", has_tag="two_subunit") +def test_simple_model_regrid( + simple_2d_grid_with_subunits, grid_data_dict: dict[str, xr.DataArray] +): + grid_data = GridData(**grid_data_dict) new_grid = simple_2d_grid_with_subunits regrid_context = RegridderWeightsCache() @@ -363,84 +372,11 @@ def test_simple_model_regrid(simple_2d_grid_with_subunits): assert np.sum(regridded_area.values) == regridded_total_area -def test_simple_model_1_subunit(fixed_format_parser): - x = [1.0, 2.0, 3.0] - y = [1.0, 2.0, 3.0] - subunit = [0] - dx = 1.0 - dy = 1.0 - # fmt: off - area = xr.DataArray( - np.array( - [ - [[0.5, 0.5, 0.5], - [0.7, 0.7, 0.7], - [1.0, 1.0, 1.0]], - - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} - ) - landuse = xr.DataArray( - np.array( - [ - [[1.0, 1.0, 1.0], - [1.0, 1.0, 1.0], - [1.0, 1.0, 1.0]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} - ) - rootzone_depth = xr.DataArray( - np.array( - [ - [[1.0, 1.0, 1.0], - [1.0, 1.0, 1.0], - [1.0, 1.0, 1.0]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} - ) - - surface_elevation = xr.DataArray( - np.array( - [[1.0, 2.0, 3.0], - [4.0, 5.0, 6.0], - [7.0, 8.0, 9.0]]), - dims=("y", "x"), - coords={"y": y, "x": x, "dx": dx, "dy": dy} - ) - - soil_physical_unit = xr.DataArray( - np.array( - [[1.0, 2.0, 3.0], - [4.0, 5.0, 6.0], - [7.0, 8.0, 9.0]]), - dims=("y", "x"), - coords={"y": y, "x": x, "dx": dx, "dy": dy} - ) - - active = xr.DataArray( - np.array( - [[False, True, False], - [False, True, False], - [False, True, False]]), - dims=("y", "x"), - coords={"y": y, "x": x} - ) - # fmt: on - - grid_data = GridData( - area, - landuse, - rootzone_depth, - surface_elevation, - soil_physical_unit, - active, - ) +@parametrize_with_cases("grid_data_dict", cases=".", has_tag="one_subunit") +def test_simple_model_1_subunit( + fixed_format_parser, grid_data_dict: dict[str, xr.DataArray] +): + grid_data = GridData(**grid_data_dict) index, svat = grid_data.generate_index_array() @@ -460,17 +396,15 @@ def test_simple_model_1_subunit(fixed_format_parser): assert_almost_equal(results["rootzone_depth"], np.array([1.0, 1.0, 1.0])) -def test_area_grid_exceeds_cell_area(): +@parametrize_with_cases("grid_data_dict", cases=".", has_tag="two_subunit") +def test_area_grid_exceeds_cell_area( + grid_data_dict: dict[str, xr.DataArray], coords_two_subunit: dict +): """ Test where provided area grid exceeds total cell area, should throw error. """ - x = [1.0, 3.0, 5.0] - y = [1.0, 2.0, 3.0] - subunit = [0, 1] - dx = 1.0 - dy = 1.0 # fmt: off - area = xr.DataArray( + grid_data_dict["area"] = xr.DataArray( np.array( [ [[0.5, 0.5, 0.5], @@ -483,167 +417,53 @@ def test_area_grid_exceeds_cell_area(): ] ), dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} - ) - landuse = xr.DataArray( - np.array( - [ - [[1.0, 1.0, 1.0], - [nan, nan, nan], - [1.0, 1.0, 1.0]], - - [[2.0, 2.0, 2.0], - [2.0, 2.0, 2.0], - [nan, nan, nan]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} - ) - rootzone_depth = xr.DataArray( - np.array( - [ - [[1.0, 1.0, 1.0], - [nan, nan, nan], - [1.0, 1.0, 1.0]], - - [[1.0, 1.0, 1.0], - [1.0, 1.0, 1.0], - [nan, nan, nan]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} - ) - - surface_elevation = xr.DataArray( - np.array( - [[1.0, 2.0, 3.0], - [4.0, 5.0, 6.0], - [7.0, 8.0, 9.0]]), - dims=("y", "x"), - coords={"y": y, "x": x, "dx": dx, "dy": dy} - ) - - soil_physical_unit = xr.DataArray( - np.array( - [[1.0, 2.0, 3.0], - [4.0, 5.0, 6.0], - [7.0, 8.0, 9.0]]), - dims=("y", "x"), - coords={"y": y, "x": x, "dx": dx, "dy": dy} - ) - - active = xr.DataArray( - np.array( - [[False, True, False], - [False, True, False], - [False, True, False]]), - dims=("y", "x"), - coords={"y": y, "x": x} + coords=coords_two_subunit ) # fmt: on with pytest.raises(ValueError): - GridData( - area, - landuse, - rootzone_depth, - surface_elevation, - soil_physical_unit, - active, - ) + GridData(**grid_data_dict) -def test_non_equidistant(): +@parametrize_with_cases("grid_data_dict", cases=".") +def test_non_equidistant(grid_data_dict: dict[str, xr.DataArray]): """ Test where provided grid is non-equidistant, should throw error. """ - x = [1.0, 2.0, 5.0] - y = [1.0, 2.0, 5.0] - subunit = [0, 1] dx = [1.0, 1.0, 5.0] dy = [1.0, 1.0, 5.0] - # fmt: off - area = xr.DataArray( - np.array( - [ - [[0.5, 0.5, 0.5], - [nan, nan, nan], - [1.0, 1.0, 1.0]], + for key, value in grid_data_dict.items(): + grid_data_dict[key] = value.assign_coords(dx=("x", dx), dy=("y", dy)) - [[0.5, 0.5, 0.5], - [1.0, 1.0, 1.0], - [nan, nan, nan]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": ("x", dx), "dy": ("y", dy)} - ) - landuse = xr.DataArray( - np.array( - [ - [[1.0, 1.0, 1.0], - [nan, nan, nan], - [1.0, 1.0, 1.0]], + with pytest.raises(ValueError): + GridData(**grid_data_dict) - [[2.0, 2.0, 2.0], - [2.0, 2.0, 2.0], - [nan, nan, nan]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": ("x", dx), "dy": ("y", dy)} - ) - rootzone_depth = xr.DataArray( - np.array( - [ - [[1.0, 1.0, 1.0], - [nan, nan, nan], - [1.0, 1.0, 1.0]], - [[1.0, 1.0, 1.0], - [1.0, 1.0, 1.0], - [nan, nan, nan]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": ("x", dx), "dy": ("y", dy)} - ) +@parametrize_with_cases("grid_data_dict", cases=".") +def test_from_imod5_data(grid_data_dict: dict[str, xr.DataArray]): + cap_data = {} + cap_data["wetted_area"] = 1 - grid_data_dict["area"].sum(dim="subunit") + like = grid_data_dict["area"].sel(subunit=0, drop=True) + top = xr.zeros_like(like, dtype=float) + cap_data["boundary"] = xr.ones_like(like, dtype=int) + cap_data["urban_area"] = xr.zeros_like(like, dtype=float) + 0.1 + cap_data["landuse"] = xr.ones_like(like, dtype=int) + cap_data["rootzone_thickness"] = xr.ones_like(like, dtype=int) + cap_data["surface_elevation"] = top + cap_data["soil_physical_unit"] = xr.ones_like(like, dtype=int) + cap_data["active"] = xr.ones_like(like, dtype=bool) - surface_elevation = xr.DataArray( - np.array( - [[1.0, 2.0, 3.0], - [4.0, 5.0, 6.0], - [7.0, 8.0, 9.0]]), - dims=("y", "x"), - coords={"y": y, "x": x, "dx": ("x", dx), "dy": ("y", dy)} - ) + imod5_data = {"cap": cap_data} - soil_physical_unit = xr.DataArray( - np.array( - [[1.0, 2.0, 3.0], - [4.0, 5.0, 6.0], - [7.0, 8.0, 9.0]]), - dims=("y", "x"), - coords={"y": y, "x": x, "dx": ("x", dx), "dy": ("y", dy)} - ) + layer = xr.DataArray([1, 1], coords={"layer": [1, 2]}, dims=("layer",)) + idomain = layer * xr.ones_like(like, dtype=int) - active = xr.DataArray( - np.array( - [[False, True, False], - [False, True, False], - [False, True, False]]), - dims=("y", "x"), - coords={"y": y, "x": x} + # Dis only needed for idomain + dis = StructuredDiscretization(top, top - 0.1, idomain, validate=False) + + griddata = GridData.from_imod5_data(imod5_data, target_dis=dis) + expected_rootzone_depth = cap_data["rootzone_thickness"] * 0.01 + xr.testing.assert_allclose( + expected_rootzone_depth, griddata["rootzone_depth"].sel(subunit=0, drop=True) ) - # fmt: on - with pytest.raises(ValueError): - GridData( - area, - landuse, - rootzone_depth, - surface_elevation, - soil_physical_unit, - active, - ) + assert (griddata["landuse"].sel(subunit=1, drop=True) == 18).all() diff --git a/imod/typing/__init__.py b/imod/typing/__init__.py index 0d17b1127..760346894 100644 --- a/imod/typing/__init__.py +++ b/imod/typing/__init__.py @@ -11,6 +11,7 @@ GridDataArray: TypeAlias = Union[xr.DataArray, xu.UgridDataArray] GridDataset: TypeAlias = Union[xr.Dataset, xu.UgridDataset] +GridDataDict: TypeAlias = dict[str, GridDataArray] ScalarAsDataArray: TypeAlias = Union[xr.DataArray, xu.UgridDataArray] ScalarAsDataset: TypeAlias = Union[xr.Dataset, xu.UgridDataset] UnstructuredData: TypeAlias = Union[xu.UgridDataset, xu.UgridDataArray] From 26918e4ab94c32e5e721fa8e51e21727323b5a8e Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Fri, 15 Nov 2024 10:59:10 +0100 Subject: [PATCH 02/13] Issue #1288 Ponding.from_imod5_data (#1289) Fixes #1288 # Description * Adds Ponding.from_imod5_data * Adds msw.utilities namespace to put useful utility functions in * Rename _concat_subunits to concat_imod5, as it is a very imod5 specific function. --- imod/msw/grid_data.py | 13 +++----- imod/msw/ponding.py | 16 +++++++++- imod/msw/utilities/common.py | 6 ++++ imod/tests/test_msw/test_ponding.py | 48 +++++++++++++++++++++++------ 4 files changed, 64 insertions(+), 19 deletions(-) create mode 100644 imod/msw/utilities/common.py diff --git a/imod/msw/grid_data.py b/imod/msw/grid_data.py index 27c08829c..3e08a69b0 100644 --- a/imod/msw/grid_data.py +++ b/imod/msw/grid_data.py @@ -13,15 +13,12 @@ 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 concat, ones_like +from imod.typing.grid import ones_like from imod.util.spatial import get_cell_area, spatial_reference -def _concat_subunits(arg1: GridDataArray, arg2: GridDataArray): - return concat([arg1, arg2], dim="subunit").assign_coords(subunit=[0, 1]) - - def get_cell_area_from_imod5_data( imod5_cap: GridDataDict, ) -> GridDataArray: @@ -45,7 +42,7 @@ def get_cell_area_from_imod5_data( ) urban_area = urban_area.where(rural_area > 0.0, other=0.0) rural_area = mf6_area - (wetted_area + urban_area) - return _concat_subunits(rural_area, urban_area) + return concat_imod5(rural_area, urban_area) def get_landuse_from_imod5_data( @@ -59,7 +56,7 @@ def get_landuse_from_imod5_data( rural_landuse = imod5_cap["landuse"] # Urban landuse = 18 urban_landuse = ones_like(rural_landuse) * 18 - return _concat_subunits(rural_landuse, urban_landuse) + return concat_imod5(rural_landuse, urban_landuse) def get_rootzone_depth_from_imod5_data( @@ -72,7 +69,7 @@ def get_rootzone_depth_from_imod5_data( """ rootzone_thickness = imod5_cap["rootzone_thickness"] * 0.01 # rootzone depth is equal for both svats. - return _concat_subunits(rootzone_thickness, rootzone_thickness) + return concat_imod5(rootzone_thickness, rootzone_thickness) class GridData(MetaSwapPackage, IRegridPackage): diff --git a/imod/msw/ponding.py b/imod/msw/ponding.py index afc1629af..89ab6af44 100644 --- a/imod/msw/ponding.py +++ b/imod/msw/ponding.py @@ -7,7 +7,8 @@ from imod.msw.fixed_format import VariableMetaData from imod.msw.pkgbase import DataDictType, MetaSwapPackage from imod.msw.regrid.regrid_schemes import PondingRegridMethod -from imod.typing import IntArray +from imod.msw.utilities.common import concat_imod5 +from imod.typing import GridDataDict, IntArray class Ponding(MetaSwapPackage, IRegridPackage): @@ -70,3 +71,16 @@ def _render(self, file: TextIO, index: IntArray, svat: xr.DataArray, *args: Any) self._check_range(dataframe) return self.write_dataframe_fixed_width(file, dataframe) + + @classmethod + def from_imod5_data(cls, imod5_data: dict[str, GridDataDict]) -> "Ponding": + """ + Concatenate ponding depths along subunits + """ + cap_data = imod5_data["cap"] + data = {} + for key in cls._with_subunit: + data_ls = [cap_data[f"{landuse}_{key}"] for landuse in ["rural", "urban"]] + data[key] = concat_imod5(*data_ls) + + return cls(**data) diff --git a/imod/msw/utilities/common.py b/imod/msw/utilities/common.py new file mode 100644 index 000000000..836902662 --- /dev/null +++ b/imod/msw/utilities/common.py @@ -0,0 +1,6 @@ +from imod.typing import GridDataArray +from imod.typing.grid import concat + + +def concat_imod5(arg1: GridDataArray, arg2: GridDataArray) -> GridDataArray: + return concat([arg1, arg2], dim="subunit").assign_coords(subunit=[0, 1]) diff --git a/imod/tests/test_msw/test_ponding.py b/imod/tests/test_msw/test_ponding.py index 626cd6412..e7bb47e06 100644 --- a/imod/tests/test_msw/test_ponding.py +++ b/imod/tests/test_msw/test_ponding.py @@ -10,11 +10,12 @@ RegridderWeightsCache, ) from imod.msw import Ponding +from imod.typing import GridDataArray, GridDataDict -def setup_ponding(): +def setup_ponding() -> tuple[GridDataDict, np.ndarray, GridDataArray]: x = [1.0, 2.0, 3.0] - y = [1.0, 2.0, 3.0] + y = [3.0, 2.0, 1.0] subunit = [0, 1] dx = 1.0 dy = 1.0 @@ -69,16 +70,17 @@ def setup_ponding(): # fmt: on index = (svat != 0).values.ravel() - ponding = Ponding( - ponding_depth=ponding_depth, - runoff_resistance=runoff_resistance, - runon_resistance=runoff_resistance, - ) - return ponding, index, svat + data_ponding = { + "ponding_depth": ponding_depth, + "runoff_resistance": runoff_resistance, + "runon_resistance": runoff_resistance, + } + return data_ponding, index, svat def test_simple_model(fixed_format_parser): - ponding, index, svat = setup_ponding() + data_ponding, index, svat = setup_ponding() + ponding = Ponding(**data_ponding) with tempfile.TemporaryDirectory() as output_dir: output_dir = Path(output_dir) ponding.write(output_dir, index, svat, None, None) @@ -95,7 +97,8 @@ def test_simple_model(fixed_format_parser): def test_regrid_ponding(simple_2d_grid_with_subunits): - ponding, _, _ = setup_ponding() + data_ponding, _, _ = setup_ponding() + ponding = Ponding(**data_ponding) new_grid = simple_2d_grid_with_subunits regrid_context = RegridderWeightsCache() @@ -104,3 +107,28 @@ def test_regrid_ponding(simple_2d_grid_with_subunits): assert np.all(regridded_ponding.dataset["x"].values == new_grid["x"].values) assert np.all(regridded_ponding.dataset["y"].values == new_grid["y"].values) + + +def test_from_imod5_data(): + data_ponding, _, _ = setup_ponding() + expected_ponding = Ponding(**data_ponding) + + # Create cap data + cap_data = {} + mapping_ls = [ + ("rural_runoff_resistance", "runoff_resistance", 0), + ("rural_runoff_resistance", "runoff_resistance", 0), + ("urban_runoff_resistance", "runoff_resistance", 1), + ("rural_runon_resistance", "runon_resistance", 0), + ("urban_runon_resistance", "runon_resistance", 1), + ("rural_ponding_depth", "ponding_depth", 0), + ("urban_ponding_depth", "ponding_depth", 1), + ] + for cap_key, pkg_key, subunit_nr in mapping_ls: + cap_data[cap_key] = data_ponding[pkg_key].sel(subunit=subunit_nr, drop=True) + + imod5_data = {"cap": cap_data} + + actual_ponding = Ponding.from_imod5_data(imod5_data) + + xr.testing.assert_equal(expected_ponding.dataset, actual_ponding.dataset) From c2f9f525f9e32791112d99caff06fb0754908335 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Fri, 15 Nov 2024 15:03:52 +0100 Subject: [PATCH 03/13] Issue #1290 infiltration from imod5 (#1291) Fixes #1290 # Description Adds from_imod5_data classmethod. Requirements: * Set runon/runoff resistance to downward/upward resistance respectively * Set extra storage coefficient and bottom resistance to 1 * If resistance smaller than 5.0, deactivate resistance by setting it to MetaSWAP nodata value. I also added the ``upward_resistance`` and ``downward_resistance`` to the list with grids which have a subunit coordinate, as a comment in @HendrikKok's script mentioned this was missing, and a required feature. --- docs/api/changelog.rst | 10 ++ imod/msw/infiltration.py | 64 +++++++- imod/tests/fixtures/msw_model_fixture.py | 4 +- imod/tests/test_msw/test_infiltration.py | 184 +++++++++++++++-------- imod/tests/test_msw/test_ponding.py | 1 - 5 files changed, 192 insertions(+), 71 deletions(-) diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst index 59642ac66..68ebc2591 100644 --- a/docs/api/changelog.rst +++ b/docs/api/changelog.rst @@ -7,6 +7,16 @@ The format is based on `Keep a Changelog`_, and this project adheres to `Semantic Versioning`_. +[Unreleased] +------------ + +Changed +~~~~~~~ + +- :class:`imod.msw.Infiltration`'s variables ``upward_resistance`` and + ``downward_resistance`` now require a ``subunit`` coordinate. + + [0.18.0] -------- diff --git a/imod/msw/infiltration.py b/imod/msw/infiltration.py index bf300761f..830365681 100644 --- a/imod/msw/infiltration.py +++ b/imod/msw/infiltration.py @@ -1,9 +1,35 @@ +from textwrap import dedent + import xarray as xr +from imod.logging import LogLevel, logger from imod.mf6.interfaces.iregridpackage import IRegridPackage from imod.msw.fixed_format import VariableMetaData from imod.msw.pkgbase import MetaSwapPackage from imod.msw.regrid.regrid_schemes import InfiltrationRegridMethod +from imod.msw.utilities.common import concat_imod5 +from imod.typing import GridDataDict +from imod.typing.grid import ones_like + + +def deactivate_small_resistances_in_data(data: GridDataDict): + """ + Deactivate cells where resistance smaller than 5 days are set to + -9999.0. + """ + message = dedent("""Detected cells with resistances smaller than 5.0 in {var}, set + to inactive""") + + for var in ["downward_resistance", "upward_resistance"]: + to_deactivate = data[var] < 5.0 + if to_deactivate.any(): + logger.log( + loglevel=LogLevel.WARNING, + message=message.format(var=var), + additional_depth=1, + ) + data[var] = data[var].where(~to_deactivate, -9999.0) + return data class Infiltration(MetaSwapPackage, IRegridPackage): @@ -19,11 +45,11 @@ class Infiltration(MetaSwapPackage, IRegridPackage): a subunit coordinate to describe different land uses. downward_resistance: array of floats (xr.DataArray) Describes the downward resisitance of SVAT units. Set to -9999.0 to make - MetaSWAP ignore this resistance. This array must not have a subunit + MetaSWAP ignore this resistance. This array must have a subunit coordinate. upward_resistance: array of floats (xr.DataArray) Describes the upward resistance of SVAT units. Set to -9999.0 to make - MetaSWAP ignore this resistance. This array must not have a subunit + MetaSWAP ignore this resistance. This array must have a subunit coordinate. bottom_resistance: array of floats (xr.DataArray) Describes the infiltration capacity of SVAT units. Set to -9999.0 to @@ -32,9 +58,6 @@ class Infiltration(MetaSwapPackage, IRegridPackage): extra_storage_coefficient: array of floats (xr.DataArray) Extra storage coefficient of phreatic layer. This array must not have a subunit coordinate. - active: array of bools (xr.DataArray) - Describes whether SVAT units are active or not. This array must not have - a subunit coordinate. """ _file_name = "infi_svat.inp" @@ -47,10 +70,12 @@ class Infiltration(MetaSwapPackage, IRegridPackage): "extra_storage_coefficient": VariableMetaData(8, 0.01, 1.0, float), } - _with_subunit = ("infiltration_capacity",) - _without_subunit = ( + _with_subunit = ( + "infiltration_capacity", "downward_resistance", "upward_resistance", + ) + _without_subunit = ( "bottom_resistance", "extra_storage_coefficient", ) @@ -74,3 +99,28 @@ def __init__( self.dataset["extra_storage_coefficient"] = extra_storage_coefficient self._pkgcheck() + + @classmethod + def from_imod5_data(cls, imod5_data: dict[str, GridDataDict]) -> "Infiltration": + cap_data = imod5_data["cap"] + data = {} + # Use runon resistance as downward resistance, and runoff for downward + # resistance + key_mapping = { + "infiltration_capacity": "infiltration_capacity", + "downward_resistance": "runon_resistance", + "upward_resistance": "runoff_resistance", + } + for var_rename, var_key in key_mapping.items(): + data_ls = [ + cap_data[f"{landuse}_{var_key}"] for landuse in ["rural", "urban"] + ] + data[var_rename] = concat_imod5(*data_ls) + + data = deactivate_small_resistances_in_data(data) + + like = data["downward_resistance"].isel(subunit=0, drop=True) + data["bottom_resistance"] = ones_like(like) + data["extra_storage_coefficient"] = ones_like(like) + + return cls(**data) diff --git a/imod/tests/fixtures/msw_model_fixture.py b/imod/tests/fixtures/msw_model_fixture.py index 3b7c73c80..1fcaf9015 100644 --- a/imod/tests/fixtures/msw_model_fixture.py +++ b/imod/tests/fixtures/msw_model_fixture.py @@ -226,8 +226,8 @@ def make_msw_model(): # %% Infiltration msw_model["infiltration"] = msw.Infiltration( infiltration_capacity=xr.full_like(area, 1.0), - downward_resistance=xr.full_like(msw_grid, -9999.0), - upward_resistance=xr.full_like(msw_grid, -9999.0), + downward_resistance=xr.full_like(area, -9999.0), + upward_resistance=xr.full_like(area, -9999.0), bottom_resistance=xr.full_like(msw_grid, -9999.0), extra_storage_coefficient=xr.full_like(msw_grid, 0.1), ) diff --git a/imod/tests/test_msw/test_infiltration.py b/imod/tests/test_msw/test_infiltration.py index d3cb0a003..b87940a88 100644 --- a/imod/tests/test_msw/test_infiltration.py +++ b/imod/tests/test_msw/test_infiltration.py @@ -1,22 +1,60 @@ import tempfile +from copy import deepcopy from pathlib import Path import numpy as np +import pytest import xarray as xr from hypothesis import given, settings from hypothesis.strategies import floats from numpy import nan from numpy.testing import assert_almost_equal, assert_equal +from pytest_cases import case, parametrize_with_cases from imod.mf6.utilities.regrid import ( RegridderWeightsCache, ) from imod.msw import Infiltration from imod.msw.fixed_format import format_fixed_width +from imod.typing import GridDataDict -def setup_infiltration_package(subunit, y, x, dy, dx): - infiltration_capacity = xr.DataArray( +@pytest.fixture(scope="function") +def coords_planar() -> dict: + x = [1.0, 2.0, 3.0] + y = [3.0, 2.0, 1.0] + dx = 1.0 + dy = 1.0 + return {"y": y, "x": x, "dx": dx, "dy": dy} + + +@pytest.fixture(scope="function") +def coords_subunit(coords_planar: dict) -> dict: + coords_subunit = deepcopy(coords_planar) + coords_subunit["subunit"] = [0, 1] + return coords_subunit + + +@pytest.fixture(scope="function") +def svat_index(coords_subunit: dict) -> tuple[xr.DataArray, np.ndarray]: + svat = xr.DataArray( + np.array( + [ + [[0, 1, 0], [0, 0, 0], [0, 2, 0]], + [[0, 3, 0], [0, 4, 0], [0, 0, 0]], + ] + ), + dims=("subunit", "y", "x"), + coords=coords_subunit, + ) + index = (svat != 0).values.ravel() + return svat, index + + +@pytest.fixture(scope="function") +def setup_infiltration_data(coords_planar, coords_subunit) -> GridDataDict: + data = {} + data["infiltration_capacity"] = xr.DataArray( np.array( [ [[0.5, 0.5, 0.5], [nan, nan, nan], [1.0, 1.0, 1.0]], @@ -24,55 +62,53 @@ def setup_infiltration_package(subunit, y, x, dy, dx): ] ), dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy}, + coords=coords_subunit, ) - - downward_resistance = xr.DataArray( - np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]), - dims=("y", "x"), - coords={"y": y, "x": x, "dx": dx, "dy": dy}, + data["downward_resistance"] = xr.DataArray( + np.array( + [ + [[1.0, 2.0, 3.0], [nan, nan, nan], [7.0, 8.0, 9.0]], + [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [nan, nan, nan]], + ] + ), + dims=("subunit", "y", "x"), + coords=coords_subunit, ) - - upward_resistance = xr.DataArray( - np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]), - dims=("y", "x"), - coords={"y": y, "x": x, "dx": dx, "dy": dy}, + data["upward_resistance"] = xr.DataArray( + np.array( + [ + [[1.0, 2.0, 3.0], [nan, nan, nan], [7.0, 8.0, 9.0]], + [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [nan, nan, nan]], + ] + ), + dims=("subunit", "y", "x"), + coords=coords_subunit, ) - - bottom_resistance = xr.DataArray( + data["bottom_resistance"] = xr.DataArray( np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]), dims=("y", "x"), - coords={"y": y, "x": x, "dx": dx, "dy": dy}, + coords=coords_planar, ) - - extra_storage_coefficient = xr.DataArray( + data["extra_storage_coefficient"] = xr.DataArray( np.array([[0.1, 0.2, 0.3], [0.4, 0.5, 0.6], [0.7, 0.8, 0.9]]), dims=("y", "x"), - coords={"y": y, "x": x, "dx": dx, "dy": dy}, + coords=coords_planar, ) - svat = xr.DataArray( - np.array( - [ - [[0, 1, 0], [0, 0, 0], [0, 2, 0]], - [[0, 3, 0], [0, 4, 0], [0, 0, 0]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy}, - ) - # fmt: on - index = (svat != 0).values.ravel() + return data - infiltration = Infiltration( - infiltration_capacity, - downward_resistance, - upward_resistance, - bottom_resistance, - extra_storage_coefficient, - ) - return infiltration, svat, index +@case(tags="r_low") +def case_low_resistance(setup_infiltration_data: GridDataDict) -> GridDataDict: + return setup_infiltration_data + + +@case(tags="r_high") +def case_high_resistance(setup_infiltration_data: GridDataDict) -> GridDataDict: + data = setup_infiltration_data + data["downward_resistance"] += 10.0 + data["upward_resistance"] += 10.0 + return data @given( @@ -108,8 +144,8 @@ def test_write( ): infiltration = Infiltration( xr.DataArray(infiltration_capacity).expand_dims(subunit=[0]), - xr.DataArray(downward_resistance), - xr.DataArray(upward_resistance), + xr.DataArray(downward_resistance).expand_dims(subunit=[0]), + xr.DataArray(upward_resistance).expand_dims(subunit=[0]), xr.DataArray(bottom_resistance), xr.DataArray(extra_storage_coefficient), ) @@ -175,14 +211,10 @@ def test_write( ) -def test_simple_model(fixed_format_parser): - x = [1.0, 2.0, 3.0] - y = [1.0, 2.0, 3.0] - subunit = [0, 1] - dx = 1.0 - dy = 1.0 - # fmt: off - infiltration, svat, index = setup_infiltration_package(subunit, y, x, dy, dx) +@parametrize_with_cases("infiltration_data", cases=".", has_tag="r_low") +def test_simple_model(fixed_format_parser, svat_index, infiltration_data): + svat, index = svat_index + infiltration = Infiltration(**infiltration_data) with tempfile.TemporaryDirectory() as output_dir: output_dir = Path(output_dir) @@ -204,28 +236,58 @@ def test_simple_model(fixed_format_parser): ) -def test_regrid(): - x = [1.0, 2.0, 3.0] - y = [3.0, 2.0, 1.0] - subunit = [0, 1] - dx = 1.0 - dy = 1.0 - - infiltration, _, _ = setup_infiltration_package(subunit, y, x, dy, dx) +@parametrize_with_cases("infiltration_data", cases=".", has_tag="r_low") +def test_regrid(infiltration_data): + infiltration = Infiltration(**infiltration_data) x = [1.0, 1.5, 2.0, 2.5, 3.0] y = [3.0, 2.5, 2.0, 1.5, 1.0] subunit = [0, 1] dx = 0.5 dy = 0.5 - # fmt: off new_grid = xr.DataArray( dims=("subunit", "y", "x"), - coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} + coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy}, ) - new_grid.values[:,:,:] = 1 + new_grid.values[:, :, :] = 1 regrid_context = RegridderWeightsCache() - regridded = infiltration.regrid_like(new_grid, regrid_context ) + regridded = infiltration.regrid_like(new_grid, regrid_context) assert_almost_equal(regridded.dataset.coords["x"].values, x) assert_almost_equal(regridded.dataset.coords["y"].values, y) + + +@parametrize_with_cases("data_infiltration", cases=".") +def test_from_imod5_data(data_infiltration): + expected_pkg = Infiltration(**data_infiltration) + # Deactivate cells which have a resistance lower than 5.0 + for var in ["upward_resistance", "downward_resistance"]: + da = expected_pkg.dataset[var] + to_deactivate = da < 5.0 + expected_pkg.dataset[var] = da.where(~to_deactivate, -9999.0) + + cap_data = {} + mapping_ls = [ + ("rural_infiltration_capacity", "infiltration_capacity", 0), + ("urban_infiltration_capacity", "infiltration_capacity", 1), + ("rural_runoff_resistance", "upward_resistance", 0), + ("urban_runoff_resistance", "upward_resistance", 1), + ("rural_runon_resistance", "downward_resistance", 0), + ("urban_runon_resistance", "downward_resistance", 1), + ] + for cap_key, pkg_key, subunit_nr in mapping_ls: + cap_data[cap_key] = data_infiltration[pkg_key].sel( + subunit=subunit_nr, drop=True + ) + + imod5_data = {"cap": cap_data} + actual_pkg = Infiltration.from_imod5_data(imod5_data) + + ones_vars = ["bottom_resistance", "extra_storage_coefficient"] + expected = expected_pkg.dataset.drop_vars(ones_vars) + actual = actual_pkg.dataset.drop_vars(ones_vars) + + xr.testing.assert_equal(actual, expected) + + for var in ones_vars: + assert (actual_pkg.dataset[var] == 1.0).all() diff --git a/imod/tests/test_msw/test_ponding.py b/imod/tests/test_msw/test_ponding.py index e7bb47e06..94f0d64d0 100644 --- a/imod/tests/test_msw/test_ponding.py +++ b/imod/tests/test_msw/test_ponding.py @@ -116,7 +116,6 @@ def test_from_imod5_data(): # Create cap data cap_data = {} mapping_ls = [ - ("rural_runoff_resistance", "runoff_resistance", 0), ("rural_runoff_resistance", "runoff_resistance", 0), ("urban_runoff_resistance", "runoff_resistance", 1), ("rural_runon_resistance", "runon_resistance", 0), From 72c0255e9bf4f933871a87f407d9eebb3148a027 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Wed, 27 Nov 2024 10:14:07 +0100 Subject: [PATCH 04/13] Issue #1292 mete grids from imod5 (#1293) Fixes #1292 # Description Context: MetaSWAP can map meteorological grids, which have a coarser grid but a very fine time resolution, to its svats. The ``mete_grid.inp`` file contains paths to meteorological information, stored in ASCII grids, and there are separate mapping files ``"svat2precgrid.inp"`` and ``"svat2etrefgrid.inp"``. Nota bene: MetaSWAP has hardcoded filenames for packages, hence why I can search for "mete_grid.inp". Adds the following: - A ``MeteoGridCopy`` class to copy "mete_grid.inp" files, to avoid having to read and write the crazy amount of meteo files that exist in existing databases. - A ``MeteoMapping.from_imod5_data`` method which looks for the first ascii grid in the ``mete_grid.inp`` file, and reads that to derive mappings from meteorological cell to metaswap grid cell. - Add ``Imod5DataDict`` type alias, which can be implemented for other ``from_imod5_data`` methods in the future. --- docs/api/changelog.rst | 6 + imod/msw/__init__.py | 2 +- imod/msw/meteo_grid.py | 41 ++- imod/msw/meteo_mapping.py | 78 ++++- imod/msw/utilities/common.py | 9 + .../test_evapotranspiration_mapping.py | 226 +++++++------- imod/tests/test_msw/test_meteo_grid.py | 46 ++- .../test_msw/test_precipitation_mapping.py | 280 ++++++++---------- imod/typing/__init__.py | 1 + 9 files changed, 411 insertions(+), 278 deletions(-) diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst index 8c7aac2ed..d0ecceaef 100644 --- a/docs/api/changelog.rst +++ b/docs/api/changelog.rst @@ -9,6 +9,12 @@ The format is based on `Keep a Changelog`_, and this project adheres to [Unreleased] ------------ +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. + Changed ~~~~~~~ diff --git a/imod/msw/__init__.py b/imod/msw/__init__.py index bcbe174aa..83f8df90b 100644 --- a/imod/msw/__init__.py +++ b/imod/msw/__init__.py @@ -9,7 +9,7 @@ InitialConditionsSavedState, ) from imod.msw.landuse import LanduseOptions -from imod.msw.meteo_grid import MeteoGrid +from imod.msw.meteo_grid import MeteoGrid, MeteoGridCopy from imod.msw.meteo_mapping import EvapotranspirationMapping, PrecipitationMapping from imod.msw.model import MetaSwapModel from imod.msw.output_control import TimeOutputControl, VariableOutputControl diff --git a/imod/msw/meteo_grid.py b/imod/msw/meteo_grid.py index d8c283252..70d3625a5 100644 --- a/imod/msw/meteo_grid.py +++ b/imod/msw/meteo_grid.py @@ -1,6 +1,7 @@ import csv from pathlib import Path -from typing import Optional, Union +from shutil import copyfile +from typing import Optional, Union, cast import numpy as np import pandas as pd @@ -11,6 +12,9 @@ from imod.msw.pkgbase import MetaSwapPackage from imod.msw.regrid.regrid_schemes import MeteoGridRegridMethod from imod.msw.timeutil import to_metaswap_timeformat +from imod.msw.utilities.common import find_in_file_list +from imod.typing import Imod5DataDict +from imod.util.regrid_method_type import EmptyRegridMethod, RegridMethodType class MeteoGrid(MetaSwapPackage, IRegridPackage): @@ -188,3 +192,38 @@ def _pkgcheck(self): f"Received excess dims {excess_dims} in {self.__class__} for " f"{varname}, please provide data with {allowed_dims}" ) + + +class MeteoGridCopy(MetaSwapPackage, IRegridPackage): + """ + Class to copy existing ``mete_grid.inp``, which contains the meteorological + grid data. Next to a MeteoGridCopy instance, instances of + PrecipitationMapping and EvapotranspirationMapping are required as well to + specify meteorological information to MetaSWAP. + + Parameters + ---------- + path: Path to mete_grid.inp file + """ + + _file_name = "mete_grid.inp" + _meteo_dirname = "meteo_grids" + + _regrid_method: RegridMethodType = EmptyRegridMethod() + + def __init__(self, path: Path | str): + super().__init__() + self.dataset["path"] = path + + def write(self, directory: Path | str, *args): + directory = Path(directory) + path_metegrid = Path(str(self.dataset["path"].values[()])) + new_path = directory / self._file_name + copyfile(path_metegrid, new_path) + + @classmethod + def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "MeteoGridCopy": + paths = cast(list[str], imod5_data["extra"]["paths"]) + filepath = find_in_file_list(cls._file_name, paths) + + return cls(filepath) diff --git a/imod/msw/meteo_mapping.py b/imod/msw/meteo_mapping.py index 0fa3c818c..1de00fa69 100644 --- a/imod/msw/meteo_mapping.py +++ b/imod/msw/meteo_mapping.py @@ -1,18 +1,68 @@ from copy import deepcopy -from typing import Any, Optional, TextIO +from pathlib import Path +from textwrap import dedent +from typing import Any, Optional, TextIO, cast import numpy as np import pandas as pd import xarray as xr +import imod from imod.mf6.utilities.regrid import RegridderWeightsCache from imod.msw.fixed_format import VariableMetaData from imod.msw.pkgbase import MetaSwapPackage +from imod.msw.utilities.common import find_in_file_list from imod.prepare import common -from imod.typing import GridDataArray, IntArray +from imod.typing import GridDataArray, Imod5DataDict, IntArray from imod.util.regrid_method_type import RegridMethodType +def _is_parsable_and_existing_path(potential_path: str, mete_grid_path: Path) -> bool: + """ + mete_grid.inp can contain values like "0.", which are converted to float by + MetaSWAP. String is converted to path and checked if existing path. + """ + try: + float(potential_path) + return False + except ValueError: + # Resolve paths relative to mete_grid.inp path. + path = mete_grid_path / ".." / Path(potential_path) + return path.is_file() + + +def open_first_meteo_grid(mete_grid_path: str | Path, column_nr: int) -> xr.DataArray: + """ + Find and open first meteo grid path in mete_grid.inp. This grid is enough to + generate meteomappings. There can be floats before in the column which + should be skipped. + """ + if column_nr not in [2, 3]: + raise ValueError("Column nr should be 2 or 3") + + mete_grid_path = Path(mete_grid_path) + with open(mete_grid_path, "r") as f: + lines = f.readlines() + + potential_paths = [line.split(",")[column_nr].replace('"', "") for line in lines] + for potential_path in potential_paths: + if _is_parsable_and_existing_path(potential_path, mete_grid_path): + resolved_path = mete_grid_path / ".." / Path(potential_path) + return imod.rasterio.open(resolved_path) + + error_message = dedent(f""" + Did not find parsable path to existing .ASC file in column {column_nr}. Got + values (printing first 10): {potential_paths[:10]}.""") + + raise ValueError(error_message) + + +def open_first_meteo_grid_from_imod5_data(imod5_data: Imod5DataDict, column_nr: int): + paths = cast(list[str], imod5_data["extra"]["paths"]) + metegrid_path = find_in_file_list("mete_grid.inp", paths) + return open_first_meteo_grid(metegrid_path, column_nr=column_nr) + + class MeteoMapping(MetaSwapPackage): """ This class provides common methods for creating mappings between @@ -125,6 +175,18 @@ def __init__( super().__init__() self.meteo = precipitation + @classmethod + def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "PrecipitationMapping": + """ + Construct precipitation mapping from imod5 data. Opens first ascii grid + in mete_grid.inp, which is used to construct mappings to svats. The + grids should not change in dimension over time. No checks are done + whether cells switch from inactive to active or vice versa. + """ + column_nr = 2 + meteo_grid = open_first_meteo_grid_from_imod5_data(imod5_data, column_nr) + return cls(meteo_grid) + class EvapotranspirationMapping(MeteoMapping): """ @@ -156,3 +218,15 @@ def __init__( ): super().__init__() self.meteo = evapotranspiration + + @classmethod + def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "EvapotranspirationMapping": + """ + Construct evapotranspiration mapping from imod5 data. Opens first ascii + grid in mete_grid.inp, which is used to construct mappings to svats. The + grids should not change in dimension over time. No checks are done + whether cells switch from inactive to active or vice versa. + """ + column_nr = 3 + meteo_grid = open_first_meteo_grid_from_imod5_data(imod5_data, column_nr) + return cls(meteo_grid) diff --git a/imod/msw/utilities/common.py b/imod/msw/utilities/common.py index 836902662..db69aab33 100644 --- a/imod/msw/utilities/common.py +++ b/imod/msw/utilities/common.py @@ -1,6 +1,15 @@ +from pathlib import Path + from imod.typing import GridDataArray from imod.typing.grid import concat def concat_imod5(arg1: GridDataArray, arg2: GridDataArray) -> GridDataArray: return concat([arg1, arg2], dim="subunit").assign_coords(subunit=[0, 1]) + + +def find_in_file_list(filename: str, paths: list[str]) -> str: + for file in paths: + if filename == Path(file[0]).name.lower(): + return file[0] + raise ValueError(f"could not find {filename} in list of paths: {paths}") diff --git a/imod/tests/test_msw/test_evapotranspiration_mapping.py b/imod/tests/test_msw/test_evapotranspiration_mapping.py index 0d5278f58..b2752a60f 100644 --- a/imod/tests/test_msw/test_evapotranspiration_mapping.py +++ b/imod/tests/test_msw/test_evapotranspiration_mapping.py @@ -1,3 +1,4 @@ +import re import tempfile from pathlib import Path @@ -5,33 +6,13 @@ import pytest import xarray as xr from numpy.testing import assert_equal +from pytest_cases import parametrize_with_cases from imod import msw -def test_evapotranspiration_mapping_simple(fixed_format_parser): - x_meteo = [-0.5, 1.5, 3.5] - y_meteo = [0.5, 2.5, 4.5] - subunit_meteo = [0, 1] - dx_meteo = 2.0 - dy_meteo = 2.0 - # fmt: off - evapotranspiration = xr.DataArray( - np.array( - [ - [[1.0, 1.0, 1.0], - [1.0, 1.0, 1.0], - [1.0, 1.0, 1.0]], - - [[1.0, 1.0, 1.0], - [1.0, 1.0, 1.0], - [1.0, 1.0, 1.0]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit_meteo, "y": y_meteo, "x": x_meteo, "dx": dx_meteo, "dy": dy_meteo} - ) - # fmt: on +@pytest.fixture(scope="function") +def svat_index() -> xr.DataArray: x_svat = [1.0, 2.0, 3.0] y_svat = [1.0, 2.0, 3.0] subunit_svat = [0, 1] @@ -56,31 +37,12 @@ def test_evapotranspiration_mapping_simple(fixed_format_parser): ) # fmt: on index = (svat != 0).values.ravel() + return svat, index - evapotranspiration_mapping = msw.EvapotranspirationMapping(evapotranspiration) - with tempfile.TemporaryDirectory() as output_dir: - output_dir = Path(output_dir) - evapotranspiration_mapping.write(output_dir, index, svat, None, None) - - results = fixed_format_parser( - output_dir / msw.EvapotranspirationMapping._file_name, - msw.EvapotranspirationMapping._metadata_dict, - ) - - assert_equal(results["svat"], np.array([1, 2, 3, 4, 5, 6])) - assert_equal(results["row"], np.array([1, 2, 1, 2, 2, 2])) - assert_equal(results["column"], np.array([2, 2, 2, 2, 2, 3])) - - -def test_evapotranspiration_mapping_negative_dx(fixed_format_parser): - x_meteo = [3.5, 1.5, -0.5] - y_meteo = [0.5, 2.5, 4.5] - subunit_meteo = [0, 1] - dx_meteo = -2.0 - dy_meteo = 2.0 +def create_meteo_grid(x, y, subunit, dx, dy): # fmt: off - evapotranspiration = xr.DataArray( + return xr.DataArray( np.array( [ [[1.0, 1.0, 1.0], @@ -93,34 +55,20 @@ def test_evapotranspiration_mapping_negative_dx(fixed_format_parser): ] ), dims=("subunit", "y", "x"), - coords={"subunit": subunit_meteo, "y": y_meteo, "x": x_meteo, "dx": dx_meteo, "dy": dy_meteo} + coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} ) # fmt: on - x_svat = [1.0, 2.0, 3.0] - y_svat = [1.0, 2.0, 3.0] - subunit_svat = [0, 1] - dx_svat = 1.0 - dy_svat = 1.0 - # fmt: off - svat = xr.DataArray( - np.array( - [ - [[0, 1, 0], - [0, 0, 0], - [0, 2, 0]], - [[0, 3, 0], - [4, 5, 6], - [0, 0, 0]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit_svat, "y": y_svat, "x": x_svat, "dx": dx_svat, "dy": dy_svat} - ) - # fmt: on - index = (svat != 0).values.ravel() +def test_evapotranspiration_mapping_simple(fixed_format_parser, svat_index): + svat, index = svat_index + x = [-0.5, 1.5, 3.5] + y = [0.5, 2.5, 4.5] + subunit = [0, 1] + dx = 2.0 + dy = 2.0 + evapotranspiration = create_meteo_grid(x, y, subunit, dx, dy) evapotranspiration_mapping = msw.EvapotranspirationMapping(evapotranspiration) with tempfile.TemporaryDirectory() as output_dir: @@ -134,58 +82,45 @@ def test_evapotranspiration_mapping_negative_dx(fixed_format_parser): assert_equal(results["svat"], np.array([1, 2, 3, 4, 5, 6])) assert_equal(results["row"], np.array([1, 2, 1, 2, 2, 2])) - assert_equal(results["column"], np.array([2, 2, 2, 2, 2, 1])) + assert_equal(results["column"], np.array([2, 2, 2, 2, 2, 3])) -def test_evapotranspiration_mapping_out_of_bound(): - x_meteo = [-2.5, -0.5, 1.5] - y_meteo = [0.5, 2.5, 4.5] - subunit_meteo = [0, 1] - dx_meteo = 2.0 - dy_meteo = 2.0 +def test_evapotranspiration_mapping_negative_dx(fixed_format_parser, svat_index): + svat, index = svat_index - # fmt: off - evapotranspiration = xr.DataArray( - np.array( - [ - [[1.0, 1.0, 1.0], - [1.0, 1.0, 1.0], - [1.0, 1.0, 1.0]], + x = [3.5, 1.5, -0.5] + y = [0.5, 2.5, 4.5] + subunit = [0, 1] + dx = -2.0 + dy = 2.0 - [[1.0, 1.0, 1.0], - [1.0, 1.0, 1.0], - [1.0, 1.0, 1.0]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit_meteo, "y": y_meteo, "x": x_meteo, "dx": dx_meteo, "dy": dy_meteo} - ) - # fmt: on + evapotranspiration = create_meteo_grid(x, y, subunit, dx, dy) + evapotranspiration_mapping = msw.EvapotranspirationMapping(evapotranspiration) - x_svat = [1.0, 2.0, 3.0] - y_svat = [1.0, 2.0, 3.0] - subunit_svat = [0, 1] - dx_svat = 1.0 - dy_svat = 1.0 - # fmt: off - svat = xr.DataArray( - np.array( - [ - [[0, 1, 0], - [0, 0, 0], - [0, 2, 0]], + with tempfile.TemporaryDirectory() as output_dir: + output_dir = Path(output_dir) + evapotranspiration_mapping.write(output_dir, index, svat, None, None) - [[0, 3, 0], - [4, 5, 6], - [0, 0, 0]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit_svat, "y": y_svat, "x": x_svat, "dx": dx_svat, "dy": dy_svat} - ) - # fmt: on - index = (svat != 0).values.ravel() + results = fixed_format_parser( + output_dir / msw.EvapotranspirationMapping._file_name, + msw.EvapotranspirationMapping._metadata_dict, + ) + assert_equal(results["svat"], np.array([1, 2, 3, 4, 5, 6])) + assert_equal(results["row"], np.array([1, 2, 1, 2, 2, 2])) + assert_equal(results["column"], np.array([2, 2, 2, 2, 2, 1])) + + +def test_evapotranspiration_mapping_out_of_bound(svat_index): + svat, index = svat_index + + x = [-2.5, -0.5, 1.5] + y = [0.5, 2.5, 4.5] + subunit = [0, 1] + dx = 2.0 + dy = 2.0 + + evapotranspiration = create_meteo_grid(x, y, subunit, dx, dy) evapotranspiration_mapping = msw.EvapotranspirationMapping(evapotranspiration) with tempfile.TemporaryDirectory() as output_dir: @@ -193,3 +128,68 @@ def test_evapotranspiration_mapping_out_of_bound(): # The grid is out of bounds, which is why we expect a ValueError to be raisen with pytest.raises(ValueError): evapotranspiration_mapping.write(output_dir, index, svat, None, None) + + +def setup_meteo_grid(datadir): + """Setup precipitation grid and write mete_grid.inp""" + # Arrange + x = [-0.5, 1.5, 3.5] + y = [4.5, 2.5, 0.5] + subunit = [0, 1] + dx = 2.0 + dy = -2.0 + + time = [np.datetime64(t) for t in ["2001-01-01", "2001-01-02", "2001-01-03"]] + time_da = xr.DataArray([1.0, 1.0, 1.0], coords={"time": time}) + + precipitation = create_meteo_grid(x, y, subunit, dx, dy).isel(subunit=0, drop=True) + precipitation_times = time_da * precipitation + mete_grid = msw.MeteoGrid(precipitation_times, precipitation_times) + mete_grid.write(datadir) + return precipitation + + +class MeteGridCases: + """ + Cases return strings to replace in mete_grid.inp, to set paths to floats. + """ + + def case_all_paths(self): + return r"nothing_to_replace" + + def case_some_paths(self): + return r'"meteo_grids\\(\w+)_20010101000000.asc"' + + def case_no_paths(self): + return r'"meteo_grids\\(\w+)_([0-9]+).asc"' + + +@parametrize_with_cases("replace_string", cases=MeteGridCases) +def test_evapotranspiration_from_imod5(tmpdir_factory, replace_string, request): + datadir = Path(tmpdir_factory.mktemp("evapotranspiration_mapping")) + evapotranspiration = setup_meteo_grid(datadir) + + paths = [["foo"], [datadir / "mete_grid.inp"], ["bar"]] + imod5_data = {"extra": {"paths": paths}} + + # Replace text in existing mete_grid.inp + with open(datadir / "mete_grid.inp", "r") as f: + text = f.read() + text_replaced = re.sub(replace_string, '"0.0"', text) + with open(datadir / "mete_grid.inp", "w") as f: + f.write(text_replaced) + + # Act + if request.node.callspec.id == "no_paths": + with pytest.raises(ValueError): + msw.EvapotranspirationMapping.from_imod5_data(imod5_data) + + else: + evapotranspiration_mapping = msw.EvapotranspirationMapping.from_imod5_data( + imod5_data + ) + actual = evapotranspiration_mapping.meteo + + # Assert + assert len(actual.coords["time"]) == 1 + xr.testing.assert_equal(evapotranspiration, actual.isel(time=0, drop=True)) diff --git a/imod/tests/test_msw/test_meteo_grid.py b/imod/tests/test_msw/test_meteo_grid.py index 8cca2cc06..b5108842e 100644 --- a/imod/tests/test_msw/test_meteo_grid.py +++ b/imod/tests/test_msw/test_meteo_grid.py @@ -1,4 +1,5 @@ import csv +import filecmp import os import tempfile from pathlib import Path @@ -11,7 +12,7 @@ from numpy.testing import assert_equal from imod.mf6.utilities.regrid import RegridderWeightsCache -from imod.msw import MeteoGrid +from imod.msw import MeteoGrid, MeteoGridCopy def setup_meteo_grid(): @@ -123,3 +124,46 @@ def test_regrid_meteo(simple_2d_grid_with_subunits): assert np.all(regridded_ponding.dataset["x"].values == new_grid["x"].values) assert np.all(regridded_ponding.dataset["y"].values == new_grid["y"].values) + + +def test_meteogridcopy_write(): + # Arrange + meteo_grid = setup_meteo_grid() + + with tempfile.TemporaryDirectory() as output_dir: + grid_dir = Path(output_dir) / "grid" + grid_dir.mkdir(exist_ok=True, parents=True) + meteo_grid.write(grid_dir) + + meteo_grid_copy = MeteoGridCopy(grid_dir / "mete_grid.inp") + copy_dir = Path(output_dir) / "copied" + copy_dir.mkdir(exist_ok=True, parents=True) + # Act + meteo_grid_copy.write(copy_dir) + # Assert + assert filecmp.cmp(grid_dir / "mete_grid.inp", copy_dir / "mete_grid.inp") + + +def test_meteogridcopy_from_imod5(): + meteo_grid = setup_meteo_grid() + + with tempfile.TemporaryDirectory() as output_dir: + grid_dir = Path(output_dir) / "grid" + grid_dir.mkdir(exist_ok=True, parents=True) + meteo_grid.write(grid_dir) + + imod5_data = {} + imod5_data["extra"] = {} + imod5_data["extra"]["paths"] = [ + ["foo"], + [(grid_dir / "mete_grid.inp").resolve()], + ["bar"], + ] + + meteo_grid_copy = MeteoGridCopy.from_imod5_data(imod5_data) + copy_dir = Path(output_dir) / "copied" + copy_dir.mkdir(exist_ok=True, parents=True) + # Act + meteo_grid_copy.write(copy_dir) + # Assert + assert filecmp.cmp(grid_dir / "mete_grid.inp", copy_dir / "mete_grid.inp") diff --git a/imod/tests/test_msw/test_precipitation_mapping.py b/imod/tests/test_msw/test_precipitation_mapping.py index 3d74ac306..e9162fdf0 100644 --- a/imod/tests/test_msw/test_precipitation_mapping.py +++ b/imod/tests/test_msw/test_precipitation_mapping.py @@ -1,3 +1,4 @@ +import re import tempfile from pathlib import Path @@ -5,34 +6,13 @@ import pytest import xarray as xr from numpy.testing import assert_equal +from pytest_cases import parametrize_with_cases from imod import msw -def test_precipitation_mapping_simple(fixed_format_parser): - x_meteo = [-0.5, 1.5, 3.5] - y_meteo = [0.5, 2.5, 4.5] - subunit_meteo = [0, 1] - dx_meteo = 2.0 - dy_meteo = 2.0 - # fmt: off - precipitation = xr.DataArray( - np.array( - [ - [[1.0, 1.0, 1.0], - [1.0, 1.0, 1.0], - [1.0, 1.0, 1.0]], - - [[1.0, 1.0, 1.0], - [1.0, 1.0, 1.0], - [1.0, 1.0, 1.0]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit_meteo, "y": y_meteo, "x": x_meteo, "dx": dx_meteo, "dy": dy_meteo} - ) - # fmt: on - +@pytest.fixture(scope="function") +def svat_index() -> xr.DataArray: x_svat = [1.0, 2.0, 3.0] y_svat = [1.0, 2.0, 3.0] subunit_svat = [0, 1] @@ -57,31 +37,12 @@ def test_precipitation_mapping_simple(fixed_format_parser): ) # fmt: on index = (svat != 0).values.ravel() - - precipitation_mapping = msw.PrecipitationMapping(precipitation) - - with tempfile.TemporaryDirectory() as output_dir: - output_dir = Path(output_dir) - precipitation_mapping.write(output_dir, index, svat, None, None) - - results = fixed_format_parser( - output_dir / msw.PrecipitationMapping._file_name, - msw.PrecipitationMapping._metadata_dict, - ) - - assert_equal(results["svat"], np.array([1, 2, 3, 4, 5, 6])) - assert_equal(results["row"], np.array([1, 2, 1, 2, 2, 2])) - assert_equal(results["column"], np.array([2, 2, 2, 2, 2, 3])) + return svat, index -def test_precipitation_mapping_negative_dy_meteo(fixed_format_parser): - x_meteo = [-0.5, 1.5, 3.5] - y_meteo = [4.5, 2.5, 0.5] - subunit_meteo = [0, 1] - dx_meteo = 2.0 - dy_meteo = -2.0 +def create_meteo_grid(x, y, subunit, dx, dy): # fmt: off - precipitation = xr.DataArray( + return xr.DataArray( np.array( [ [[1.0, 1.0, 1.0], @@ -94,35 +55,21 @@ def test_precipitation_mapping_negative_dy_meteo(fixed_format_parser): ] ), dims=("subunit", "y", "x"), - coords={"subunit": subunit_meteo, "y": y_meteo, "x": x_meteo, "dx": dx_meteo, "dy": dy_meteo} + coords={"subunit": subunit, "y": y, "x": x, "dx": dx, "dy": dy} ) # fmt: on - x_svat = [1.0, 2.0, 3.0] - y_svat = [1.0, 2.0, 3.0] - subunit_svat = [0, 1] - dx_svat = 1.0 - dy_svat = 1.0 - # fmt: off - svat = xr.DataArray( - np.array( - [ - [[0, 1, 0], - [0, 0, 0], - [0, 2, 0]], +def test_precipitation_mapping_simple(fixed_format_parser, svat_index): + svat, index = svat_index - [[0, 3, 0], - [4, 5, 6], - [0, 0, 0]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit_svat, "y": y_svat, "x": x_svat, "dx": dx_svat, "dy": dy_svat} - ) - # fmt: on - index = (svat != 0).values.ravel() + x = [-0.5, 1.5, 3.5] + y = [0.5, 2.5, 4.5] + subunit = [0, 1] + dx = 2.0 + dy = 2.0 + precipitation = create_meteo_grid(x, y, subunit, dx, dy) precipitation_mapping = msw.PrecipitationMapping(precipitation) with tempfile.TemporaryDirectory() as output_dir: @@ -135,58 +82,22 @@ def test_precipitation_mapping_negative_dy_meteo(fixed_format_parser): ) assert_equal(results["svat"], np.array([1, 2, 3, 4, 5, 6])) - assert_equal(results["row"], np.array([3, 2, 3, 2, 2, 2])) + assert_equal(results["row"], np.array([1, 2, 1, 2, 2, 2])) assert_equal(results["column"], np.array([2, 2, 2, 2, 2, 3])) -def test_precipitation_mapping_negative_dy_meteo_svat(fixed_format_parser): - x_meteo = [-0.5, 1.5, 3.5] - y_meteo = [4.5, 2.5, 0.5] - subunit_meteo = [0, 1] - dx_meteo = 2.0 - dy_meteo = -2.0 - # fmt: off - precipitation = xr.DataArray( - np.array( - [ - [[1.0, 1.0, 1.0], - [1.0, 1.0, 1.0], - [1.0, 1.0, 1.0]], +def test_precipitation_mapping_negative_dy_meteo(fixed_format_parser, svat_index): + svat, index = svat_index - [[1.0, 1.0, 1.0], - [1.0, 1.0, 1.0], - [1.0, 1.0, 1.0]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit_meteo, "y": y_meteo, "x": x_meteo, "dx": dx_meteo, "dy": dy_meteo} - ) - # fmt: on - x_svat = [1.0, 2.0, 3.0] - y_svat = [3.0, 2.0, 1.0] - subunit_svat = [0, 1] - dx_svat = 1.0 - dy_svat = 1.0 + x = [-0.5, 1.5, 3.5] + y = [4.5, 2.5, 0.5] + subunit = [0, 1] + dx = 2.0 + dy = -2.0 - # fmt: off - svat = xr.DataArray( - np.array( - [ - [[0, 1, 0], - [0, 0, 0], - [0, 2, 0]], - - [[0, 3, 0], - [4, 5, 6], - [0, 0, 0]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit_svat, "y": y_svat, "x": x_svat, "dx": dx_svat, "dy": dy_svat} - ) - # fmt: on - index = (svat != 0).values.ravel() + precipitation = create_meteo_grid(x, y, subunit, dx, dy) precipitation_mapping = msw.PrecipitationMapping(precipitation) + with tempfile.TemporaryDirectory() as output_dir: output_dir = Path(output_dir) precipitation_mapping.write(output_dir, index, svat, None, None) @@ -197,60 +108,46 @@ def test_precipitation_mapping_negative_dy_meteo_svat(fixed_format_parser): ) assert_equal(results["svat"], np.array([1, 2, 3, 4, 5, 6])) - assert_equal(results["row"], np.array([2, 3, 2, 2, 2, 2])) + assert_equal(results["row"], np.array([3, 2, 3, 2, 2, 2])) assert_equal(results["column"], np.array([2, 2, 2, 2, 2, 3])) -def test_precipitation_mapping_out_of_bound(): - x_meteo = [-0.5, 1.5, 3.5] - y_meteo = [2.5, 4.5, 6.5] - subunit_meteo = [0, 1] - dx_meteo = 2.0 - dy_meteo = 2.0 +def test_precipitation_mapping_negative_dy_meteo_svat(fixed_format_parser, svat_index): + svat, index = svat_index + svat = svat.assign_coords(y=[3.0, 2.0, 1.0], dy=-1.0) - # fmt: off - precipitation = xr.DataArray( - np.array( - [ - [[1.0, 1.0, 1.0], - [1.0, 1.0, 1.0], - [1.0, 1.0, 1.0]], + x = [-0.5, 1.5, 3.5] + y = [4.5, 2.5, 0.5] + subunit = [0, 1] + dx = 2.0 + dy = -2.0 - [[1.0, 1.0, 1.0], - [1.0, 1.0, 1.0], - [1.0, 1.0, 1.0]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit_meteo, "y": y_meteo, "x": x_meteo, "dx": dx_meteo, "dy": dy_meteo} - ) - # fmt: on + precipitation = create_meteo_grid(x, y, subunit, dx, dy) + precipitation_mapping = msw.PrecipitationMapping(precipitation) + with tempfile.TemporaryDirectory() as output_dir: + output_dir = Path(output_dir) + precipitation_mapping.write(output_dir, index, svat, None, None) - x_svat = [1.0, 2.0, 3.0] - y_svat = [1.0, 2.0, 3.0] - subunit_svat = [0, 1] - dx_svat = 1.0 - dy_svat = 1.0 + results = fixed_format_parser( + output_dir / msw.PrecipitationMapping._file_name, + msw.PrecipitationMapping._metadata_dict, + ) - # fmt: off - svat = xr.DataArray( - np.array( - [ - [[0, 1, 0], - [0, 0, 0], - [0, 2, 0]], + assert_equal(results["svat"], np.array([1, 2, 3, 4, 5, 6])) + assert_equal(results["row"], np.array([2, 3, 2, 2, 2, 2])) + assert_equal(results["column"], np.array([2, 2, 2, 2, 2, 3])) - [[0, 3, 0], - [4, 5, 6], - [0, 0, 0]], - ] - ), - dims=("subunit", "y", "x"), - coords={"subunit": subunit_svat, "y": y_svat, "x": x_svat, "dx": dx_svat, "dy": dy_svat} - ) - # fmt: on - index = (svat != 0).values.ravel() +def test_precipitation_mapping_out_of_bound(svat_index): + svat, index = svat_index + + x = [-0.5, 1.5, 3.5] + y = [2.5, 4.5, 6.5] + subunit = [0, 1] + dx = 2.0 + dy = 2.0 + + precipitation = create_meteo_grid(x, y, subunit, dx, dy) precipitation_mapping = msw.PrecipitationMapping(precipitation) with tempfile.TemporaryDirectory() as output_dir: @@ -258,3 +155,66 @@ def test_precipitation_mapping_out_of_bound(): # The grid is out of bounds, which is why we expect a ValueError to be raisen with pytest.raises(ValueError): precipitation_mapping.write(output_dir, index, svat, None, None) + + +def setup_meteo_grid(datadir): + """Setup precipitation grid and write mete_grid.inp""" + # Arrange + x = [-0.5, 1.5, 3.5] + y = [4.5, 2.5, 0.5] + subunit = [0, 1] + dx = 2.0 + dy = -2.0 + + time = [np.datetime64(t) for t in ["2001-01-01", "2001-01-02", "2001-01-03"]] + time_da = xr.DataArray([1.0, 1.0, 1.0], coords={"time": time}) + + precipitation = create_meteo_grid(x, y, subunit, dx, dy).isel(subunit=0, drop=True) + precipitation_times = time_da * precipitation + mete_grid = msw.MeteoGrid(precipitation_times, precipitation_times) + mete_grid.write(datadir) + return precipitation + + +class MeteGridCases: + """ + Cases return strings to replace in mete_grid.inp, to set paths to floats. + """ + + def case_all_paths(self): + return r"nothing_to_replace" + + def case_some_paths(self): + return r'"meteo_grids\\(\w+)_20010101000000.asc"' + + def case_no_paths(self): + return r'"meteo_grids\\(\w+)_([0-9]+).asc"' + + +@parametrize_with_cases("replace_string", cases=MeteGridCases) +def test_precipitation_from_imod5(tmpdir_factory, replace_string, request): + # Arrange + datadir = Path(tmpdir_factory.mktemp("precipitation_mapping")) + precipitation = setup_meteo_grid(datadir) + paths = [["foo"], [datadir / "mete_grid.inp"], ["bar"]] + imod5_data = {"extra": {"paths": paths}} + + # Replace text in existing mete_grid.inp + with open(datadir / "mete_grid.inp", "r") as f: + text = f.read() + text_replaced = re.sub(replace_string, '"0.0"', text) + with open(datadir / "mete_grid.inp", "w") as f: + f.write(text_replaced) + + # Act + if request.node.callspec.id == "no_paths": + with pytest.raises(ValueError): + msw.PrecipitationMapping.from_imod5_data(imod5_data) + + else: + precipitation_mapping = msw.PrecipitationMapping.from_imod5_data(imod5_data) + actual = precipitation_mapping.meteo + + # Assert + assert len(actual.coords["time"]) == 1 + xr.testing.assert_equal(precipitation, actual.isel(time=0, drop=True)) diff --git a/imod/typing/__init__.py b/imod/typing/__init__.py index 760346894..8beb7fd13 100644 --- a/imod/typing/__init__.py +++ b/imod/typing/__init__.py @@ -12,6 +12,7 @@ GridDataArray: TypeAlias = Union[xr.DataArray, xu.UgridDataArray] GridDataset: TypeAlias = Union[xr.Dataset, xu.UgridDataset] GridDataDict: TypeAlias = dict[str, GridDataArray] +Imod5DataDict: TypeAlias = dict[str, GridDataDict | dict[str, list[str]]] ScalarAsDataArray: TypeAlias = Union[xr.DataArray, xu.UgridDataArray] ScalarAsDataset: TypeAlias = Union[xr.Dataset, xu.UgridDataset] UnstructuredData: TypeAlias = Union[xu.UgridDataset, xu.UgridDataArray] From 8461d634bb877c274ceadd6c9c7c85fb10b5263f Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Tue, 3 Dec 2024 11:20:02 +0100 Subject: [PATCH 05/13] Issue #1313 recharge from imod5 cap data (#1315) Fixes #1313 # Description Construct an rch-package from iMOD5 data in the CAP package, loaded with the ``open_projectfile_data`` function. Package is used to couple MODFLOW6 to MetaSWAP models. Active cells will have a recharge rate of 0.0. At the moment, MetaSWAP can only be coupled to the first layer, as this is also the case for ``primod``. iMOD Coupler these days supports coupling to other layers as well, but ``primod`` doesn't. Picking this up for iMOD Python and ``primod`` is worthy a separate story. In detail: - Add ``Recharge.from_imod5_cap_data`` class method, to construct an empty Recharge package with 0.0 rate. - Minor refactor: Put ``GridData`` helper functions to ``msw/utilities/imod5_converter.py``, so that they can be reused. # Checklist - [x] Links to correct issue - [x] Update changelog, if changes affect users - [x] PR title starts with ``Issue #nr``, e.g. ``Issue #737`` - [x] Unit tests were added - [ ] **If feature added**: Added/extended example --- docs/api/changelog.rst | 2 + imod/mf6/model_gwf.py | 3 + imod/mf6/rch.py | 39 ++++++++++-- imod/msw/grid_data.py | 71 +++------------------- imod/msw/utilities/imod5_converter.py | 87 +++++++++++++++++++++++++++ imod/tests/test_mf6/test_mf6_rch.py | 27 +++++++++ 6 files changed, 163 insertions(+), 66 deletions(-) create mode 100644 imod/msw/utilities/imod5_converter.py 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 ~~~~~ 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"] diff --git a/imod/mf6/rch.py b/imod/mf6/rch.py index ba917cbd8..392cc9d56 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 GridDataArray, GridDataDict, Imod5DataDict from imod.typing.grid import ( enforce_dim_order, is_planar_grid, @@ -218,7 +223,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, ) @@ -227,7 +232,33 @@ 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: 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. 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"]) + + msw_area = get_cell_area_from_imod5_data(cap_data) + active, _ = is_msw_active_cell(target_dis, cap_data, msw_area) + + data = {} + 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) @classmethod def get_regrid_methods(cls) -> RechargeRegridMethod: 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 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 68b5ef2b6daaf3230ef0593fb4b59081cf3a5636 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Tue, 3 Dec 2024 16:49:29 +0100 Subject: [PATCH 06/13] Issue #1311 sprinkling from imod5 data (#1312) Fixes #1311 # Description Adds the following: - ``imod.mf6.LayeredWell.from_imod5_cap_data``, to construct LayeredWells for MODLOW6 from the CAP package in iMOD5 Data. Throws an error if data is provided as point data (IPF, which is read as pandas DataFrame). This has rate 0.0, actual rates will be inserted by in the coupling scheme of iMOD Coupler. - ``imod.mf6.GroundwaterFlowModel.from_imod5_data`` adds Sprinkling well if CAP package present. This does not affect MODFLOW6 model if not coupled to MetaSWAP: The wells have rates 0.0, so will not affect the simulation results. - ``Sprinkling.from_imod5_data`` to construct a Sprinkling package from iMOD5 data. # Checklist - [x] Links to correct issue - [x] Update changelog, if changes affect users - [x] PR title starts with ``Issue #nr``, e.g. ``Issue #737`` - [x] Unit tests were added - [ ] **If feature added**: Added/extended example --- docs/api/changelog.rst | 6 ++ docs/api/mf6.rst | 1 + imod/mf6/model_gwf.py | 4 +- imod/mf6/utilities/imod5_converter.py | 73 ++++++++++++++++++- imod/mf6/wel.py | 43 +++++++++++- imod/msw/sprinkling.py | 97 +++++++++++++++++++++++++- imod/tests/conftest.py | 4 ++ imod/tests/fixtures/imod5_cap_data.py | 60 ++++++++++++++++ imod/tests/test_mf6/test_mf6_wel.py | 22 ++++++ imod/tests/test_msw/test_sprinkling.py | 37 ++++++++++ imod/typing/__init__.py | 8 ++- 11 files changed, 348 insertions(+), 7 deletions(-) create mode 100644 imod/tests/fixtures/imod5_cap_data.py 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 From 29886af18e53f8a1fc5b9b6abbed2a75bb617356 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Wed, 4 Dec 2024 13:00:19 +0100 Subject: [PATCH 07/13] Issue #1320 Add CopyFiles package (#1321) Fixes #1320 # Description The iMOD5 projectfile contains an "EXTRA" section, with files to copy. These files usually contain settings and lookup tables for MetaSWAP, which have no direct relation to cell data. This PR adds CopyFiles package to copy these files. # Checklist - [x] Links to correct issue - [x] Update changelog, if changes affect users - [x] PR title starts with ``Issue #nr``, e.g. ``Issue #737`` - [x] Unit tests were added - [ ] **If feature added**: Added/extended example --- docs/api/changelog.rst | 2 + imod/msw/__init__.py | 1 + imod/msw/copy_files.py | 55 ++++++++++++++++++++ imod/tests/test_msw/test_copy_files.py | 69 ++++++++++++++++++++++++++ 4 files changed, 127 insertions(+) create mode 100644 imod/msw/copy_files.py create mode 100644 imod/tests/test_msw/test_copy_files.py diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst index 531a88116..1f073881c 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. +- :class:`imod.msw.CopyFiles` to copy settings and lookup tables in existing + ``.inp`` files. - :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. diff --git a/imod/msw/__init__.py b/imod/msw/__init__.py index 83f8df90b..d8d053c49 100644 --- a/imod/msw/__init__.py +++ b/imod/msw/__init__.py @@ -1,3 +1,4 @@ +from imod.msw.copy_files import FileCopier from imod.msw.coupler_mapping import CouplerMapping from imod.msw.grid_data import GridData from imod.msw.idf_mapping import IdfMapping diff --git a/imod/msw/copy_files.py b/imod/msw/copy_files.py new file mode 100644 index 000000000..c5615b5fb --- /dev/null +++ b/imod/msw/copy_files.py @@ -0,0 +1,55 @@ +from pathlib import Path +from shutil import copy2 +from typing import cast + +import numpy as np +import xarray as xr + +from imod.logging import logger +from imod.logging.loglevel import LogLevel +from imod.msw.pkgbase import MetaSwapPackage +from imod.typing import Imod5DataDict + +_LOG_MESSAGE_TEMPLATE = """\ +Will not copy files {filtered}, these will be generated by iMOD Python +instead.""" + + +class FileCopier(MetaSwapPackage): + def __init__(self, paths: list[str]): + super().__init__() + paths_da = xr.DataArray( + paths, coords={"file_nr": np.arange(len(paths))}, dims=("file_nr",) + ) + self.dataset["paths"] = paths_da + + @classmethod + def from_imod5_data(cls, imod5_data: Imod5DataDict): + paths = cast(list[list[str]], imod5_data["extra"]["paths"]) + paths_unpacked = {Path(p[0]) for p in paths} + files_to_filter = ( + "mete_grid.inp", + "para_sim.inp", + "svat2precgrid.inp", + "svat2etrefgrid.inp", + ) + paths_included = [ + str(p) for p in paths_unpacked if p.name.lower() not in files_to_filter + ] + paths_excluded = {str(p) for p in paths_unpacked} - set(paths_included) + if paths_excluded: + log_message = _LOG_MESSAGE_TEMPLATE.format(filtered=paths_excluded) + logger.log( + loglevel=LogLevel.INFO, + message=log_message, + ) + return cls(paths_included) + + def write(self, directory: str | Path, *_): + directory = Path(directory) + + src_paths = [Path(p) for p in self.dataset["paths"].to_numpy()] + dst_paths = [directory / p.name for p in src_paths] + + for src_path, dst_path in zip(src_paths, dst_paths): + copy2(src_path, dst_path) diff --git a/imod/tests/test_msw/test_copy_files.py b/imod/tests/test_msw/test_copy_files.py new file mode 100644 index 000000000..80d8ceb08 --- /dev/null +++ b/imod/tests/test_msw/test_copy_files.py @@ -0,0 +1,69 @@ +from pytest_cases import parametrize_with_cases + +from imod.msw.copy_files import FileCopier + + +def write_test_files(directory, filenames): + paths = [directory / filename for filename in filenames] + for p in paths: + with open(p, mode="w") as f: + f.write("test") + return paths + + +def case_simple_files(tmp_path_factory): + directory = tmp_path_factory.mktemp("simple_files") + filenames = [ + "a.inp", + "b.inp", + "c.inp", + ] + return write_test_files(directory, filenames) + + +def case_imod5_extra_files(tmp_path_factory): + directory = tmp_path_factory.mktemp("imod5_extra_files") + filenames = [ + "a.inp", + "b.inp", + "c.inp", + "mete_grid.inp", + "para_sim.inp", + "svat2precgrid.inp", + "svat2etrefgrid.inp", + ] + return write_test_files(directory, filenames) + + +@parametrize_with_cases("src_files", cases=".") +def test_copyfile_init(src_files): + # Act + copyfiles = FileCopier(src_files) + # Arrange + assert "paths" in copyfiles.dataset.keys() + assert len(copyfiles.dataset["paths"]) == len(src_files) + + +@parametrize_with_cases("src_files", cases=".") +def test_copyfile_write(src_files, tmp_path): + # Arrange + expected_filenames = {f.name for f in src_files} + # Act + copyfiles = FileCopier(src_files) + copyfiles.write(tmp_path) + # Assert + actual_filepaths = tmp_path.glob("*.inp") + actual_filenames = {f.name for f in actual_filepaths} + diff = expected_filenames ^ actual_filenames + assert len(diff) == 0 + + +@parametrize_with_cases("src_files", cases=".") +def test_from_imod5_data(src_files): + # Arrange + imod5_ls = [[p] for p in src_files] + imod5_data = {"extra": {"paths": imod5_ls}} + # Act + copyfiles = FileCopier.from_imod5_data(imod5_data) + # Assert + len(copyfiles.dataset["paths"]) == 3 From 8e55dd5db24f06d607efad80a3babf0a4b06afb0 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Thu, 5 Dec 2024 11:50:38 +0100 Subject: [PATCH 08/13] Add from_imod5_data method and test (#1325) Fixes #1324 # Description Adds a ``ScalingFactors.from_imod5_data`` method. I think, the final missing constructor method. Scaling factors for urban areas are set 1.0, not sure if iMOD5 does this as well @HendrikKok? # Checklist - [x] Links to correct issue - [ ] Update changelog, if changes affect users - [x] PR title starts with ``Issue #nr``, e.g. ``Issue #737`` - [x] Unit tests were added - [ ] **If feature added**: Added/extended example --- imod/msw/scaling_factors.py | 26 +++++++++ imod/tests/test_msw/test_scaling_factors.py | 60 +++++++++++++++++---- 2 files changed, 76 insertions(+), 10 deletions(-) diff --git a/imod/msw/scaling_factors.py b/imod/msw/scaling_factors.py index 6e544d1c5..87dccaf61 100644 --- a/imod/msw/scaling_factors.py +++ b/imod/msw/scaling_factors.py @@ -1,7 +1,12 @@ +from typing import cast + from imod.mf6.interfaces.iregridpackage import IRegridPackage from imod.msw.fixed_format import VariableMetaData from imod.msw.pkgbase import MetaSwapPackage from imod.msw.regrid.regrid_schemes import ScalingRegridMethod +from imod.msw.utilities.common import concat_imod5 +from imod.typing import GridDataDict, Imod5DataDict +from imod.typing.grid import ones_like class ScalingFactors(MetaSwapPackage, IRegridPackage): @@ -76,3 +81,24 @@ def __init__( self.dataset["depth_perched_water_table"] = depth_perched_water_table self._pkgcheck() + + @classmethod + def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "ScalingFactors": + """ + Import ScalingFactors from iMOD5 data. Pressure head factor is set to + one for all factors, as well as all factors for urban areas. + """ + cap_data = cast(GridDataDict, imod5_data["cap"]) + grid_ones = ones_like(cap_data["boundary"]) + + data = {} + data["scale_soil_moisture"] = concat_imod5( + cap_data["soil_moisture_fraction"], grid_ones + ) + data["scale_hydraulic_conductivity"] = concat_imod5( + cap_data["conductivitiy_factor"], grid_ones + ) + data["scale_pressure_head"] = concat_imod5(grid_ones, grid_ones) + data["depth_perched_water_table"] = cap_data["perched_water_table_level"] + + return cls(**data) diff --git a/imod/tests/test_msw/test_scaling_factors.py b/imod/tests/test_msw/test_scaling_factors.py index d04f071de..9b9f8c36f 100644 --- a/imod/tests/test_msw/test_scaling_factors.py +++ b/imod/tests/test_msw/test_scaling_factors.py @@ -10,11 +10,12 @@ RegridderWeightsCache, ) from imod.msw import ScalingFactors +from imod.typing.grid import ones_like -def setup_scaling_factor(): +def setup_scaling_factor_grids(): x = [1.0, 2.0, 3.0] - y = [1.0, 2.0, 3.0] + y = [3.0, 2.0, 1.0] subunit = [0, 1] dx = 1.0 dy = 1.0 @@ -63,6 +64,12 @@ def setup_scaling_factor(): # fmt: on index = (svat != 0).values.ravel() + return scale, depth_perched_water_table, index, svat + + +def test_simple_model(fixed_format_parser): + scale, depth_perched_water_table, index, svat = setup_scaling_factor_grids() + scaling_factors = ScalingFactors( scale_soil_moisture=scale, scale_hydraulic_conductivity=scale, @@ -70,12 +77,6 @@ def setup_scaling_factor(): depth_perched_water_table=depth_perched_water_table, ) - return scaling_factors, index, svat - - -def test_simple_model(fixed_format_parser): - scaling_factors, index, svat = setup_scaling_factor() - with tempfile.TemporaryDirectory() as output_dir: output_dir = Path(output_dir) scaling_factors.write(output_dir, index, svat, None, None) @@ -95,8 +96,16 @@ def test_simple_model(fixed_format_parser): ) -def test_regrid_scaling_factor(fixed_format_parser, simple_2d_grid_with_subunits): - scaling_factors, _, _ = setup_scaling_factor() +def test_regrid_scaling_factor(simple_2d_grid_with_subunits): + scale, depth_perched_water_table, _, _ = setup_scaling_factor_grids() + + scaling_factors = ScalingFactors( + scale_soil_moisture=scale, + scale_hydraulic_conductivity=scale, + scale_pressure_head=scale, + depth_perched_water_table=depth_perched_water_table, + ) + new_grid = simple_2d_grid_with_subunits regrid_context = RegridderWeightsCache() @@ -105,3 +114,34 @@ def test_regrid_scaling_factor(fixed_format_parser, simple_2d_grid_with_subunits assert np.all(regridded_scaling_factor.dataset["x"].values == new_grid["x"].values) assert np.all(regridded_scaling_factor.dataset["y"].values == new_grid["y"].values) + + +def test_from_imod5_data(fixed_format_parser): + scale, depth_perched_water_table, index, svat = setup_scaling_factor_grids() + + imod5_data = {"cap": {}} + scale_rural = scale.sel(subunit=0, drop=True) + imod5_data["cap"]["boundary"] = ones_like(scale_rural) + imod5_data["cap"]["soil_moisture_fraction"] = scale_rural + imod5_data["cap"]["conductivitiy_factor"] = scale_rural + imod5_data["cap"]["perched_water_table_level"] = depth_perched_water_table + + scaling_factors = ScalingFactors.from_imod5_data(imod5_data) + + with tempfile.TemporaryDirectory() as output_dir: + output_dir = Path(output_dir) + scaling_factors.write(output_dir, index, svat, None, None) + + results = fixed_format_parser( + output_dir / ScalingFactors._file_name, ScalingFactors._metadata_dict + ) + + assert_equal(results["svat"], np.array([1, 2, 3, 4])) + assert_almost_equal(results["scale_soil_moisture"], np.array([0.5, 1.0, 1.0, 1.0])) + assert_almost_equal( + results["scale_hydraulic_conductivity"], np.array([0.5, 1.0, 1.0, 1.0]) + ) + assert_almost_equal(results["scale_pressure_head"], np.array([1.0, 1.0, 1.0, 1.0])) + assert_almost_equal( + results["depth_perched_water_table"], np.array([0.5, 1.0, 0.5, 0.7]) + ) From bec1bebc4a1187ba0a2478079975ec589afd9994 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Thu, 12 Dec 2024 17:06:14 +0100 Subject: [PATCH 09/13] Issue #1316 msw model from imod5 data (#1328) Fixes #1316 # Description Adds ``MetaSwapModel.from_imod5_data``, currently only models can be written by turning off validation, that is still something to pick up. - Adds ``MetaSwapModel.from_imod5_data`` - Converts ``Imod5DataDict`` TypeAlias to TypedDict, with which we can avoid a lot of type casting. - Add masking and broadcasting utilities for MetaSWAP iMOD5 data. I could not easily make calls to the MODFLOW6 utilities for masking, as these assumed presence of a layer coordinate and had logic to deal with that, which is irrelevant for MetaSWAP. - Adds ``MetaSwapActive`` dataclass to store masks - Allow user-defined ``settings`` upon initialization, instead of only supporting default settings. - Add ``validate`` argument to ``MetaSwapModel.write`` to turn off model validation upon writing. # Checklist - [x] Links to correct issue - [x] Update changelog, if changes affect users - [x] PR title starts with ``Issue #nr``, e.g. ``Issue #737`` - [x] Unit tests were added - [ ] **If feature added**: Added/extended example --- docs/api/changelog.rst | 6 + docs/api/msw.rst | 11 + imod/mf6/model_gwf.py | 8 +- imod/mf6/rch.py | 3 +- imod/mf6/simulation.py | 2 + imod/mf6/wel.py | 2 + imod/msw/grid_data.py | 63 +++-- imod/msw/infiltration.py | 9 +- imod/msw/meteo_grid.py | 9 +- imod/msw/meteo_mapping.py | 4 +- imod/msw/model.py | 102 +++++++- imod/msw/ponding.py | 4 +- imod/msw/scaling_factors.py | 6 +- imod/msw/sprinkling.py | 6 +- imod/msw/utilities/imod5_converter.py | 40 ++- imod/msw/utilities/mask.py | 60 +++++ imod/msw/utilities/parse.py | 38 +++ imod/msw/utilities/select.py | 16 ++ imod/prepare/__init__.py | 2 + imod/tests/conftest.py | 2 + imod/tests/fixtures/msw_imod5_cap_fixture.py | 237 ++++++++++++++++++ imod/tests/fixtures/msw_meteo_fixture.py | 43 ++++ imod/tests/test_msw/test_grid_data.py | 2 +- imod/tests/test_msw/test_infiltration.py | 10 +- imod/tests/test_msw/test_meteo_grid.py | 144 ++++------- imod/tests/test_msw/test_model.py | 157 ++++++++++++ .../test_msw/test_utilities/test_mask.py | 131 ++++++++++ .../test_msw/test_utilities/test_parse.py | 28 +++ imod/typing/__init__.py | 6 +- 29 files changed, 985 insertions(+), 166 deletions(-) create mode 100644 imod/msw/utilities/mask.py create mode 100644 imod/msw/utilities/parse.py create mode 100644 imod/msw/utilities/select.py create mode 100644 imod/tests/fixtures/msw_imod5_cap_fixture.py create mode 100644 imod/tests/fixtures/msw_meteo_fixture.py create mode 100644 imod/tests/test_msw/test_utilities/test_mask.py create mode 100644 imod/tests/test_msw/test_utilities/test_parse.py diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst index 1f073881c..f3927430e 100644 --- a/docs/api/changelog.rst +++ b/docs/api/changelog.rst @@ -21,6 +21,12 @@ Added 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. +- :meth:`imod.msw.MetaSwapModel.from_imod5_data` to construct a MetaSWAP model + from data in an iMOD5 projectfile. +- :meth:`imod.msw.MetaSwapModel.write` has a ``validate`` argument, which can be + used to turn off validation upon writing, use at your own risk! +- :class:`imod.msw.MetaSwapModel` got ``settings`` argument to set simulation + settings. Fixed ~~~~~ diff --git a/docs/api/msw.rst b/docs/api/msw.rst index f2260ed0f..067710445 100644 --- a/docs/api/msw.rst +++ b/docs/api/msw.rst @@ -7,10 +7,15 @@ MetaSWAP :toctree: generated/msw GridData + GridData.from_imod5_data Infiltration + Infiltration.from_imod5_data Ponding + Ponding.from_imod5_data ScalingFactors + ScalingFactors.from_imod5_data Sprinkling + Sprinkling.from_imod5_data IdfMapping TimeOutputControl @@ -25,9 +30,15 @@ MetaSWAP AnnualCropFactors MeteoGrid + MeteoGridCopy + MeteoGridCopy.from_imod5_data EvapotranspirationMapping + EvapotranspirationMapping.from_imod5_data PrecipitationMapping + PrecipitationMapping.from_imod5_data CouplerMapping MetaSwapModel + MetaSwapModel.write + MetaSwapModel.from_imod5_data diff --git a/imod/mf6/model_gwf.py b/imod/mf6/model_gwf.py index e83ee5e1c..ebc878876 100644 --- a/imod/mf6/model_gwf.py +++ b/imod/mf6/model_gwf.py @@ -218,9 +218,11 @@ def from_imod5_data( times: list[datetime] Time discretization of the simulation. These times are used for the following: - * Times of wells with associated timeseries are resampled to these times - * Start- and end times in the list are used to repeat the stresses - of periodic data (e.g. river stages in iMOD5 for "summer", "winter") + + * Times of wells with associated timeseries are resampled to these times + * Start- and end times in the list are used to repeat the stresses + of periodic data (e.g. river stages in iMOD5 for "summer", "winter") + allocation_options: SimulationAllocationOptions object containing the allocation options per package type. If you want a package to have a different allocation option, diff --git a/imod/mf6/rch.py b/imod/mf6/rch.py index 392cc9d56..3fd892c90 100644 --- a/imod/mf6/rch.py +++ b/imod/mf6/rch.py @@ -249,7 +249,8 @@ def from_imod5_cap_data( 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) + msw_active = is_msw_active_cell(target_dis, cap_data, msw_area) + active = msw_active.all data = {} layer_da = xr.full_like( diff --git a/imod/mf6/simulation.py b/imod/mf6/simulation.py index 5bc2469f1..4f324f31b 100644 --- a/imod/mf6/simulation.py +++ b/imod/mf6/simulation.py @@ -1361,11 +1361,13 @@ def from_imod5_data( times: list[datetime] time discretization of the model to be imported. This is used for the following: + * Times of wells with associated timeseries are resampled to these times * Start and end times in the list are used to repeat the stresses of periodic data (e.g. river stages in iMOD5 for "summer", "winter") * The simulation is discretized to these times, using :meth:`imod.mf6.Modflow6Simulation.create_time_discretization` + allocation_options: SimulationAllocationOptions, optional object containing the allocation options per package type. If you want a specific package to have a different allocation option, then diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py index d30b58d30..11334b2cf 100644 --- a/imod/mf6/wel.py +++ b/imod/mf6/wel.py @@ -1312,9 +1312,11 @@ def from_imod5_cap_data(cls, imod5_data: Imod5DataDict): (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. diff --git a/imod/msw/grid_data.py b/imod/msw/grid_data.py index ad7299206..d4c387bd3 100644 --- a/imod/msw/grid_data.py +++ b/imod/msw/grid_data.py @@ -1,14 +1,8 @@ -from typing import Optional - import numpy as np import xarray as xr from imod.mf6 import StructuredDiscretization from imod.mf6.interfaces.iregridpackage import IRegridPackage -from imod.mf6.regrid.regrid_schemes import ( - RegridMethodType, -) -from imod.mf6.utilities.regrid import RegridderWeightsCache from imod.msw.fixed_format import VariableMetaData from imod.msw.pkgbase import MetaSwapPackage from imod.msw.regrid.regrid_schemes import GridDataRegridMethod @@ -18,7 +12,8 @@ get_rootzone_depth_from_imod5_data, is_msw_active_cell, ) -from imod.typing import GridDataDict +from imod.msw.utilities.mask import MetaSwapActive, mask_and_broadcast_pkg_data +from imod.typing import Imod5DataDict from imod.util.spatial import get_cell_area, spatial_reference @@ -129,12 +124,38 @@ def _pkgcheck(self): @classmethod def from_imod5_data( cls, - imod5_data: dict[str, GridDataDict], + imod5_data: Imod5DataDict, target_dis: StructuredDiscretization, - regridder_types: Optional[RegridMethodType] = None, - regrid_cache: RegridderWeightsCache = RegridderWeightsCache(), - ) -> "GridData": - # Get iMOD5 capillary zone data + ) -> tuple["GridData", MetaSwapActive]: + """ + Construct a MetaSWAP GridData package from iMOD5 data in the CAP + package, loaded with the :func:`imod.formats.prj.open_projectfile_data` + function. + + Method does the following things: + + - Computes rural area from the wetted area and urban area grids. + - Sets a second subunit for urban landuse (== 18). + - Rootzone depth is converted from cm's to m's. + - Determines where MetaSWAP should be active based on area grids, + boundary grid, and MODFLOW6 idomain. + + Parameters + ---------- + imod5_data: Imod5DataDict + iMOD5 data as returned by + :func:`imod.formats.prj.open_projectfile_data` + target_dis: imod.mf6.StructuredDiscretization + MODFLOW6 discretization to couple the MetaSWAP model to. + + Returns + ------- + imod.msw.GridData + GridData package + MetaSwapActive + Dataclass containing where MetaSwAP is active, per subunit, as well + as aggregated over subunits. + """ imod5_cap = imod5_data["cap"] data = {} @@ -142,17 +163,9 @@ def from_imod5_data( data["landuse"] = get_landuse_from_imod5_data(imod5_cap) data["rootzone_depth"] = get_rootzone_depth_from_imod5_data(imod5_cap) data["surface_elevation"] = imod5_cap["surface_elevation"] - data["soil_physical_unit"] = imod5_cap["soil_physical_unit"] - - active, subunit_active = is_msw_active_cell(target_dis, imod5_cap, data["area"]) + data["soil_physical_unit"] = imod5_cap["soil_physical_unit"].astype(int) - data_active = { - key: ( - griddata.where(subunit_active) - if key in cls._with_subunit - else griddata.where(active) - ) - for key, griddata in data.items() - } - data_active["active"] = active - return cls(**data_active) + msw_active = is_msw_active_cell(target_dis, imod5_cap, data["area"]) + data_active = mask_and_broadcast_pkg_data(cls, data, msw_active) + data_active["active"] = msw_active.all + return cls(**data_active), msw_active diff --git a/imod/msw/infiltration.py b/imod/msw/infiltration.py index 830365681..7e5472c7b 100644 --- a/imod/msw/infiltration.py +++ b/imod/msw/infiltration.py @@ -8,7 +8,8 @@ from imod.msw.pkgbase import MetaSwapPackage from imod.msw.regrid.regrid_schemes import InfiltrationRegridMethod from imod.msw.utilities.common import concat_imod5 -from imod.typing import GridDataDict +from imod.msw.utilities.mask import MaskValues +from imod.typing import GridDataDict, Imod5DataDict from imod.typing.grid import ones_like @@ -28,7 +29,7 @@ def deactivate_small_resistances_in_data(data: GridDataDict): message=message.format(var=var), additional_depth=1, ) - data[var] = data[var].where(~to_deactivate, -9999.0) + data[var] = data[var].where(~to_deactivate, MaskValues.default) return data @@ -101,7 +102,7 @@ def __init__( self._pkgcheck() @classmethod - def from_imod5_data(cls, imod5_data: dict[str, GridDataDict]) -> "Infiltration": + def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "Infiltration": cap_data = imod5_data["cap"] data = {} # Use runon resistance as downward resistance, and runoff for downward @@ -120,7 +121,7 @@ def from_imod5_data(cls, imod5_data: dict[str, GridDataDict]) -> "Infiltration": data = deactivate_small_resistances_in_data(data) like = data["downward_resistance"].isel(subunit=0, drop=True) - data["bottom_resistance"] = ones_like(like) + data["bottom_resistance"] = ones_like(like) * MaskValues.default data["extra_storage_coefficient"] = ones_like(like) return cls(**data) diff --git a/imod/msw/meteo_grid.py b/imod/msw/meteo_grid.py index 70d3625a5..3974eb838 100644 --- a/imod/msw/meteo_grid.py +++ b/imod/msw/meteo_grid.py @@ -1,7 +1,7 @@ import csv from pathlib import Path from shutil import copyfile -from typing import Optional, Union, cast +from typing import Optional, Union import numpy as np import pandas as pd @@ -13,6 +13,7 @@ from imod.msw.regrid.regrid_schemes import MeteoGridRegridMethod from imod.msw.timeutil import to_metaswap_timeformat from imod.msw.utilities.common import find_in_file_list +from imod.msw.utilities.mask import MaskValues from imod.typing import Imod5DataDict from imod.util.regrid_method_type import EmptyRegridMethod, RegridMethodType @@ -176,7 +177,9 @@ def write(self, directory: Union[str, Path], *args): path = (directory / self._meteo_dirname / str(varname)).with_suffix( ".asc" ) - imod.rasterio.save(path, self.dataset[str(varname)], nodata=-9999.0) + imod.rasterio.save( + path, self.dataset[str(varname)], nodata=MaskValues.default + ) def _pkgcheck(self): for varname in self.dataset.data_vars: @@ -223,7 +226,7 @@ def write(self, directory: Path | str, *args): @classmethod def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "MeteoGridCopy": - paths = cast(list[str], imod5_data["extra"]["paths"]) + paths = imod5_data["extra"]["paths"] filepath = find_in_file_list(cls._file_name, paths) return cls(filepath) diff --git a/imod/msw/meteo_mapping.py b/imod/msw/meteo_mapping.py index 1de00fa69..b728bb497 100644 --- a/imod/msw/meteo_mapping.py +++ b/imod/msw/meteo_mapping.py @@ -1,7 +1,7 @@ from copy import deepcopy from pathlib import Path from textwrap import dedent -from typing import Any, Optional, TextIO, cast +from typing import Any, Optional, TextIO import numpy as np import pandas as pd @@ -58,7 +58,7 @@ def open_first_meteo_grid(mete_grid_path: str | Path, column_nr: int) -> xr.Data def open_first_meteo_grid_from_imod5_data(imod5_data: Imod5DataDict, column_nr: int): - paths = cast(list[str], imod5_data["extra"]["paths"]) + paths = imod5_data["extra"]["paths"] metegrid_path = find_in_file_list("mete_grid.inp", paths) return open_first_meteo_grid(metegrid_path, column_nr=column_nr) diff --git a/imod/msw/model.py b/imod/msw/model.py index f349021bb..6075ca59a 100644 --- a/imod/msw/model.py +++ b/imod/msw/model.py @@ -1,14 +1,17 @@ import collections from copy import copy +from datetime import datetime from pathlib import Path -from typing import Optional, Tuple, Union +from typing import Any, Optional, Tuple, Union, cast import jinja2 import numpy as np +import xarray as xr from imod.mf6.dis import StructuredDiscretization from imod.mf6.mf6_wel_adapter import Mf6Wel from imod.mf6.utilities.regrid import RegridderWeightsCache +from imod.msw.copy_files import FileCopier from imod.msw.coupler_mapping import CouplerMapping from imod.msw.grid_data import GridData from imod.msw.idf_mapping import IdfMapping @@ -20,11 +23,20 @@ InitialConditionsSavedState, ) from imod.msw.landuse import LanduseOptions -from imod.msw.meteo_grid import MeteoGrid +from imod.msw.meteo_grid import MeteoGrid, MeteoGridCopy from imod.msw.meteo_mapping import EvapotranspirationMapping, PrecipitationMapping from imod.msw.output_control import TimeOutputControl +from imod.msw.ponding import Ponding +from imod.msw.scaling_factors import ScalingFactors +from imod.msw.sprinkling import Sprinkling from imod.msw.timeutil import to_metaswap_timeformat +from imod.msw.utilities.common import find_in_file_list +from imod.msw.utilities.imod5_converter import has_active_scaling_factor +from imod.msw.utilities.mask import MaskValues, mask_and_broadcast_cap_data +from imod.msw.utilities.parse import read_para_sim +from imod.msw.utilities.select import drop_layer_dim_cap_data from imod.msw.vegetation import AnnualCropFactors +from imod.typing import Imod5DataDict from imod.util.regrid_method_type import RegridderType REQUIRED_PACKAGES = ( @@ -88,6 +100,7 @@ class MetaSwapModel(Model): ---------- unsaturated_database: Path-like or str Path to the MetaSWAP soil physical database folder. + settings: dict """ _pkg_id = "model" @@ -99,10 +112,18 @@ class MetaSwapModel(Model): "{%endfor%}" ) - def __init__(self, unsaturated_database): + def __init__( + self, + unsaturated_database: Path | str, + settings: Optional[dict[str, Any]] = None, + ): super().__init__() - self.simulation_settings = copy(DEFAULT_SETTINGS) + if settings is None: + self.simulation_settings = copy(DEFAULT_SETTINGS) + else: + self.simulation_settings = settings + self.simulation_settings["unsa_svat_path"] = ( self._render_unsaturated_database_path(unsaturated_database) ) @@ -210,6 +231,7 @@ def write( directory: Union[str, Path], mf6_dis: StructuredDiscretization, mf6_wel: Mf6Wel, + validate: bool = True, ): """ Write packages and simulation settings (para_sim.inp). @@ -221,9 +243,10 @@ def write( """ # Model checks - self._check_required_packages() - self._check_vegetation_indices_in_annual_crop_factors() - self._check_landuse_indices_in_lookup_options() + if validate: + self._check_required_packages() + self._check_vegetation_indices_in_annual_crop_factors() + self._check_landuse_indices_in_lookup_options() # Force to Path directory = Path(directory) @@ -258,7 +281,7 @@ def regrid_like( regrid_context: Optional[RegridderWeightsCache] = None, regridder_types: Optional[dict[str, Tuple[RegridderType, str]]] = None, ) -> "MetaSwapModel": - unsat_database = self.simulation_settings["unsa_svat_path"] + unsat_database = cast(str, self.simulation_settings["unsa_svat_path"]) regridded_model = MetaSwapModel(unsat_database) target_grid = mf6_regridded_dis["idomain"] @@ -282,3 +305,66 @@ def regrid_like( regridded_model[mod2svat_name] = CouplerMapping() return regridded_model + + @classmethod + def from_imod5_data( + cls, + imod5_data: Imod5DataDict, + target_dis: StructuredDiscretization, + times: list[datetime], + ): + """ + Construct a MetaSWAP model from iMOD5 data in the CAP package, loaded + with the :func:`imod.formats.prj.open_projectfile_data` function. + + Parameters + ---------- + imod5_data: dict + Dictionary with iMOD5 data. This can be constructed from the + :func:`imod.formats.prj.open_projectfile_data` method. + target_dis: imod.mf6.StructuredDiscretization + Target discretization, cells where MODLOW6 is inactive will be + inactive in MetaSWAP as well. + times: list[datetime] + List of datetimes, will be used to set the output control times. + Is also used to infer the starttime of the simulation. + + Returns + ------- + MetaSwapModel + MetaSWAP model imported from imod5 data. + """ + extra_paths = imod5_data["extra"]["paths"] + path_to_parasim = find_in_file_list("para_sim.inp", extra_paths) + parasim_settings = read_para_sim(path_to_parasim) + unsa_svat_path = cast(str, parasim_settings["unsa_svat_path"]) + # Drop layer coord + imod5_cap_no_layer = drop_layer_dim_cap_data(imod5_data) + model = cls(unsa_svat_path, parasim_settings) + model["grid"], msw_active = GridData.from_imod5_data( + imod5_cap_no_layer, target_dis + ) + cap_data_masked = mask_and_broadcast_cap_data( + imod5_cap_no_layer["cap"], msw_active + ) + imod5_masked: Imod5DataDict = { + "cap": cap_data_masked, + "extra": {"paths": extra_paths}, + } + model["infiltration"] = Infiltration.from_imod5_data(imod5_masked) + model["ponding"] = Ponding.from_imod5_data(imod5_masked) + model["sprinkling"] = Sprinkling.from_imod5_data(imod5_masked) + model["meteo_grid"] = MeteoGridCopy.from_imod5_data(imod5_masked) + model["prec_mapping"] = PrecipitationMapping.from_imod5_data(imod5_masked) + model["evt_mapping"] = EvapotranspirationMapping.from_imod5_data(imod5_masked) + if has_active_scaling_factor(imod5_cap_no_layer["cap"]): + model["scaling_factor"] = ScalingFactors.from_imod5_data(imod5_masked) + area = model["grid"]["area"].isel(subunit=0, drop=True) + model["idf_mapping"] = IdfMapping(area, MaskValues.default) + model["coupling"] = CouplerMapping() + model["extra_files"] = FileCopier.from_imod5_data(imod5_masked) + + times_da = xr.DataArray(times, coords={"time": times}, dims=("time",)) + model["time_oc"] = TimeOutputControl(times_da) + + return model diff --git a/imod/msw/ponding.py b/imod/msw/ponding.py index 89ab6af44..86a7c5cff 100644 --- a/imod/msw/ponding.py +++ b/imod/msw/ponding.py @@ -8,7 +8,7 @@ from imod.msw.pkgbase import DataDictType, MetaSwapPackage from imod.msw.regrid.regrid_schemes import PondingRegridMethod from imod.msw.utilities.common import concat_imod5 -from imod.typing import GridDataDict, IntArray +from imod.typing import Imod5DataDict, IntArray class Ponding(MetaSwapPackage, IRegridPackage): @@ -73,7 +73,7 @@ def _render(self, file: TextIO, index: IntArray, svat: xr.DataArray, *args: Any) return self.write_dataframe_fixed_width(file, dataframe) @classmethod - def from_imod5_data(cls, imod5_data: dict[str, GridDataDict]) -> "Ponding": + def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "Ponding": """ Concatenate ponding depths along subunits """ diff --git a/imod/msw/scaling_factors.py b/imod/msw/scaling_factors.py index 87dccaf61..91984dd3d 100644 --- a/imod/msw/scaling_factors.py +++ b/imod/msw/scaling_factors.py @@ -1,11 +1,9 @@ -from typing import cast - from imod.mf6.interfaces.iregridpackage import IRegridPackage from imod.msw.fixed_format import VariableMetaData from imod.msw.pkgbase import MetaSwapPackage from imod.msw.regrid.regrid_schemes import ScalingRegridMethod from imod.msw.utilities.common import concat_imod5 -from imod.typing import GridDataDict, Imod5DataDict +from imod.typing import Imod5DataDict from imod.typing.grid import ones_like @@ -88,7 +86,7 @@ def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "ScalingFactors": Import ScalingFactors from iMOD5 data. Pressure head factor is set to one for all factors, as well as all factors for urban areas. """ - cap_data = cast(GridDataDict, imod5_data["cap"]) + cap_data = imod5_data["cap"] grid_ones = ones_like(cap_data["boundary"]) data = {} diff --git a/imod/msw/sprinkling.py b/imod/msw/sprinkling.py index ea9925112..b5fa8cc69 100644 --- a/imod/msw/sprinkling.py +++ b/imod/msw/sprinkling.py @@ -1,4 +1,4 @@ -from typing import TextIO, cast +from typing import TextIO import numpy as np import pandas as pd @@ -180,9 +180,11 @@ def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "Sprinkling": (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. @@ -208,7 +210,7 @@ def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "Sprinkling": ------- Sprinkling package """ - cap_data = cast(GridDataDict, imod5_data["cap"]) + cap_data = imod5_data["cap"] if isinstance(cap_data["artificial_recharge_layer"], pd.DataFrame): data = _sprinkling_data_from_imod5_ipf(cap_data) else: diff --git a/imod/msw/utilities/imod5_converter.py b/imod/msw/utilities/imod5_converter.py index 70f88c982..0c4d21604 100644 --- a/imod/msw/utilities/imod5_converter.py +++ b/imod/msw/utilities/imod5_converter.py @@ -1,6 +1,9 @@ +from xarray.core.utils import is_scalar + from imod.logging import LogLevel, logger from imod.mf6 import StructuredDiscretization from imod.msw.utilities.common import concat_imod5 +from imod.msw.utilities.mask import MaskValues, MetaSwapActive from imod.typing import GridDataArray, GridDataDict from imod.typing.grid import ones_like from imod.util.spatial import get_cell_area @@ -10,7 +13,7 @@ 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"]) + mf6_area = get_cell_area(imod5_cap["boundary"]) wetted_area = imod5_cap["wetted_area"] urban_area = imod5_cap["urban_area"] rural_area = mf6_area - (wetted_area + urban_area) @@ -43,7 +46,7 @@ def get_landuse_from_imod5_data( rural_landuse = imod5_cap["landuse"] # Urban landuse = 18 urban_landuse = ones_like(rural_landuse) * 18 - return concat_imod5(rural_landuse, urban_landuse) + return concat_imod5(rural_landuse, urban_landuse).astype(int) def get_rootzone_depth_from_imod5_data( @@ -63,7 +66,7 @@ def is_msw_active_cell( target_dis: StructuredDiscretization, imod5_cap: GridDataDict, msw_area: GridDataArray, -) -> tuple[GridDataArray, GridDataArray]: +) -> MetaSwapActive: """ Return grid of cells that are active in the coupled computation, based on following criteria: @@ -84,4 +87,33 @@ def is_msw_active_cell( (imod5_cap["boundary"] == 1) & (msw_area > 0) & (mf6_top_active >= 1) ) active = subunit_active.any(dim="subunit") - return active, subunit_active + return MetaSwapActive(active, subunit_active) + + +def _is_equal_scalar_value(da, value): + """ + Helper function to guarantee that the check in ``has_active_scaling_factor`` + can shortcut after is_scalar returns False. + """ + return da.to_numpy()[()] == value + + +def has_active_scaling_factor(imod5_cap: GridDataDict): + """ + Check if scaling factor grids are active. Carefully checks if data is + provided as constant (scalar) and if it matches an inactivity value. The + function shortcuts if data is provided as constant. + """ + variable_inactive_mapping = { + "perched_water_table_level": MaskValues.default, + "soil_moisture_fraction": 1.0, + "conductivitiy_factor": 1.0, + } + scaling_factor_inactive = True + for var, inactive_value in variable_inactive_mapping.items(): + da = imod5_cap[var] + scaling_factor_inactive &= is_scalar(da) and _is_equal_scalar_value( + da, inactive_value + ) + + return not scaling_factor_inactive diff --git a/imod/msw/utilities/mask.py b/imod/msw/utilities/mask.py new file mode 100644 index 000000000..6a26c9b9d --- /dev/null +++ b/imod/msw/utilities/mask.py @@ -0,0 +1,60 @@ +import numbers +from dataclasses import dataclass + +from imod.msw.pkgbase import MetaSwapPackage +from imod.typing import GridDataArray, GridDataDict + + +@dataclass +class MetaSwapActive: + all: GridDataArray + per_subunit: GridDataArray + + +@dataclass +class MaskValues: + """Stores sentinel values for nodata. Most cases use -9999.0, but exceptions + can be added here.""" + + default = -9999.0 + integer = 0 + + +def mask_and_broadcast_cap_data( + cap_data: GridDataDict, msw_active: MetaSwapActive +) -> GridDataDict: + """ + Mask and broadcast cap data, always mask with "all" of MetaSwapActive. + """ + return { + key: _mask_spatial_var(grid, msw_active.all) for key, grid in cap_data.items() + } + + +def mask_and_broadcast_pkg_data( + package: type[MetaSwapPackage], grid_data: GridDataDict, msw_active: MetaSwapActive +) -> GridDataDict: + """ + Mask and broadcast grid data, carefully mask per subunit if variable needs + to contain subunits. + """ + + return { + key: ( + _mask_spatial_var(grid, msw_active.per_subunit) + if key in package._with_subunit + else _mask_spatial_var(grid, msw_active.all) + ) + for key, grid in grid_data.items() + } + + +def _mask_spatial_var(da: GridDataArray, active: GridDataArray) -> GridDataArray: + if issubclass(da.dtype.type, numbers.Integral): + return da.where(active, other=MaskValues.integer) + elif issubclass(da.dtype.type, numbers.Real): + return da.where(active) + else: + raise TypeError( + f"Expected dtype float or integer. Received instead: {da.dtype}" + ) diff --git a/imod/msw/utilities/parse.py b/imod/msw/utilities/parse.py new file mode 100644 index 000000000..efe7f4428 --- /dev/null +++ b/imod/msw/utilities/parse.py @@ -0,0 +1,38 @@ +from pathlib import Path +from typing import TypeAlias + +ScalarType: TypeAlias = str | float | int + + +def _try_parsing_string_to_number(s: str) -> ScalarType: + """ + Convert string to number: + + "1" -> 1 + "1.0" -> 1.0 + "a" -> "a" + """ + try: + return int(s) + except ValueError: + pass + try: + return float(s) + except ValueError: + pass + return s + + +def read_para_sim(file: Path | str) -> dict[str, ScalarType]: + with open(file, "r") as f: + lines = f.readlines() + out = {} + for line in lines: + # Strip comments starting with "!" and split keys from values at the + # equals sign. + key_values = line[0 : line.find("!")].split("=") + if len(key_values) > 1: + key = key_values[0].strip() + value = _try_parsing_string_to_number(key_values[1].strip()) + out[key] = value + return out diff --git a/imod/msw/utilities/select.py b/imod/msw/utilities/select.py new file mode 100644 index 000000000..abd28818f --- /dev/null +++ b/imod/msw/utilities/select.py @@ -0,0 +1,16 @@ +from imod.typing import Imod5DataDict, SelSettingsType + +_DROP_LAYER_KWARGS: SelSettingsType = { + "layer": 0, + "drop": True, + "missing_dims": "ignore", +} + + +def drop_layer_dim_cap_data(imod5_data: Imod5DataDict) -> Imod5DataDict: + cap_data = imod5_data["cap"] + return { + "cap": { + key: da.isel(**_DROP_LAYER_KWARGS).compute() for key, da in cap_data.items() + } + } diff --git a/imod/prepare/__init__.py b/imod/prepare/__init__.py index 03342c87e..41e16f46d 100644 --- a/imod/prepare/__init__.py +++ b/imod/prepare/__init__.py @@ -43,6 +43,8 @@ from imod.prepare.topsystem import ( ALLOCATION_OPTION, DISTRIBUTING_OPTION, + SimulationAllocationOptions, + SimulationDistributingOptions, allocate_drn_cells, allocate_ghb_cells, allocate_rch_cells, diff --git a/imod/tests/conftest.py b/imod/tests/conftest.py index 0677d4b4f..b0f400192 100644 --- a/imod/tests/conftest.py +++ b/imod/tests/conftest.py @@ -92,4 +92,6 @@ well_test_data_transient, ) from .fixtures.msw_fixture import fixed_format_parser, simple_2d_grid_with_subunits +from .fixtures.msw_imod5_cap_fixture import imod5_cap_data +from .fixtures.msw_meteo_fixture import meteo_grids from .fixtures.msw_model_fixture import coupled_mf6_model, coupled_mf6wel, msw_model diff --git a/imod/tests/fixtures/msw_imod5_cap_fixture.py b/imod/tests/fixtures/msw_imod5_cap_fixture.py new file mode 100644 index 000000000..c57925015 --- /dev/null +++ b/imod/tests/fixtures/msw_imod5_cap_fixture.py @@ -0,0 +1,237 @@ +import numpy as np +import pytest +import xarray as xr +from numpy import nan + +from imod.typing import GridDataDict + + +@pytest.fixture(scope="function") +def imod5_cap_data() -> GridDataDict: + x = [1.0, 2.0, 3.0] + y = [3.0, 2.0, 1.0] + layer = [1] + dx = 1.0 + dy = -1.0 + + da_kwargs = {} + da_kwargs["dims"] = ("layer", "y", "x") + da_kwargs["coords"] = {"layer": layer, "y": y, "x": x, "dx": dx, "dy": dy} + + imod5_data = {} + d = {} + + # fmt: off + d["boundary"] = xr.DataArray( + np.array([ + [ + [1, 1, 1], + [0, 0, 0], + [1, 1, 0], + ]], + dtype=int), + **da_kwargs + ) + d["landuse"] = xr.DataArray( + np.array([ + [ + [1, 2, 3], + [0, 0, 0], + [3, 2, 1], + ]], + dtype=int), + **da_kwargs + ) + d["rootzone_thickness"] = xr.DataArray( + np.array([ + [ + [0.1, 0.1, 0.1], + [nan, nan, nan], + [0.2, 0.2, 0.2], + ]] + ), + **da_kwargs + ) + d["soil_physical_unit"] = xr.DataArray( + np.array([ + [ + [1, 1, 1], + [0, 0, 0], + [2, 1, 1], + ]], + dtype=int), + **da_kwargs + ) + d["surface_elevation"] = xr.DataArray( + np.array([ + [ + [1.1, 1.2, 1.3], + [nan, nan, nan], + [1.4, 1.5, 1.6], + ]] + ), + **da_kwargs + ) + d["artificial_recharge"] = xr.DataArray( + np.array([ + [ + [0.2, 0.2, 0.2], + [nan, nan, nan], + [0.3, 0.3, 0.3], + ]] + ), + **da_kwargs + ) + d["artificial_recharge_layer"] = xr.DataArray( + np.array([ + [ + [1, 2, 3], + [0, 0, 0], + [3, 2, 1], + ]], + dtype=int), + **da_kwargs + ) + d["artificial_recharge_capacity"] = xr.DataArray( + np.array([ + [ + [0.4, 0.4, 0.4], + [nan, nan, nan], + [0.6, 0.6, 0.6], + ]] + ), + **da_kwargs + ) + d["wetted_area"] = xr.DataArray( + np.array([ + [ + [0.1, 0.1, 0.1], + [nan, nan, nan], + [0.2, 0.2, 0.2], + ]] + ), + **da_kwargs + ) + d["urban_area"] = xr.DataArray( + np.array([ + [ + [0.2, 0.2, 0.2], + [nan, nan, nan], + [0.3, 0.3, 0.3], + ]] + ), + **da_kwargs + ) + d["urban_ponding_depth"] = xr.DataArray( + np.array([ + [ + [2.2, 2.2, 2.2], + [nan, nan, nan], + [2.3, 2.3, 2.3], + ]] + ), + **da_kwargs + ) + d["rural_ponding_depth"] = xr.DataArray( + np.array([ + [ + [1.2, 1.2, 1.2], + [nan, nan, nan], + [1.3, 1.3, 1.3], + ]] + ), + **da_kwargs + ) + d["urban_runoff_resistance"] = xr.DataArray( + np.array([ + [ + [3.2, 3.2, 3.2], + [nan, nan, nan], + [3.3, 3.3, 3.3], + ]] + ), + **da_kwargs + ) + d["rural_runoff_resistance"] = xr.DataArray( + np.array([ + [ + [3.6, 3.6, 3.6], + [nan, nan, nan], + [3.7, 3.7, 3.7], + ]] + ), + **da_kwargs + ) + d["urban_runon_resistance"] = xr.DataArray( + np.array([ + [ + [5.2, 5.2, 5.2], + [nan, nan, nan], + [5.3, 5.3, 5.3], + ]] + ), + **da_kwargs + ) + d["rural_runon_resistance"] = xr.DataArray( + np.array([ + [ + [5.6, 5.6, 5.6], + [nan, nan, nan], + [5.7, 5.7, 5.7], + ]] + ), + **da_kwargs + ) + d["urban_infiltration_capacity"] = xr.DataArray( + np.array([ + [ + [10.2, 10.2, 10.2], + [nan, nan, nan], + [10.3, 10.3, 10.3], + ]] + ), + **da_kwargs + ) + d["rural_infiltration_capacity"] = xr.DataArray( + np.array([ + [ + [20.2, 20.2, 20.2], + [nan, nan, nan], + [20.3, 20.3, 20.3], + ]] + ), + **da_kwargs + ) + d["perched_water_table_level"]= xr.DataArray( + np.array([ + [ + [2.0, 2.0, 2.0], + [nan, nan, nan], + [2.0, 2.0, 2.0], + ]] + ), + **da_kwargs + ) + d["soil_moisture_fraction"]= xr.DataArray( + np.array([ + [ + [1.5, 1.5, 1.5], + [nan, nan, nan], + [1.5, 1.5, 1.5], + ]] + ), + **da_kwargs + ) + d["conductivitiy_factor"]= xr.DataArray( + np.array([ + [ + [2.5, 2.5, 2.5], + [nan, nan, nan], + [2.5, 2.5, 2.5], + ]] + ), + **da_kwargs + ) + # fmt: on + imod5_data["cap"] = d + return imod5_data diff --git a/imod/tests/fixtures/msw_meteo_fixture.py b/imod/tests/fixtures/msw_meteo_fixture.py new file mode 100644 index 000000000..ebfdee449 --- /dev/null +++ b/imod/tests/fixtures/msw_meteo_fixture.py @@ -0,0 +1,43 @@ +import numpy as np +import pandas as pd +import pytest +import xarray as xr +from numpy import nan + +from imod.typing import GridDataArray + + +@pytest.fixture(scope="function") +def meteo_grids() -> tuple[GridDataArray, GridDataArray]: + x = [1.0, 2.0, 3.0] + y = [1.0, 2.0, 3.0] + time = pd.date_range(start="2000-01-01", end="2000-01-02", freq="D") + + dx = 1.0 + dy = -1.0 + # fmt: off + precipitation = xr.DataArray( + np.array( + [ + [[1.0, 1.0, 1.0], + [nan, nan, nan], + [1.0, 1.0, 1.0]], + + [[2.0, 2.0, 1.0], + [nan, nan, nan], + [1.0, 2.0, 1.0]], + ] + ), + dims=("time", "y", "x"), + coords={"time": time, "y": y, "x": x, "dx": dx, "dy": dy} + ) + + evapotranspiration = xr.DataArray( + np.array( + [1.0, 3.0] + ), + dims=("time",), + coords={"time": time} + ) + # fmt: on + return precipitation, evapotranspiration diff --git a/imod/tests/test_msw/test_grid_data.py b/imod/tests/test_msw/test_grid_data.py index 0d0fbf0c3..a775a93a9 100644 --- a/imod/tests/test_msw/test_grid_data.py +++ b/imod/tests/test_msw/test_grid_data.py @@ -461,7 +461,7 @@ def test_from_imod5_data(grid_data_dict: dict[str, xr.DataArray]): # Dis only needed for idomain dis = StructuredDiscretization(top, top - 0.1, idomain, validate=False) - griddata = GridData.from_imod5_data(imod5_data, target_dis=dis) + griddata, _ = GridData.from_imod5_data(imod5_data, target_dis=dis) expected_rootzone_depth = cap_data["rootzone_thickness"] * 0.01 xr.testing.assert_allclose( expected_rootzone_depth, griddata["rootzone_depth"].sel(subunit=0, drop=True) diff --git a/imod/tests/test_msw/test_infiltration.py b/imod/tests/test_msw/test_infiltration.py index b87940a88..47990e3a6 100644 --- a/imod/tests/test_msw/test_infiltration.py +++ b/imod/tests/test_msw/test_infiltration.py @@ -283,11 +283,11 @@ def test_from_imod5_data(data_infiltration): imod5_data = {"cap": cap_data} actual_pkg = Infiltration.from_imod5_data(imod5_data) - ones_vars = ["bottom_resistance", "extra_storage_coefficient"] - expected = expected_pkg.dataset.drop_vars(ones_vars) - actual = actual_pkg.dataset.drop_vars(ones_vars) + hardcoded_vars = {"bottom_resistance": -9999.0, "extra_storage_coefficient": 1.0} + expected = expected_pkg.dataset.drop_vars(hardcoded_vars.keys()) + actual = actual_pkg.dataset.drop_vars(hardcoded_vars.keys()) xr.testing.assert_equal(actual, expected) - for var in ones_vars: - assert (actual_pkg.dataset[var] == 1.0).all() + for var, value in hardcoded_vars.items(): + assert (actual_pkg.dataset[var] == value).all() diff --git a/imod/tests/test_msw/test_meteo_grid.py b/imod/tests/test_msw/test_meteo_grid.py index b5108842e..66d7fe95e 100644 --- a/imod/tests/test_msw/test_meteo_grid.py +++ b/imod/tests/test_msw/test_meteo_grid.py @@ -7,54 +7,21 @@ import numpy as np import pandas as pd import pytest -import xarray as xr -from numpy import nan from numpy.testing import assert_equal from imod.mf6.utilities.regrid import RegridderWeightsCache from imod.msw import MeteoGrid, MeteoGridCopy -def setup_meteo_grid(): - x = [1.0, 2.0, 3.0] - y = [1.0, 2.0, 3.0] - time = pd.date_range(start="2000-01-01", end="2000-01-02", freq="D") - - dx = 1.0 - dy = -1.0 - # fmt: off - precipitation = xr.DataArray( - np.array( - [ - [[1.0, 1.0, 1.0], - [nan, nan, nan], - [1.0, 1.0, 1.0]], - - [[2.0, 2.0, 1.0], - [nan, nan, nan], - [1.0, 2.0, 1.0]], - ] - ), - dims=("time", "y", "x"), - coords={"time": time, "y": y, "x": x, "dx": dx, "dy": dy} - ) - - evapotranspiration = xr.DataArray( - np.array( - [1.0, 3.0] - ), - dims=("time",), - coords={"time": time} - ) - # fmt: on - - meteo_grid = MeteoGrid(precipitation, evapotranspiration) +def test_meteo_grid_init(meteo_grids): + meteo_grids = MeteoGrid(*meteo_grids) - return meteo_grid + assert meteo_grids.dataset["precipitation"].dims == ("time", "y", "x") + assert meteo_grids.dataset["evapotranspiration"].dims == ("time",) -def test_meteo_grid(): - meteo_grid = setup_meteo_grid() +def test_meteo_grid_write(meteo_grids): + meteo_grid = MeteoGrid(*meteo_grids) with tempfile.TemporaryDirectory() as output_dir: output_dir = Path(output_dir) @@ -83,29 +50,8 @@ def test_meteo_grid(): assert gridnames == expected_filenames -def test_meteo_no_time_grid(): - x = [1.0, 2.0, 3.0] - y = [1.0, 2.0, 3.0] - time = pd.date_range(start="2000-01-01", end="2000-01-02", freq="D") - - dx = 1.0 - dy = -1.0 - # fmt: off - precipitation = xr.DataArray( - np.array( - [ - [[1.0, 1.0, 1.0], - [nan, nan, nan], - [1.0, 1.0, 1.0]], - - [[2.0, 2.0, 1.0], - [nan, nan, nan], - [1.0, 2.0, 1.0]], - ] - ), - dims=("time", "y", "x"), - coords={"time": time, "y": y, "x": x, "dx": dx, "dy": dy} - ) +def test_meteo_no_time_grid(meteo_grids): + precipitation, _ = meteo_grids evapotranspiration = 3.0 # fmt: on @@ -114,8 +60,8 @@ def test_meteo_no_time_grid(): MeteoGrid(precipitation, evapotranspiration) -def test_regrid_meteo(simple_2d_grid_with_subunits): - meteo = setup_meteo_grid() +def test_regrid_meteo(meteo_grids, simple_2d_grid_with_subunits): + meteo = MeteoGrid(*meteo_grids) new_grid = simple_2d_grid_with_subunits regrid_context = RegridderWeightsCache() @@ -126,44 +72,40 @@ def test_regrid_meteo(simple_2d_grid_with_subunits): assert np.all(regridded_ponding.dataset["y"].values == new_grid["y"].values) -def test_meteogridcopy_write(): - # Arrange - meteo_grid = setup_meteo_grid() - - with tempfile.TemporaryDirectory() as output_dir: - grid_dir = Path(output_dir) / "grid" - grid_dir.mkdir(exist_ok=True, parents=True) - meteo_grid.write(grid_dir) +def setup_written_meteo_grids(meteo_grids, tmp_path) -> Path: + meteo_grid = MeteoGrid(*meteo_grids) + grid_dir = tmp_path / "grid" + grid_dir.mkdir(exist_ok=True, parents=True) + meteo_grid.write(grid_dir) + return grid_dir - meteo_grid_copy = MeteoGridCopy(grid_dir / "mete_grid.inp") - copy_dir = Path(output_dir) / "copied" - copy_dir.mkdir(exist_ok=True, parents=True) - # Act - meteo_grid_copy.write(copy_dir) - # Assert - assert filecmp.cmp(grid_dir / "mete_grid.inp", copy_dir / "mete_grid.inp") +def test_meteogridcopy_write(meteo_grids, tmp_path): + # Arrange + grid_dir = setup_written_meteo_grids(meteo_grids, tmp_path) + meteo_grid_copy = MeteoGridCopy(grid_dir / "mete_grid.inp") + copy_dir = tmp_path / "copied" + copy_dir.mkdir(exist_ok=True, parents=True) + # Act + meteo_grid_copy.write(copy_dir) + # Assert + assert filecmp.cmp(grid_dir / "mete_grid.inp", copy_dir / "mete_grid.inp") -def test_meteogridcopy_from_imod5(): - meteo_grid = setup_meteo_grid() - with tempfile.TemporaryDirectory() as output_dir: - grid_dir = Path(output_dir) / "grid" - grid_dir.mkdir(exist_ok=True, parents=True) - meteo_grid.write(grid_dir) - - imod5_data = {} - imod5_data["extra"] = {} - imod5_data["extra"]["paths"] = [ - ["foo"], - [(grid_dir / "mete_grid.inp").resolve()], - ["bar"], - ] - - meteo_grid_copy = MeteoGridCopy.from_imod5_data(imod5_data) - copy_dir = Path(output_dir) / "copied" - copy_dir.mkdir(exist_ok=True, parents=True) - # Act - meteo_grid_copy.write(copy_dir) - # Assert - assert filecmp.cmp(grid_dir / "mete_grid.inp", copy_dir / "mete_grid.inp") +def test_meteogridcopy_from_imod5(meteo_grids, tmp_path): + # Arrange + grid_dir = setup_written_meteo_grids(meteo_grids, tmp_path) + imod5_data = {} + imod5_data["extra"] = {} + imod5_data["extra"]["paths"] = [ + ["foo"], + [(grid_dir / "mete_grid.inp").resolve()], + ["bar"], + ] + meteo_grid_copy = MeteoGridCopy.from_imod5_data(imod5_data) + copy_dir = tmp_path / "copied" + copy_dir.mkdir(exist_ok=True, parents=True) + # Act + meteo_grid_copy.write(copy_dir) + # Assert + assert filecmp.cmp(grid_dir / "mete_grid.inp", copy_dir / "mete_grid.inp") diff --git a/imod/tests/test_msw/test_model.py b/imod/tests/test_msw/test_model.py index 0e97a59c5..f0904e96f 100644 --- a/imod/tests/test_msw/test_model.py +++ b/imod/tests/test_msw/test_model.py @@ -1,11 +1,15 @@ +from copy import copy from pathlib import Path import pytest import xarray as xr from numpy.testing import assert_almost_equal, assert_equal +from pytest_cases import parametrize_with_cases from imod import mf6, msw from imod.mf6.utilities.regrid import RegridderWeightsCache +from imod.msw.model import DEFAULT_SETTINGS +from imod.typing import GridDataArray, Imod5DataDict def test_msw_model_write(msw_model, coupled_mf6_model, coupled_mf6wel, tmp_path): @@ -138,3 +142,156 @@ def test_coupled_model_regrid(msw_model, coupled_mf6_model, tmp_path): ) regridded_msw_model.write(tmp_path, mf6_discretization, regridded_mf6_wel) + + +def setup_written_meteo_grids( + meteo_grids: tuple[GridDataArray], tmp_path: Path +) -> Path: + precipitation, _ = meteo_grids + meteo_grid = msw.MeteoGrid(precipitation, precipitation) + grid_dir = tmp_path / "grid" + grid_dir.mkdir(exist_ok=True, parents=True) + meteo_grid.write(grid_dir) + return grid_dir + + +def setup_parasim_inp(directory: Path): + settings = copy(DEFAULT_SETTINGS) + settings["iybg"] = 2000 + settings["tdbg"] = 200.45 + settings["unsa_svat_path"] = str(directory) + + filename = directory / "para_sim.inp" + with open(filename, "w") as f: + rendered = msw.MetaSwapModel._template.render(settings=settings) + f.write(rendered) + return filename + + +def write_test_files(directory: Path, filenames: list[str]) -> list[Path]: + paths = [directory / filename for filename in filenames] + for p in paths: + with open(p, mode="w") as f: + f.write("test") + return paths + + +def setup_extra_files(meteo_grids: tuple[GridDataArray], directory: Path): + grid_dir = setup_written_meteo_grids(meteo_grids, directory) + setup_parasim_inp(grid_dir) + write_test_files(grid_dir, ["a.inp", "b.inp"]) + return { + "paths": [ + [str(grid_dir / fn)] + for fn in ["a.inp", "mete_grid.inp", "para_sim.inp", "b.inp"] + ] + } + + +class Imod5DataCases: + def case_grid(self, imod5_cap_data: Imod5DataDict) -> tuple[Imod5DataDict, bool]: + has_scaling_factor = True + return imod5_cap_data, has_scaling_factor + + def case_no_scaling_factors( + self, imod5_cap_data: Imod5DataDict + ) -> tuple[Imod5DataDict, bool]: + has_scaling_factor = False + cap_data = imod5_cap_data["cap"] + # open_projectfile_data adds layer kwargs to constants + layer_kwargs = {"coords": {"layer": [1]}, "dims": ("layer",)} + cap_data["perched_water_table_level"] = xr.DataArray([-9999.0], **layer_kwargs) + cap_data["soil_moisture_fraction"] = xr.DataArray([1.0], **layer_kwargs) + cap_data["conductivitiy_factor"] = xr.DataArray([1.0], **layer_kwargs) + return imod5_cap_data, has_scaling_factor + + def case_constants( + self, imod5_cap_data: Imod5DataDict + ) -> tuple[Imod5DataDict, bool]: + has_scaling_factor = False + cap_data = imod5_cap_data["cap"] + # open_projectfile_data adds layer kwargs to constants + layer_kwargs = {"coords": {"layer": [1]}, "dims": ("layer",)} + cap_data["perched_water_table_level"] = xr.DataArray([-9999.0], **layer_kwargs) + cap_data["soil_moisture_fraction"] = xr.DataArray([1.0], **layer_kwargs) + cap_data["conductivitiy_factor"] = xr.DataArray([1.0], **layer_kwargs) + cap_data["urban_ponding_depth"] = xr.DataArray([1.0], **layer_kwargs) + cap_data["rural_ponding_depth"] = xr.DataArray([1.0], **layer_kwargs) + cap_data["urban_runoff_resistance"] = xr.DataArray([1.0], **layer_kwargs) + cap_data["rural_runoff_resistance"] = xr.DataArray([1.0], **layer_kwargs) + cap_data["urban_runon_resistance"] = xr.DataArray([1.0], **layer_kwargs) + cap_data["rural_runon_resistance"] = xr.DataArray([1.0], **layer_kwargs) + cap_data["urban_infiltration_capacity"] = xr.DataArray([1.0], **layer_kwargs) + cap_data["rural_infiltration_capacity"] = xr.DataArray([1.0], **layer_kwargs) + return imod5_cap_data, has_scaling_factor + + +@parametrize_with_cases("imod5_data, has_scaling_factor", cases=Imod5DataCases) +def test_import_from_imod5( + imod5_data: Imod5DataDict, + has_scaling_factor: bool, + meteo_grids: tuple[GridDataArray], + coupled_mf6_model: mf6.Modflow6Simulation, + tmp_path: Path, +): + # Arrange + imod5_data["extra"] = setup_extra_files(meteo_grids, tmp_path) + times = coupled_mf6_model["time_discretization"].dataset.coords["time"] + dis_pkg = coupled_mf6_model["GWF_1"]["dis"] + # Act + model = msw.MetaSwapModel.from_imod5_data(imod5_data, dis_pkg, times) + # Assert + grid_packages = { + "grid", + "infiltration", + "ponding", + "sprinkling", + "idf_mapping", + } + expected_keys = grid_packages | { + "meteo_grid", + "prec_mapping", + "evt_mapping", + "coupling", + "extra_files", + "time_oc", + } + missing_keys = expected_keys - set(model.keys()) + assert not missing_keys + for pkgname in grid_packages: + # Test if all grid packages broadcasted to grid. + missing_dims = {"y", "x"} - set(model[pkgname].dataset.dims.keys()) + assert not missing_dims + assert "time" in model["time_oc"].dataset.dims.keys() + assert len(model["meteo_grid"].dataset.dims) == 0 + assert ("scaling_factor" in model.keys()) == has_scaling_factor + + +@parametrize_with_cases("imod5_data, has_scaling_factor", cases=Imod5DataCases) +def test_import_from_imod5_and_write( + imod5_data: Imod5DataDict, + has_scaling_factor: bool, + meteo_grids: tuple[GridDataArray], + coupled_mf6_model: mf6.Modflow6Simulation, + tmp_path: Path, +): + # Arrange + imod5_data["extra"] = setup_extra_files(meteo_grids, tmp_path) + times = coupled_mf6_model["time_discretization"].dataset.coords["time"] + dis_pkg = coupled_mf6_model["GWF_1"]["dis"] + npf_pkg = coupled_mf6_model["GWF_1"]["npf"] + active = dis_pkg["idomain"] == 1 + modeldir = tmp_path / "modeldir" + # Act + model = msw.MetaSwapModel.from_imod5_data(imod5_data, dis_pkg, times) + well_pkg = mf6.LayeredWell.from_imod5_cap_data(imod5_data) + mf6_wel_pkg = well_pkg.to_mf6_pkg( + active, dis_pkg["top"], dis_pkg["bottom"], npf_pkg["k"] + ) + model.write(modeldir, dis_pkg, mf6_wel_pkg, validate=False) + + # Assert + expected_n_files = 13 + if has_scaling_factor: + expected_n_files += 1 + assert len(list(modeldir.rglob(r"*.inp"))) == expected_n_files diff --git a/imod/tests/test_msw/test_utilities/test_mask.py b/imod/tests/test_msw/test_utilities/test_mask.py new file mode 100644 index 000000000..424f7dbc7 --- /dev/null +++ b/imod/tests/test_msw/test_utilities/test_mask.py @@ -0,0 +1,131 @@ +from copy import deepcopy +from typing import Callable + +import numpy as np +import pytest +import xarray as xr +from pytest_cases import parametrize_with_cases + +from imod.msw import GridData +from imod.msw.utilities.mask import ( + MetaSwapActive, + mask_and_broadcast_cap_data, + mask_and_broadcast_pkg_data, +) +from imod.typing import GridDataDict + + +@pytest.fixture(scope="function") +def coords_planar() -> dict: + x = [1.0, 2.0, 3.0] + y = [3.0, 2.0, 1.0] + dx = 1.0 + dy = 1.0 + return {"y": y, "x": x, "dx": dx, "dy": dy} + + +@pytest.fixture(scope="function") +def coords_subunit(coords_planar: dict) -> dict: + coords_subunit = deepcopy(coords_planar) + coords_subunit["subunit"] = [0, 1] + return coords_subunit + + +@pytest.fixture(scope="function") +def mask_fixture(coords_subunit: dict) -> MetaSwapActive: + mask_per_subunit = xr.DataArray( + np.array( + [ + [[0, 0, 0], [0, 1, 1], [0, 0, 0]], + [[1, 1, 1], [0, 1, 1], [0, 0, 0]], + ] + ).astype(bool), + dims=("subunit", "y", "x"), + coords=coords_subunit, + ) + mask_all = mask_per_subunit.any(dim="subunit") + + return MetaSwapActive(mask_all, mask_per_subunit) + + +@pytest.fixture(scope="function") +def data_subunit(coords_subunit: dict): + return xr.DataArray( + np.array( + [ + [[2.0, 2.0, 2.0], [2.0, 3.0, 3.0], [2.0, 2.0, 2.0]], + [[3.0, 3.0, 3.0], [2.0, 3.0, 3.0], [2.0, 2.0, 2.0]], + ] + ), + dims=("subunit", "y", "x"), + coords=coords_subunit, + ) + + +@pytest.fixture(scope="function") +def data_planar(data_subunit: xr.DataArray) -> xr.DataArray: + return data_subunit.sel(subunit=0, drop=True) + + +def integer_nan(da): + return da == 0 + + +def case_grid_float( + data_subunit: xr.DataArray, data_planar: xr.DataArray +) -> tuple[GridDataDict, Callable]: + return {"landuse": data_subunit, "elevation": data_planar}, np.isnan + + +def case_grid_integer( + data_subunit: xr.DataArray, data_planar: xr.DataArray +) -> tuple[GridDataDict, Callable]: + return { + "landuse": data_subunit.astype(int), + "elevation": data_planar.astype(int), + }, integer_nan + + +def case_constant_float(): + data = xr.DataArray(2.0) + return {"landuse": data, "elevation": data}, np.isnan + + +def case_constant_integer(): + data = xr.DataArray(2) + return {"landuse": data, "elevation": data}, integer_nan + + +@parametrize_with_cases("imod5_grid_data, isnan", cases=".") +def test_mask_and_broadcast_pkg_data( + imod5_grid_data: GridDataDict, isnan: Callable, mask_fixture: MetaSwapActive +): + # Act + masked_data = mask_and_broadcast_pkg_data(GridData, imod5_grid_data, mask_fixture) + + # Assert + assert isnan(masked_data["landuse"]).to_numpy().sum() == 11 + assert isnan(masked_data["elevation"]).to_numpy().sum() == 4 + + +@parametrize_with_cases("imod5_grid_data, isnan", cases=".") +def test_mask_and_broadcast_cap_data( + imod5_grid_data: GridDataDict, + isnan: Callable, + mask_fixture: MetaSwapActive, + coords_planar: dict, +): + # Arrange, cap data doesn't have subunit coords + imod5_grid_data["landuse"] = imod5_grid_data["landuse"].isel( + subunit=0, drop=True, missing_dims="ignore" + ) + + # Act + masked_data = mask_and_broadcast_cap_data(imod5_grid_data, mask_fixture) + + # Assert + coords_expected = xr.Coordinates(coords_planar) + xr.testing.assert_equal(masked_data["landuse"].coords, coords_expected) + xr.testing.assert_equal(masked_data["elevation"].coords, coords_expected) + assert isnan(masked_data["landuse"]).to_numpy().sum() == 4 + assert isnan(masked_data["elevation"]).to_numpy().sum() == 4 diff --git a/imod/tests/test_msw/test_utilities/test_parse.py b/imod/tests/test_msw/test_utilities/test_parse.py new file mode 100644 index 000000000..a4293252a --- /dev/null +++ b/imod/tests/test_msw/test_utilities/test_parse.py @@ -0,0 +1,28 @@ +from pytest_cases import parametrize_with_cases + +from imod.msw.utilities.parse import _try_parsing_string_to_number + + +class ParseCases: + def case_int(self): + return "1", int + + def case_float(self): + return "1.0", float + + def case_string(self): + return "a", str + + def case_aterisk(self): + return "*", str + + def case_exclamation(self): + return "!", str + + +@parametrize_with_cases(["s", "expected_type"], cases=ParseCases) +def test_try_parsing_string_to_value(s, expected_type): + # Act + parsed = _try_parsing_string_to_number(s) + # Assert + assert type(parsed) is expected_type diff --git a/imod/typing/__init__.py b/imod/typing/__init__.py index 0217f9782..da8bb6c19 100644 --- a/imod/typing/__init__.py +++ b/imod/typing/__init__.py @@ -12,7 +12,6 @@ GridDataArray: TypeAlias = Union[xr.DataArray, xu.UgridDataArray] GridDataset: TypeAlias = Union[xr.Dataset, xu.UgridDataset] GridDataDict: TypeAlias = dict[str, GridDataArray] -Imod5DataDict: TypeAlias = dict[str, GridDataDict | dict[str, list[str]]] ScalarAsDataArray: TypeAlias = Union[xr.DataArray, xu.UgridDataArray] ScalarAsDataset: TypeAlias = Union[xr.Dataset, xu.UgridDataset] UnstructuredData: TypeAlias = Union[xu.UgridDataset, xu.UgridDataArray] @@ -26,6 +25,11 @@ class SelSettingsType(TypedDict, total=False): missing_dims: Literal["raise", "warn", "ignore"] +class Imod5DataDict(TypedDict, total=False): + cap: GridDataDict + extra: dict[str, list[str]] + + # Types for optional dependencies. if TYPE_CHECKING: import geopandas as gpd From d3451cd483aad16aa63de60f27e6879dac66b949 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Thu, 12 Dec 2024 17:39:27 +0100 Subject: [PATCH 10/13] Finish docstrings for methods --- imod/msw/infiltration.py | 21 ++++++++++++++++++ imod/msw/meteo_grid.py | 16 ++++++++++++++ imod/msw/meteo_mapping.py | 44 ++++++++++++++++++++++++++++++------- imod/msw/ponding.py | 18 ++++++++++++++- imod/msw/scaling_factors.py | 19 ++++++++++++++-- 5 files changed, 107 insertions(+), 11 deletions(-) diff --git a/imod/msw/infiltration.py b/imod/msw/infiltration.py index 7e5472c7b..08a7ef18e 100644 --- a/imod/msw/infiltration.py +++ b/imod/msw/infiltration.py @@ -103,6 +103,27 @@ def __init__( @classmethod def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "Infiltration": + """ + Construct a MetaSWAP Infiltration package from iMOD5 data in the CAP + package, loaded with the :func:`imod.formats.prj.open_projectfile_data` + function. + + Concatenates infiltration_capacity, runon_resistance, and + runoff_resistance along the subunit dimension. 0 = rural landuse, 1 = + urban landuse. Resistances smaller than 5 days are deactivated for + transparency, as MetaSWAP also does this internally. + + Parameters + ---------- + imod5_data: Imod5DataDict + iMOD5 data as returned by + :func:`imod.formats.prj.open_projectfile_data` + + Returns + ------- + imod.msw.Infiltration + """ + cap_data = imod5_data["cap"] data = {} # Use runon resistance as downward resistance, and runoff for downward diff --git a/imod/msw/meteo_grid.py b/imod/msw/meteo_grid.py index 3974eb838..79c179c5c 100644 --- a/imod/msw/meteo_grid.py +++ b/imod/msw/meteo_grid.py @@ -226,6 +226,22 @@ def write(self, directory: Path | str, *args): @classmethod def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "MeteoGridCopy": + """ + Construct a MetaSWAP MeteoGridCopy package from iMOD5 data in the CAP + package, loaded with the :func:`imod.formats.prj.open_projectfile_data` + function. + + Parameters + ---------- + imod5_data: Imod5DataDict + iMOD5 data as returned by + :func:`imod.formats.prj.open_projectfile_data` + + Returns + ------- + imod.msw.MeteoGridCopy + """ + paths = imod5_data["extra"]["paths"] filepath = find_in_file_list(cls._file_name, paths) diff --git a/imod/msw/meteo_mapping.py b/imod/msw/meteo_mapping.py index b728bb497..b39fb8418 100644 --- a/imod/msw/meteo_mapping.py +++ b/imod/msw/meteo_mapping.py @@ -178,10 +178,24 @@ def __init__( @classmethod def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "PrecipitationMapping": """ - Construct precipitation mapping from imod5 data. Opens first ascii grid - in mete_grid.inp, which is used to construct mappings to svats. The - grids should not change in dimension over time. No checks are done - whether cells switch from inactive to active or vice versa. + Construct a MetaSWAP PrecipitationMapping package from iMOD5 data in the + CAP package, loaded with the + :func:`imod.formats.prj.open_projectfile_data` function. + + Opens first ascii grid in mete_grid.inp, which is used to construct + mappings to svats. The grids should not change in dimension over time. + No checks are done whether cells switch from inactive to active or vice + versa. + + Parameters + ---------- + imod5_data: Imod5DataDict + iMOD5 data as returned by + :func:`imod.formats.prj.open_projectfile_data` + + Returns + ------- + imod.msw.PrecipitationMapping """ column_nr = 2 meteo_grid = open_first_meteo_grid_from_imod5_data(imod5_data, column_nr) @@ -222,10 +236,24 @@ def __init__( @classmethod def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "EvapotranspirationMapping": """ - Construct evapotranspiration mapping from imod5 data. Opens first ascii - grid in mete_grid.inp, which is used to construct mappings to svats. The - grids should not change in dimension over time. No checks are done - whether cells switch from inactive to active or vice versa. + Construct a MetaSWAP EvapotranspirationMapping package from iMOD5 data + in the CAP package, loaded with the + :func:`imod.formats.prj.open_projectfile_data` function. + + Opens first ascii grid in mete_grid.inp, which is used to construct + mappings to svats. The grids should not change in dimension over time. + No checks are done whether cells switch from inactive to active or vice + versa. + + Parameters + ---------- + imod5_data: Imod5DataDict + iMOD5 data as returned by + :func:`imod.formats.prj.open_projectfile_data` + + Returns + ------- + imod.msw.EvapotranspirationMapping """ column_nr = 3 meteo_grid = open_first_meteo_grid_from_imod5_data(imod5_data, column_nr) diff --git a/imod/msw/ponding.py b/imod/msw/ponding.py index 86a7c5cff..9248eb1dd 100644 --- a/imod/msw/ponding.py +++ b/imod/msw/ponding.py @@ -75,7 +75,23 @@ def _render(self, file: TextIO, index: IntArray, svat: xr.DataArray, *args: Any) @classmethod def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "Ponding": """ - Concatenate ponding depths along subunits + Construct a MetaSWAP Ponding package from iMOD5 data in the CAP + package, loaded with the :func:`imod.formats.prj.open_projectfile_data` + function. + + Method concatenates ponding depths, runon resistance, and runoff + resistance along two subunits. Subunit 0 for rural, and 1 for urban + landuse. + + Parameters + ---------- + imod5_data: Imod5DataDict + iMOD5 data as returned by + :func:`imod.formats.prj.open_projectfile_data` + + Returns + ------- + imod.msw.Ponding """ cap_data = imod5_data["cap"] data = {} diff --git a/imod/msw/scaling_factors.py b/imod/msw/scaling_factors.py index 91984dd3d..7e0053eb6 100644 --- a/imod/msw/scaling_factors.py +++ b/imod/msw/scaling_factors.py @@ -83,8 +83,23 @@ def __init__( @classmethod def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "ScalingFactors": """ - Import ScalingFactors from iMOD5 data. Pressure head factor is set to - one for all factors, as well as all factors for urban areas. + Construct a MetaSWAP Ponding package from iMOD5 data in the CAP + package, loaded with the :func:`imod.formats.prj.open_projectfile_data` + function. + + Pressure head factor is set to one for all subunits. All urban areas + (subunit=1) are also set to one for all variables, except the perced + water table level. + + Parameters + ---------- + imod5_data: Imod5DataDict + iMOD5 data as returned by + :func:`imod.formats.prj.open_projectfile_data` + + Returns + ------- + imod.msw.ScalingFactors """ cap_data = imod5_data["cap"] grid_ones = ones_like(cap_data["boundary"]) From 41bbfe2e1f629e6db588ba89fcf525457ffe04da Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Fri, 13 Dec 2024 17:32:05 +0100 Subject: [PATCH 11/13] format --- imod/msw/infiltration.py | 2 +- imod/msw/meteo_grid.py | 2 +- imod/msw/meteo_mapping.py | 8 ++++---- imod/msw/ponding.py | 2 +- imod/msw/scaling_factors.py | 4 ++-- 5 files changed, 9 insertions(+), 9 deletions(-) diff --git a/imod/msw/infiltration.py b/imod/msw/infiltration.py index 08a7ef18e..a826eb502 100644 --- a/imod/msw/infiltration.py +++ b/imod/msw/infiltration.py @@ -118,7 +118,7 @@ def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "Infiltration": imod5_data: Imod5DataDict iMOD5 data as returned by :func:`imod.formats.prj.open_projectfile_data` - + Returns ------- imod.msw.Infiltration diff --git a/imod/msw/meteo_grid.py b/imod/msw/meteo_grid.py index 79c179c5c..805619425 100644 --- a/imod/msw/meteo_grid.py +++ b/imod/msw/meteo_grid.py @@ -236,7 +236,7 @@ def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "MeteoGridCopy": imod5_data: Imod5DataDict iMOD5 data as returned by :func:`imod.formats.prj.open_projectfile_data` - + Returns ------- imod.msw.MeteoGridCopy diff --git a/imod/msw/meteo_mapping.py b/imod/msw/meteo_mapping.py index b39fb8418..a8bbd1bc2 100644 --- a/imod/msw/meteo_mapping.py +++ b/imod/msw/meteo_mapping.py @@ -181,7 +181,7 @@ def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "PrecipitationMapping": Construct a MetaSWAP PrecipitationMapping package from iMOD5 data in the CAP package, loaded with the :func:`imod.formats.prj.open_projectfile_data` function. - + Opens first ascii grid in mete_grid.inp, which is used to construct mappings to svats. The grids should not change in dimension over time. No checks are done whether cells switch from inactive to active or vice @@ -192,7 +192,7 @@ def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "PrecipitationMapping": imod5_data: Imod5DataDict iMOD5 data as returned by :func:`imod.formats.prj.open_projectfile_data` - + Returns ------- imod.msw.PrecipitationMapping @@ -239,7 +239,7 @@ def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "EvapotranspirationMappin Construct a MetaSWAP EvapotranspirationMapping package from iMOD5 data in the CAP package, loaded with the :func:`imod.formats.prj.open_projectfile_data` function. - + Opens first ascii grid in mete_grid.inp, which is used to construct mappings to svats. The grids should not change in dimension over time. No checks are done whether cells switch from inactive to active or vice @@ -250,7 +250,7 @@ def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "EvapotranspirationMappin imod5_data: Imod5DataDict iMOD5 data as returned by :func:`imod.formats.prj.open_projectfile_data` - + Returns ------- imod.msw.EvapotranspirationMapping diff --git a/imod/msw/ponding.py b/imod/msw/ponding.py index 9248eb1dd..a7e37a429 100644 --- a/imod/msw/ponding.py +++ b/imod/msw/ponding.py @@ -88,7 +88,7 @@ def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "Ponding": imod5_data: Imod5DataDict iMOD5 data as returned by :func:`imod.formats.prj.open_projectfile_data` - + Returns ------- imod.msw.Ponding diff --git a/imod/msw/scaling_factors.py b/imod/msw/scaling_factors.py index 7e0053eb6..e5d088e50 100644 --- a/imod/msw/scaling_factors.py +++ b/imod/msw/scaling_factors.py @@ -90,13 +90,13 @@ def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "ScalingFactors": Pressure head factor is set to one for all subunits. All urban areas (subunit=1) are also set to one for all variables, except the perced water table level. - + Parameters ---------- imod5_data: Imod5DataDict iMOD5 data as returned by :func:`imod.formats.prj.open_projectfile_data` - + Returns ------- imod.msw.ScalingFactors From 006ff0480735e6322b49f1b56b6f11bf1698ce0d Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Fri, 13 Dec 2024 17:43:03 +0100 Subject: [PATCH 12/13] Add missing assertion statement --- imod/tests/test_msw/test_copy_files.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/imod/tests/test_msw/test_copy_files.py b/imod/tests/test_msw/test_copy_files.py index 80d8ceb08..89e136286 100644 --- a/imod/tests/test_msw/test_copy_files.py +++ b/imod/tests/test_msw/test_copy_files.py @@ -66,4 +66,4 @@ def test_from_imod5_data(src_files): # Act copyfiles = FileCopier.from_imod5_data(imod5_data) # Assert - len(copyfiles.dataset["paths"]) == 3 + assert len(copyfiles.dataset["paths"]) == 3 From 4c4a6f4a5448f00e645a0ef45dad46ba538b9baf Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Fri, 13 Dec 2024 17:48:31 +0100 Subject: [PATCH 13/13] Remove unreachable return statements --- imod/mf6/utilities/imod5_converter.py | 1 - imod/msw/sprinkling.py | 1 - 2 files changed, 2 deletions(-) diff --git a/imod/mf6/utilities/imod5_converter.py b/imod/mf6/utilities/imod5_converter.py index c3301430f..d7fb73bc4 100644 --- a/imod/mf6/utilities/imod5_converter.py +++ b/imod/mf6/utilities/imod5_converter.py @@ -56,7 +56,6 @@ def _well_from_imod5_cap_point_data(cap_data: GridDataDict) -> dict[str, np.ndar 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]: diff --git a/imod/msw/sprinkling.py b/imod/msw/sprinkling.py index b5fa8cc69..2687b5802 100644 --- a/imod/msw/sprinkling.py +++ b/imod/msw/sprinkling.py @@ -26,7 +26,6 @@ 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: