diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst index 7f7d4b2cd..76ce99985 100644 --- a/docs/api/changelog.rst +++ b/docs/api/changelog.rst @@ -12,6 +12,21 @@ The format is based on `Keep a Changelog`_, and this project adheres to 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. +- :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. - :func:`imod.data.tutorial_03` to load data for the iMOD Documentation tutorial. @@ -35,8 +50,13 @@ Fixed Changed ~~~~~~~ +- :class:`imod.msw.Infiltration`'s variables ``upward_resistance`` and + ``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/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 3197cfa91..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, @@ -303,9 +305,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] @@ -330,8 +332,11 @@ def from_imod5_data( wel_key, imod5_data, times ) + 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/rch.py b/imod/mf6/rch.py index ba917cbd8..3fd892c90 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,34 @@ 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) + msw_active = is_msw_active_cell(target_dis, cap_data, msw_area) + active = msw_active.all + + 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/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/utilities/imod5_converter.py b/imod/mf6/utilities/imod5_converter.py index e599c8522..d7fb73bc4 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,71 @@ 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." + ) + + +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..11334b2cf 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,48 @@ 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/__init__.py b/imod/msw/__init__.py index bcbe174aa..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 @@ -9,7 +10,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/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/msw/grid_data.py b/imod/msw/grid_data.py index 25003716b..d4c387bd3 100644 --- a/imod/msw/grid_data.py +++ b/imod/msw/grid_data.py @@ -1,10 +1,19 @@ import numpy as np import xarray as xr +from imod.mf6 import StructuredDiscretization 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 GridDataRegridMethod +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.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 @@ -111,3 +120,52 @@ def _pkgcheck(self): raise ValueError( "Provided area grid with total areas larger than cell area" ) + + @classmethod + def from_imod5_data( + cls, + imod5_data: Imod5DataDict, + target_dis: StructuredDiscretization, + ) -> 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 = {} + 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"].astype(int) + + 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 bf300761f..a826eb502 100644 --- a/imod/msw/infiltration.py +++ b/imod/msw/infiltration.py @@ -1,9 +1,36 @@ +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.msw.utilities.mask import MaskValues +from imod.typing import GridDataDict, Imod5DataDict +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, MaskValues.default) + return data class Infiltration(MetaSwapPackage, IRegridPackage): @@ -19,11 +46,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 +59,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 +71,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 +100,49 @@ def __init__( self.dataset["extra_storage_coefficient"] = extra_storage_coefficient self._pkgcheck() + + @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 + # 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) * 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 d8c283252..805619425 100644 --- a/imod/msw/meteo_grid.py +++ b/imod/msw/meteo_grid.py @@ -1,5 +1,6 @@ import csv from pathlib import Path +from shutil import copyfile from typing import Optional, Union import numpy as np @@ -11,6 +12,10 @@ 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.msw.utilities.mask import MaskValues +from imod.typing import Imod5DataDict +from imod.util.regrid_method_type import EmptyRegridMethod, RegridMethodType class MeteoGrid(MetaSwapPackage, IRegridPackage): @@ -172,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: @@ -188,3 +195,54 @@ 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": + """ + 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) + + return cls(filepath) diff --git a/imod/msw/meteo_mapping.py b/imod/msw/meteo_mapping.py index 0fa3c818c..a8bbd1bc2 100644 --- a/imod/msw/meteo_mapping.py +++ b/imod/msw/meteo_mapping.py @@ -1,18 +1,68 @@ from copy import deepcopy +from pathlib import Path +from textwrap import dedent from typing import Any, Optional, TextIO 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 = 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,32 @@ def __init__( super().__init__() self.meteo = precipitation + @classmethod + 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 + 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) + return cls(meteo_grid) + class EvapotranspirationMapping(MeteoMapping): """ @@ -156,3 +232,29 @@ def __init__( ): super().__init__() self.meteo = evapotranspiration + + @classmethod + def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "EvapotranspirationMapping": + """ + 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) + return cls(meteo_grid) 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/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/msw/ponding.py b/imod/msw/ponding.py index afc1629af..a7e37a429 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 Imod5DataDict, IntArray class Ponding(MetaSwapPackage, IRegridPackage): @@ -70,3 +71,32 @@ 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: Imod5DataDict) -> "Ponding": + """ + 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 = {} + 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/scaling_factors.py b/imod/msw/scaling_factors.py index 6e544d1c5..e5d088e50 100644 --- a/imod/msw/scaling_factors.py +++ b/imod/msw/scaling_factors.py @@ -2,6 +2,9 @@ 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 Imod5DataDict +from imod.typing.grid import ones_like class ScalingFactors(MetaSwapPackage, IRegridPackage): @@ -76,3 +79,39 @@ def __init__( self.dataset["depth_perched_water_table"] = depth_perched_water_table self._pkgcheck() + + @classmethod + def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "ScalingFactors": + """ + 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"]) + + 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/msw/sprinkling.py b/imod/msw/sprinkling.py index 5d8bd20f1..2687b5802 100644 --- a/imod/msw/sprinkling.py +++ b/imod/msw/sprinkling.py @@ -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,48 @@ 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." + ) + + +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 +166,53 @@ 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 = 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/msw/utilities/common.py b/imod/msw/utilities/common.py new file mode 100644 index 000000000..db69aab33 --- /dev/null +++ b/imod/msw/utilities/common.py @@ -0,0 +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/msw/utilities/imod5_converter.py b/imod/msw/utilities/imod5_converter.py new file mode 100644 index 000000000..0c4d21604 --- /dev/null +++ b/imod/msw/utilities/imod5_converter.py @@ -0,0 +1,119 @@ +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 + + +def get_cell_area_from_imod5_data( + imod5_cap: GridDataDict, +) -> GridDataArray: + # area's per type of svats + 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) + 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).astype(int) + + +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, +) -> MetaSwapActive: + """ + 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 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 593ffb5bb..b0f400192 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, @@ -88,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/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/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/fixtures/msw_model_fixture.py b/imod/tests/fixtures/msw_model_fixture.py index ceb33e35a..33b9e8723 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_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() 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_copy_files.py b/imod/tests/test_msw/test_copy_files.py new file mode 100644 index 000000000..89e136286 --- /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 + assert len(copyfiles.dataset["paths"]) == 3 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_grid_data.py b/imod/tests/test_msw/test_grid_data.py index a3b62a64b..a775a93a9 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/tests/test_msw/test_infiltration.py b/imod/tests/test_msw/test_infiltration.py index d3cb0a003..47990e3a6 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) + + 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, 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 8cca2cc06..66d7fe95e 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 @@ -6,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 - - -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 +from imod.msw import MeteoGrid, MeteoGridCopy + - 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) @@ -82,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 @@ -113,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() @@ -123,3 +70,42 @@ 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 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 + + +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_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_ponding.py b/imod/tests/test_msw/test_ponding.py index 626cd6412..94f0d64d0 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,27 @@ 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), + ("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) 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/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]) + ) 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/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 0d17b1127..da8bb6c19 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 @@ -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] @@ -18,6 +19,17 @@ IntArray: TypeAlias = NDArray[np.int_] +class SelSettingsType(TypedDict, total=False): + layer: int + drop: bool + 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