diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst index 1f073881c..f3927430e 100644 --- a/docs/api/changelog.rst +++ b/docs/api/changelog.rst @@ -21,6 +21,12 @@ Added MetaSWAP). Currently only griddata (IDF) is supported. - :meth:`imod.mf6.Recharge.from_imod5_cap_data` to construct a recharge package for coupling a MODFLOW6 model to MetaSWAP. +- :meth:`imod.msw.MetaSwapModel.from_imod5_data` to construct a MetaSWAP model + from data in an iMOD5 projectfile. +- :meth:`imod.msw.MetaSwapModel.write` has a ``validate`` argument, which can be + used to turn off validation upon writing, use at your own risk! +- :class:`imod.msw.MetaSwapModel` got ``settings`` argument to set simulation + settings. Fixed ~~~~~ diff --git a/docs/api/msw.rst b/docs/api/msw.rst index f2260ed0f..067710445 100644 --- a/docs/api/msw.rst +++ b/docs/api/msw.rst @@ -7,10 +7,15 @@ MetaSWAP :toctree: generated/msw GridData + GridData.from_imod5_data Infiltration + Infiltration.from_imod5_data Ponding + Ponding.from_imod5_data ScalingFactors + ScalingFactors.from_imod5_data Sprinkling + Sprinkling.from_imod5_data IdfMapping TimeOutputControl @@ -25,9 +30,15 @@ MetaSWAP AnnualCropFactors MeteoGrid + MeteoGridCopy + MeteoGridCopy.from_imod5_data EvapotranspirationMapping + EvapotranspirationMapping.from_imod5_data PrecipitationMapping + PrecipitationMapping.from_imod5_data CouplerMapping MetaSwapModel + MetaSwapModel.write + MetaSwapModel.from_imod5_data diff --git a/imod/mf6/model_gwf.py b/imod/mf6/model_gwf.py index e83ee5e1c..ebc878876 100644 --- a/imod/mf6/model_gwf.py +++ b/imod/mf6/model_gwf.py @@ -218,9 +218,11 @@ def from_imod5_data( times: list[datetime] Time discretization of the simulation. These times are used for the following: - * Times of wells with associated timeseries are resampled to these times - * Start- and end times in the list are used to repeat the stresses - of periodic data (e.g. river stages in iMOD5 for "summer", "winter") + + * Times of wells with associated timeseries are resampled to these times + * Start- and end times in the list are used to repeat the stresses + of periodic data (e.g. river stages in iMOD5 for "summer", "winter") + allocation_options: SimulationAllocationOptions object containing the allocation options per package type. If you want a package to have a different allocation option, diff --git a/imod/mf6/rch.py b/imod/mf6/rch.py index 392cc9d56..3fd892c90 100644 --- a/imod/mf6/rch.py +++ b/imod/mf6/rch.py @@ -249,7 +249,8 @@ def from_imod5_cap_data( cap_data = cast(GridDataDict, imod5_data["cap"]) msw_area = get_cell_area_from_imod5_data(cap_data) - active, _ = is_msw_active_cell(target_dis, cap_data, msw_area) + msw_active = is_msw_active_cell(target_dis, cap_data, msw_area) + active = msw_active.all data = {} layer_da = xr.full_like( diff --git a/imod/mf6/simulation.py b/imod/mf6/simulation.py index 5bc2469f1..4f324f31b 100644 --- a/imod/mf6/simulation.py +++ b/imod/mf6/simulation.py @@ -1361,11 +1361,13 @@ def from_imod5_data( times: list[datetime] time discretization of the model to be imported. This is used for the following: + * Times of wells with associated timeseries are resampled to these times * Start and end times in the list are used to repeat the stresses of periodic data (e.g. river stages in iMOD5 for "summer", "winter") * The simulation is discretized to these times, using :meth:`imod.mf6.Modflow6Simulation.create_time_discretization` + allocation_options: SimulationAllocationOptions, optional object containing the allocation options per package type. If you want a specific package to have a different allocation option, then diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py index d30b58d30..11334b2cf 100644 --- a/imod/mf6/wel.py +++ b/imod/mf6/wel.py @@ -1312,9 +1312,11 @@ def from_imod5_cap_data(cls, imod5_data: Imod5DataDict): (IDF), this grid defines in which layer a groundwater abstraction well should be placed. The ``"artificial_recharge"`` grid contains types which point to the type of abstraction: + * 0: no abstraction * 1: groundwater abstraction * 2: surfacewater abstraction + The ``"artificial_recharge_capacity"`` grid/constant defines the capacity of each groundwater or surfacewater abstraction. This is an ``1:1`` mapping: Each grid cell maps to a separate well. diff --git a/imod/msw/grid_data.py b/imod/msw/grid_data.py index ad7299206..d4c387bd3 100644 --- a/imod/msw/grid_data.py +++ b/imod/msw/grid_data.py @@ -1,14 +1,8 @@ -from typing import Optional - import numpy as np import xarray as xr from imod.mf6 import StructuredDiscretization from imod.mf6.interfaces.iregridpackage import IRegridPackage -from imod.mf6.regrid.regrid_schemes import ( - RegridMethodType, -) -from imod.mf6.utilities.regrid import RegridderWeightsCache from imod.msw.fixed_format import VariableMetaData from imod.msw.pkgbase import MetaSwapPackage from imod.msw.regrid.regrid_schemes import GridDataRegridMethod @@ -18,7 +12,8 @@ get_rootzone_depth_from_imod5_data, is_msw_active_cell, ) -from imod.typing import GridDataDict +from imod.msw.utilities.mask import MetaSwapActive, mask_and_broadcast_pkg_data +from imod.typing import Imod5DataDict from imod.util.spatial import get_cell_area, spatial_reference @@ -129,12 +124,38 @@ def _pkgcheck(self): @classmethod def from_imod5_data( cls, - imod5_data: dict[str, GridDataDict], + imod5_data: Imod5DataDict, target_dis: StructuredDiscretization, - regridder_types: Optional[RegridMethodType] = None, - regrid_cache: RegridderWeightsCache = RegridderWeightsCache(), - ) -> "GridData": - # Get iMOD5 capillary zone data + ) -> tuple["GridData", MetaSwapActive]: + """ + Construct a MetaSWAP GridData package from iMOD5 data in the CAP + package, loaded with the :func:`imod.formats.prj.open_projectfile_data` + function. + + Method does the following things: + + - Computes rural area from the wetted area and urban area grids. + - Sets a second subunit for urban landuse (== 18). + - Rootzone depth is converted from cm's to m's. + - Determines where MetaSWAP should be active based on area grids, + boundary grid, and MODFLOW6 idomain. + + Parameters + ---------- + imod5_data: Imod5DataDict + iMOD5 data as returned by + :func:`imod.formats.prj.open_projectfile_data` + target_dis: imod.mf6.StructuredDiscretization + MODFLOW6 discretization to couple the MetaSWAP model to. + + Returns + ------- + imod.msw.GridData + GridData package + MetaSwapActive + Dataclass containing where MetaSwAP is active, per subunit, as well + as aggregated over subunits. + """ imod5_cap = imod5_data["cap"] data = {} @@ -142,17 +163,9 @@ def from_imod5_data( data["landuse"] = get_landuse_from_imod5_data(imod5_cap) data["rootzone_depth"] = get_rootzone_depth_from_imod5_data(imod5_cap) data["surface_elevation"] = imod5_cap["surface_elevation"] - data["soil_physical_unit"] = imod5_cap["soil_physical_unit"] - - active, subunit_active = is_msw_active_cell(target_dis, imod5_cap, data["area"]) + data["soil_physical_unit"] = imod5_cap["soil_physical_unit"].astype(int) - data_active = { - key: ( - griddata.where(subunit_active) - if key in cls._with_subunit - else griddata.where(active) - ) - for key, griddata in data.items() - } - data_active["active"] = active - return cls(**data_active) + msw_active = is_msw_active_cell(target_dis, imod5_cap, data["area"]) + data_active = mask_and_broadcast_pkg_data(cls, data, msw_active) + data_active["active"] = msw_active.all + return cls(**data_active), msw_active diff --git a/imod/msw/infiltration.py b/imod/msw/infiltration.py index 830365681..7e5472c7b 100644 --- a/imod/msw/infiltration.py +++ b/imod/msw/infiltration.py @@ -8,7 +8,8 @@ from imod.msw.pkgbase import MetaSwapPackage from imod.msw.regrid.regrid_schemes import InfiltrationRegridMethod from imod.msw.utilities.common import concat_imod5 -from imod.typing import GridDataDict +from imod.msw.utilities.mask import MaskValues +from imod.typing import GridDataDict, Imod5DataDict from imod.typing.grid import ones_like @@ -28,7 +29,7 @@ def deactivate_small_resistances_in_data(data: GridDataDict): message=message.format(var=var), additional_depth=1, ) - data[var] = data[var].where(~to_deactivate, -9999.0) + data[var] = data[var].where(~to_deactivate, MaskValues.default) return data @@ -101,7 +102,7 @@ def __init__( self._pkgcheck() @classmethod - def from_imod5_data(cls, imod5_data: dict[str, GridDataDict]) -> "Infiltration": + def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "Infiltration": cap_data = imod5_data["cap"] data = {} # Use runon resistance as downward resistance, and runoff for downward @@ -120,7 +121,7 @@ def from_imod5_data(cls, imod5_data: dict[str, GridDataDict]) -> "Infiltration": data = deactivate_small_resistances_in_data(data) like = data["downward_resistance"].isel(subunit=0, drop=True) - data["bottom_resistance"] = ones_like(like) + data["bottom_resistance"] = ones_like(like) * MaskValues.default data["extra_storage_coefficient"] = ones_like(like) return cls(**data) diff --git a/imod/msw/meteo_grid.py b/imod/msw/meteo_grid.py index 70d3625a5..3974eb838 100644 --- a/imod/msw/meteo_grid.py +++ b/imod/msw/meteo_grid.py @@ -1,7 +1,7 @@ import csv from pathlib import Path from shutil import copyfile -from typing import Optional, Union, cast +from typing import Optional, Union import numpy as np import pandas as pd @@ -13,6 +13,7 @@ from imod.msw.regrid.regrid_schemes import MeteoGridRegridMethod from imod.msw.timeutil import to_metaswap_timeformat from imod.msw.utilities.common import find_in_file_list +from imod.msw.utilities.mask import MaskValues from imod.typing import Imod5DataDict from imod.util.regrid_method_type import EmptyRegridMethod, RegridMethodType @@ -176,7 +177,9 @@ def write(self, directory: Union[str, Path], *args): path = (directory / self._meteo_dirname / str(varname)).with_suffix( ".asc" ) - imod.rasterio.save(path, self.dataset[str(varname)], nodata=-9999.0) + imod.rasterio.save( + path, self.dataset[str(varname)], nodata=MaskValues.default + ) def _pkgcheck(self): for varname in self.dataset.data_vars: @@ -223,7 +226,7 @@ def write(self, directory: Path | str, *args): @classmethod def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "MeteoGridCopy": - paths = cast(list[str], imod5_data["extra"]["paths"]) + paths = imod5_data["extra"]["paths"] filepath = find_in_file_list(cls._file_name, paths) return cls(filepath) diff --git a/imod/msw/meteo_mapping.py b/imod/msw/meteo_mapping.py index 1de00fa69..b728bb497 100644 --- a/imod/msw/meteo_mapping.py +++ b/imod/msw/meteo_mapping.py @@ -1,7 +1,7 @@ from copy import deepcopy from pathlib import Path from textwrap import dedent -from typing import Any, Optional, TextIO, cast +from typing import Any, Optional, TextIO import numpy as np import pandas as pd @@ -58,7 +58,7 @@ def open_first_meteo_grid(mete_grid_path: str | Path, column_nr: int) -> xr.Data def open_first_meteo_grid_from_imod5_data(imod5_data: Imod5DataDict, column_nr: int): - paths = cast(list[str], imod5_data["extra"]["paths"]) + paths = imod5_data["extra"]["paths"] metegrid_path = find_in_file_list("mete_grid.inp", paths) return open_first_meteo_grid(metegrid_path, column_nr=column_nr) diff --git a/imod/msw/model.py b/imod/msw/model.py index f349021bb..6075ca59a 100644 --- a/imod/msw/model.py +++ b/imod/msw/model.py @@ -1,14 +1,17 @@ import collections from copy import copy +from datetime import datetime from pathlib import Path -from typing import Optional, Tuple, Union +from typing import Any, Optional, Tuple, Union, cast import jinja2 import numpy as np +import xarray as xr from imod.mf6.dis import StructuredDiscretization from imod.mf6.mf6_wel_adapter import Mf6Wel from imod.mf6.utilities.regrid import RegridderWeightsCache +from imod.msw.copy_files import FileCopier from imod.msw.coupler_mapping import CouplerMapping from imod.msw.grid_data import GridData from imod.msw.idf_mapping import IdfMapping @@ -20,11 +23,20 @@ InitialConditionsSavedState, ) from imod.msw.landuse import LanduseOptions -from imod.msw.meteo_grid import MeteoGrid +from imod.msw.meteo_grid import MeteoGrid, MeteoGridCopy from imod.msw.meteo_mapping import EvapotranspirationMapping, PrecipitationMapping from imod.msw.output_control import TimeOutputControl +from imod.msw.ponding import Ponding +from imod.msw.scaling_factors import ScalingFactors +from imod.msw.sprinkling import Sprinkling from imod.msw.timeutil import to_metaswap_timeformat +from imod.msw.utilities.common import find_in_file_list +from imod.msw.utilities.imod5_converter import has_active_scaling_factor +from imod.msw.utilities.mask import MaskValues, mask_and_broadcast_cap_data +from imod.msw.utilities.parse import read_para_sim +from imod.msw.utilities.select import drop_layer_dim_cap_data from imod.msw.vegetation import AnnualCropFactors +from imod.typing import Imod5DataDict from imod.util.regrid_method_type import RegridderType REQUIRED_PACKAGES = ( @@ -88,6 +100,7 @@ class MetaSwapModel(Model): ---------- unsaturated_database: Path-like or str Path to the MetaSWAP soil physical database folder. + settings: dict """ _pkg_id = "model" @@ -99,10 +112,18 @@ class MetaSwapModel(Model): "{%endfor%}" ) - def __init__(self, unsaturated_database): + def __init__( + self, + unsaturated_database: Path | str, + settings: Optional[dict[str, Any]] = None, + ): super().__init__() - self.simulation_settings = copy(DEFAULT_SETTINGS) + if settings is None: + self.simulation_settings = copy(DEFAULT_SETTINGS) + else: + self.simulation_settings = settings + self.simulation_settings["unsa_svat_path"] = ( self._render_unsaturated_database_path(unsaturated_database) ) @@ -210,6 +231,7 @@ def write( directory: Union[str, Path], mf6_dis: StructuredDiscretization, mf6_wel: Mf6Wel, + validate: bool = True, ): """ Write packages and simulation settings (para_sim.inp). @@ -221,9 +243,10 @@ def write( """ # Model checks - self._check_required_packages() - self._check_vegetation_indices_in_annual_crop_factors() - self._check_landuse_indices_in_lookup_options() + if validate: + self._check_required_packages() + self._check_vegetation_indices_in_annual_crop_factors() + self._check_landuse_indices_in_lookup_options() # Force to Path directory = Path(directory) @@ -258,7 +281,7 @@ def regrid_like( regrid_context: Optional[RegridderWeightsCache] = None, regridder_types: Optional[dict[str, Tuple[RegridderType, str]]] = None, ) -> "MetaSwapModel": - unsat_database = self.simulation_settings["unsa_svat_path"] + unsat_database = cast(str, self.simulation_settings["unsa_svat_path"]) regridded_model = MetaSwapModel(unsat_database) target_grid = mf6_regridded_dis["idomain"] @@ -282,3 +305,66 @@ def regrid_like( regridded_model[mod2svat_name] = CouplerMapping() return regridded_model + + @classmethod + def from_imod5_data( + cls, + imod5_data: Imod5DataDict, + target_dis: StructuredDiscretization, + times: list[datetime], + ): + """ + Construct a MetaSWAP model from iMOD5 data in the CAP package, loaded + with the :func:`imod.formats.prj.open_projectfile_data` function. + + Parameters + ---------- + imod5_data: dict + Dictionary with iMOD5 data. This can be constructed from the + :func:`imod.formats.prj.open_projectfile_data` method. + target_dis: imod.mf6.StructuredDiscretization + Target discretization, cells where MODLOW6 is inactive will be + inactive in MetaSWAP as well. + times: list[datetime] + List of datetimes, will be used to set the output control times. + Is also used to infer the starttime of the simulation. + + Returns + ------- + MetaSwapModel + MetaSWAP model imported from imod5 data. + """ + extra_paths = imod5_data["extra"]["paths"] + path_to_parasim = find_in_file_list("para_sim.inp", extra_paths) + parasim_settings = read_para_sim(path_to_parasim) + unsa_svat_path = cast(str, parasim_settings["unsa_svat_path"]) + # Drop layer coord + imod5_cap_no_layer = drop_layer_dim_cap_data(imod5_data) + model = cls(unsa_svat_path, parasim_settings) + model["grid"], msw_active = GridData.from_imod5_data( + imod5_cap_no_layer, target_dis + ) + cap_data_masked = mask_and_broadcast_cap_data( + imod5_cap_no_layer["cap"], msw_active + ) + imod5_masked: Imod5DataDict = { + "cap": cap_data_masked, + "extra": {"paths": extra_paths}, + } + model["infiltration"] = Infiltration.from_imod5_data(imod5_masked) + model["ponding"] = Ponding.from_imod5_data(imod5_masked) + model["sprinkling"] = Sprinkling.from_imod5_data(imod5_masked) + model["meteo_grid"] = MeteoGridCopy.from_imod5_data(imod5_masked) + model["prec_mapping"] = PrecipitationMapping.from_imod5_data(imod5_masked) + model["evt_mapping"] = EvapotranspirationMapping.from_imod5_data(imod5_masked) + if has_active_scaling_factor(imod5_cap_no_layer["cap"]): + model["scaling_factor"] = ScalingFactors.from_imod5_data(imod5_masked) + area = model["grid"]["area"].isel(subunit=0, drop=True) + model["idf_mapping"] = IdfMapping(area, MaskValues.default) + model["coupling"] = CouplerMapping() + model["extra_files"] = FileCopier.from_imod5_data(imod5_masked) + + times_da = xr.DataArray(times, coords={"time": times}, dims=("time",)) + model["time_oc"] = TimeOutputControl(times_da) + + return model diff --git a/imod/msw/ponding.py b/imod/msw/ponding.py index 89ab6af44..86a7c5cff 100644 --- a/imod/msw/ponding.py +++ b/imod/msw/ponding.py @@ -8,7 +8,7 @@ from imod.msw.pkgbase import DataDictType, MetaSwapPackage from imod.msw.regrid.regrid_schemes import PondingRegridMethod from imod.msw.utilities.common import concat_imod5 -from imod.typing import GridDataDict, IntArray +from imod.typing import Imod5DataDict, IntArray class Ponding(MetaSwapPackage, IRegridPackage): @@ -73,7 +73,7 @@ def _render(self, file: TextIO, index: IntArray, svat: xr.DataArray, *args: Any) return self.write_dataframe_fixed_width(file, dataframe) @classmethod - def from_imod5_data(cls, imod5_data: dict[str, GridDataDict]) -> "Ponding": + def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "Ponding": """ Concatenate ponding depths along subunits """ diff --git a/imod/msw/scaling_factors.py b/imod/msw/scaling_factors.py index 87dccaf61..91984dd3d 100644 --- a/imod/msw/scaling_factors.py +++ b/imod/msw/scaling_factors.py @@ -1,11 +1,9 @@ -from typing import cast - from imod.mf6.interfaces.iregridpackage import IRegridPackage from imod.msw.fixed_format import VariableMetaData from imod.msw.pkgbase import MetaSwapPackage from imod.msw.regrid.regrid_schemes import ScalingRegridMethod from imod.msw.utilities.common import concat_imod5 -from imod.typing import GridDataDict, Imod5DataDict +from imod.typing import Imod5DataDict from imod.typing.grid import ones_like @@ -88,7 +86,7 @@ def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "ScalingFactors": Import ScalingFactors from iMOD5 data. Pressure head factor is set to one for all factors, as well as all factors for urban areas. """ - cap_data = cast(GridDataDict, imod5_data["cap"]) + cap_data = imod5_data["cap"] grid_ones = ones_like(cap_data["boundary"]) data = {} diff --git a/imod/msw/sprinkling.py b/imod/msw/sprinkling.py index ea9925112..b5fa8cc69 100644 --- a/imod/msw/sprinkling.py +++ b/imod/msw/sprinkling.py @@ -1,4 +1,4 @@ -from typing import TextIO, cast +from typing import TextIO import numpy as np import pandas as pd @@ -180,9 +180,11 @@ def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "Sprinkling": (IDF), this grid defines in which layer a groundwater abstraction well should be placed. The ``"artificial_recharge"`` grid contains types which point to the type of abstraction: + * 0: no abstraction * 1: groundwater abstraction * 2: surfacewater abstraction + The ``"artificial_recharge_capacity"`` grid/constant defines the capacity of each groundwater or surfacewater abstraction. This is an ``1:1`` mapping: Each grid cell maps to a separate well. @@ -208,7 +210,7 @@ def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "Sprinkling": ------- Sprinkling package """ - cap_data = cast(GridDataDict, imod5_data["cap"]) + cap_data = imod5_data["cap"] if isinstance(cap_data["artificial_recharge_layer"], pd.DataFrame): data = _sprinkling_data_from_imod5_ipf(cap_data) else: diff --git a/imod/msw/utilities/imod5_converter.py b/imod/msw/utilities/imod5_converter.py index 70f88c982..0c4d21604 100644 --- a/imod/msw/utilities/imod5_converter.py +++ b/imod/msw/utilities/imod5_converter.py @@ -1,6 +1,9 @@ +from xarray.core.utils import is_scalar + from imod.logging import LogLevel, logger from imod.mf6 import StructuredDiscretization from imod.msw.utilities.common import concat_imod5 +from imod.msw.utilities.mask import MaskValues, MetaSwapActive from imod.typing import GridDataArray, GridDataDict from imod.typing.grid import ones_like from imod.util.spatial import get_cell_area @@ -10,7 +13,7 @@ def get_cell_area_from_imod5_data( imod5_cap: GridDataDict, ) -> GridDataArray: # area's per type of svats - mf6_area = get_cell_area(imod5_cap["wetted_area"]) + mf6_area = get_cell_area(imod5_cap["boundary"]) wetted_area = imod5_cap["wetted_area"] urban_area = imod5_cap["urban_area"] rural_area = mf6_area - (wetted_area + urban_area) @@ -43,7 +46,7 @@ def get_landuse_from_imod5_data( rural_landuse = imod5_cap["landuse"] # Urban landuse = 18 urban_landuse = ones_like(rural_landuse) * 18 - return concat_imod5(rural_landuse, urban_landuse) + return concat_imod5(rural_landuse, urban_landuse).astype(int) def get_rootzone_depth_from_imod5_data( @@ -63,7 +66,7 @@ def is_msw_active_cell( target_dis: StructuredDiscretization, imod5_cap: GridDataDict, msw_area: GridDataArray, -) -> tuple[GridDataArray, GridDataArray]: +) -> MetaSwapActive: """ Return grid of cells that are active in the coupled computation, based on following criteria: @@ -84,4 +87,33 @@ def is_msw_active_cell( (imod5_cap["boundary"] == 1) & (msw_area > 0) & (mf6_top_active >= 1) ) active = subunit_active.any(dim="subunit") - return active, subunit_active + return MetaSwapActive(active, subunit_active) + + +def _is_equal_scalar_value(da, value): + """ + Helper function to guarantee that the check in ``has_active_scaling_factor`` + can shortcut after is_scalar returns False. + """ + return da.to_numpy()[()] == value + + +def has_active_scaling_factor(imod5_cap: GridDataDict): + """ + Check if scaling factor grids are active. Carefully checks if data is + provided as constant (scalar) and if it matches an inactivity value. The + function shortcuts if data is provided as constant. + """ + variable_inactive_mapping = { + "perched_water_table_level": MaskValues.default, + "soil_moisture_fraction": 1.0, + "conductivitiy_factor": 1.0, + } + scaling_factor_inactive = True + for var, inactive_value in variable_inactive_mapping.items(): + da = imod5_cap[var] + scaling_factor_inactive &= is_scalar(da) and _is_equal_scalar_value( + da, inactive_value + ) + + return not scaling_factor_inactive diff --git a/imod/msw/utilities/mask.py b/imod/msw/utilities/mask.py new file mode 100644 index 000000000..6a26c9b9d --- /dev/null +++ b/imod/msw/utilities/mask.py @@ -0,0 +1,60 @@ +import numbers +from dataclasses import dataclass + +from imod.msw.pkgbase import MetaSwapPackage +from imod.typing import GridDataArray, GridDataDict + + +@dataclass +class MetaSwapActive: + all: GridDataArray + per_subunit: GridDataArray + + +@dataclass +class MaskValues: + """Stores sentinel values for nodata. Most cases use -9999.0, but exceptions + can be added here.""" + + default = -9999.0 + integer = 0 + + +def mask_and_broadcast_cap_data( + cap_data: GridDataDict, msw_active: MetaSwapActive +) -> GridDataDict: + """ + Mask and broadcast cap data, always mask with "all" of MetaSwapActive. + """ + return { + key: _mask_spatial_var(grid, msw_active.all) for key, grid in cap_data.items() + } + + +def mask_and_broadcast_pkg_data( + package: type[MetaSwapPackage], grid_data: GridDataDict, msw_active: MetaSwapActive +) -> GridDataDict: + """ + Mask and broadcast grid data, carefully mask per subunit if variable needs + to contain subunits. + """ + + return { + key: ( + _mask_spatial_var(grid, msw_active.per_subunit) + if key in package._with_subunit + else _mask_spatial_var(grid, msw_active.all) + ) + for key, grid in grid_data.items() + } + + +def _mask_spatial_var(da: GridDataArray, active: GridDataArray) -> GridDataArray: + if issubclass(da.dtype.type, numbers.Integral): + return da.where(active, other=MaskValues.integer) + elif issubclass(da.dtype.type, numbers.Real): + return da.where(active) + else: + raise TypeError( + f"Expected dtype float or integer. Received instead: {da.dtype}" + ) diff --git a/imod/msw/utilities/parse.py b/imod/msw/utilities/parse.py new file mode 100644 index 000000000..efe7f4428 --- /dev/null +++ b/imod/msw/utilities/parse.py @@ -0,0 +1,38 @@ +from pathlib import Path +from typing import TypeAlias + +ScalarType: TypeAlias = str | float | int + + +def _try_parsing_string_to_number(s: str) -> ScalarType: + """ + Convert string to number: + + "1" -> 1 + "1.0" -> 1.0 + "a" -> "a" + """ + try: + return int(s) + except ValueError: + pass + try: + return float(s) + except ValueError: + pass + return s + + +def read_para_sim(file: Path | str) -> dict[str, ScalarType]: + with open(file, "r") as f: + lines = f.readlines() + out = {} + for line in lines: + # Strip comments starting with "!" and split keys from values at the + # equals sign. + key_values = line[0 : line.find("!")].split("=") + if len(key_values) > 1: + key = key_values[0].strip() + value = _try_parsing_string_to_number(key_values[1].strip()) + out[key] = value + return out diff --git a/imod/msw/utilities/select.py b/imod/msw/utilities/select.py new file mode 100644 index 000000000..abd28818f --- /dev/null +++ b/imod/msw/utilities/select.py @@ -0,0 +1,16 @@ +from imod.typing import Imod5DataDict, SelSettingsType + +_DROP_LAYER_KWARGS: SelSettingsType = { + "layer": 0, + "drop": True, + "missing_dims": "ignore", +} + + +def drop_layer_dim_cap_data(imod5_data: Imod5DataDict) -> Imod5DataDict: + cap_data = imod5_data["cap"] + return { + "cap": { + key: da.isel(**_DROP_LAYER_KWARGS).compute() for key, da in cap_data.items() + } + } diff --git a/imod/prepare/__init__.py b/imod/prepare/__init__.py index 03342c87e..41e16f46d 100644 --- a/imod/prepare/__init__.py +++ b/imod/prepare/__init__.py @@ -43,6 +43,8 @@ from imod.prepare.topsystem import ( ALLOCATION_OPTION, DISTRIBUTING_OPTION, + SimulationAllocationOptions, + SimulationDistributingOptions, allocate_drn_cells, allocate_ghb_cells, allocate_rch_cells, diff --git a/imod/tests/conftest.py b/imod/tests/conftest.py index 0677d4b4f..b0f400192 100644 --- a/imod/tests/conftest.py +++ b/imod/tests/conftest.py @@ -92,4 +92,6 @@ well_test_data_transient, ) from .fixtures.msw_fixture import fixed_format_parser, simple_2d_grid_with_subunits +from .fixtures.msw_imod5_cap_fixture import imod5_cap_data +from .fixtures.msw_meteo_fixture import meteo_grids from .fixtures.msw_model_fixture import coupled_mf6_model, coupled_mf6wel, msw_model diff --git a/imod/tests/fixtures/msw_imod5_cap_fixture.py b/imod/tests/fixtures/msw_imod5_cap_fixture.py new file mode 100644 index 000000000..c57925015 --- /dev/null +++ b/imod/tests/fixtures/msw_imod5_cap_fixture.py @@ -0,0 +1,237 @@ +import numpy as np +import pytest +import xarray as xr +from numpy import nan + +from imod.typing import GridDataDict + + +@pytest.fixture(scope="function") +def imod5_cap_data() -> GridDataDict: + x = [1.0, 2.0, 3.0] + y = [3.0, 2.0, 1.0] + layer = [1] + dx = 1.0 + dy = -1.0 + + da_kwargs = {} + da_kwargs["dims"] = ("layer", "y", "x") + da_kwargs["coords"] = {"layer": layer, "y": y, "x": x, "dx": dx, "dy": dy} + + imod5_data = {} + d = {} + + # fmt: off + d["boundary"] = xr.DataArray( + np.array([ + [ + [1, 1, 1], + [0, 0, 0], + [1, 1, 0], + ]], + dtype=int), + **da_kwargs + ) + d["landuse"] = xr.DataArray( + np.array([ + [ + [1, 2, 3], + [0, 0, 0], + [3, 2, 1], + ]], + dtype=int), + **da_kwargs + ) + d["rootzone_thickness"] = xr.DataArray( + np.array([ + [ + [0.1, 0.1, 0.1], + [nan, nan, nan], + [0.2, 0.2, 0.2], + ]] + ), + **da_kwargs + ) + d["soil_physical_unit"] = xr.DataArray( + np.array([ + [ + [1, 1, 1], + [0, 0, 0], + [2, 1, 1], + ]], + dtype=int), + **da_kwargs + ) + d["surface_elevation"] = xr.DataArray( + np.array([ + [ + [1.1, 1.2, 1.3], + [nan, nan, nan], + [1.4, 1.5, 1.6], + ]] + ), + **da_kwargs + ) + d["artificial_recharge"] = xr.DataArray( + np.array([ + [ + [0.2, 0.2, 0.2], + [nan, nan, nan], + [0.3, 0.3, 0.3], + ]] + ), + **da_kwargs + ) + d["artificial_recharge_layer"] = xr.DataArray( + np.array([ + [ + [1, 2, 3], + [0, 0, 0], + [3, 2, 1], + ]], + dtype=int), + **da_kwargs + ) + d["artificial_recharge_capacity"] = xr.DataArray( + np.array([ + [ + [0.4, 0.4, 0.4], + [nan, nan, nan], + [0.6, 0.6, 0.6], + ]] + ), + **da_kwargs + ) + d["wetted_area"] = xr.DataArray( + np.array([ + [ + [0.1, 0.1, 0.1], + [nan, nan, nan], + [0.2, 0.2, 0.2], + ]] + ), + **da_kwargs + ) + d["urban_area"] = xr.DataArray( + np.array([ + [ + [0.2, 0.2, 0.2], + [nan, nan, nan], + [0.3, 0.3, 0.3], + ]] + ), + **da_kwargs + ) + d["urban_ponding_depth"] = xr.DataArray( + np.array([ + [ + [2.2, 2.2, 2.2], + [nan, nan, nan], + [2.3, 2.3, 2.3], + ]] + ), + **da_kwargs + ) + d["rural_ponding_depth"] = xr.DataArray( + np.array([ + [ + [1.2, 1.2, 1.2], + [nan, nan, nan], + [1.3, 1.3, 1.3], + ]] + ), + **da_kwargs + ) + d["urban_runoff_resistance"] = xr.DataArray( + np.array([ + [ + [3.2, 3.2, 3.2], + [nan, nan, nan], + [3.3, 3.3, 3.3], + ]] + ), + **da_kwargs + ) + d["rural_runoff_resistance"] = xr.DataArray( + np.array([ + [ + [3.6, 3.6, 3.6], + [nan, nan, nan], + [3.7, 3.7, 3.7], + ]] + ), + **da_kwargs + ) + d["urban_runon_resistance"] = xr.DataArray( + np.array([ + [ + [5.2, 5.2, 5.2], + [nan, nan, nan], + [5.3, 5.3, 5.3], + ]] + ), + **da_kwargs + ) + d["rural_runon_resistance"] = xr.DataArray( + np.array([ + [ + [5.6, 5.6, 5.6], + [nan, nan, nan], + [5.7, 5.7, 5.7], + ]] + ), + **da_kwargs + ) + d["urban_infiltration_capacity"] = xr.DataArray( + np.array([ + [ + [10.2, 10.2, 10.2], + [nan, nan, nan], + [10.3, 10.3, 10.3], + ]] + ), + **da_kwargs + ) + d["rural_infiltration_capacity"] = xr.DataArray( + np.array([ + [ + [20.2, 20.2, 20.2], + [nan, nan, nan], + [20.3, 20.3, 20.3], + ]] + ), + **da_kwargs + ) + d["perched_water_table_level"]= xr.DataArray( + np.array([ + [ + [2.0, 2.0, 2.0], + [nan, nan, nan], + [2.0, 2.0, 2.0], + ]] + ), + **da_kwargs + ) + d["soil_moisture_fraction"]= xr.DataArray( + np.array([ + [ + [1.5, 1.5, 1.5], + [nan, nan, nan], + [1.5, 1.5, 1.5], + ]] + ), + **da_kwargs + ) + d["conductivitiy_factor"]= xr.DataArray( + np.array([ + [ + [2.5, 2.5, 2.5], + [nan, nan, nan], + [2.5, 2.5, 2.5], + ]] + ), + **da_kwargs + ) + # fmt: on + imod5_data["cap"] = d + return imod5_data diff --git a/imod/tests/fixtures/msw_meteo_fixture.py b/imod/tests/fixtures/msw_meteo_fixture.py new file mode 100644 index 000000000..ebfdee449 --- /dev/null +++ b/imod/tests/fixtures/msw_meteo_fixture.py @@ -0,0 +1,43 @@ +import numpy as np +import pandas as pd +import pytest +import xarray as xr +from numpy import nan + +from imod.typing import GridDataArray + + +@pytest.fixture(scope="function") +def meteo_grids() -> tuple[GridDataArray, GridDataArray]: + x = [1.0, 2.0, 3.0] + y = [1.0, 2.0, 3.0] + time = pd.date_range(start="2000-01-01", end="2000-01-02", freq="D") + + dx = 1.0 + dy = -1.0 + # fmt: off + precipitation = xr.DataArray( + np.array( + [ + [[1.0, 1.0, 1.0], + [nan, nan, nan], + [1.0, 1.0, 1.0]], + + [[2.0, 2.0, 1.0], + [nan, nan, nan], + [1.0, 2.0, 1.0]], + ] + ), + dims=("time", "y", "x"), + coords={"time": time, "y": y, "x": x, "dx": dx, "dy": dy} + ) + + evapotranspiration = xr.DataArray( + np.array( + [1.0, 3.0] + ), + dims=("time",), + coords={"time": time} + ) + # fmt: on + return precipitation, evapotranspiration diff --git a/imod/tests/test_msw/test_grid_data.py b/imod/tests/test_msw/test_grid_data.py index 0d0fbf0c3..a775a93a9 100644 --- a/imod/tests/test_msw/test_grid_data.py +++ b/imod/tests/test_msw/test_grid_data.py @@ -461,7 +461,7 @@ def test_from_imod5_data(grid_data_dict: dict[str, xr.DataArray]): # Dis only needed for idomain dis = StructuredDiscretization(top, top - 0.1, idomain, validate=False) - griddata = GridData.from_imod5_data(imod5_data, target_dis=dis) + griddata, _ = GridData.from_imod5_data(imod5_data, target_dis=dis) expected_rootzone_depth = cap_data["rootzone_thickness"] * 0.01 xr.testing.assert_allclose( expected_rootzone_depth, griddata["rootzone_depth"].sel(subunit=0, drop=True) diff --git a/imod/tests/test_msw/test_infiltration.py b/imod/tests/test_msw/test_infiltration.py index b87940a88..47990e3a6 100644 --- a/imod/tests/test_msw/test_infiltration.py +++ b/imod/tests/test_msw/test_infiltration.py @@ -283,11 +283,11 @@ def test_from_imod5_data(data_infiltration): imod5_data = {"cap": cap_data} actual_pkg = Infiltration.from_imod5_data(imod5_data) - ones_vars = ["bottom_resistance", "extra_storage_coefficient"] - expected = expected_pkg.dataset.drop_vars(ones_vars) - actual = actual_pkg.dataset.drop_vars(ones_vars) + hardcoded_vars = {"bottom_resistance": -9999.0, "extra_storage_coefficient": 1.0} + expected = expected_pkg.dataset.drop_vars(hardcoded_vars.keys()) + actual = actual_pkg.dataset.drop_vars(hardcoded_vars.keys()) xr.testing.assert_equal(actual, expected) - for var in ones_vars: - assert (actual_pkg.dataset[var] == 1.0).all() + for var, value in hardcoded_vars.items(): + assert (actual_pkg.dataset[var] == value).all() diff --git a/imod/tests/test_msw/test_meteo_grid.py b/imod/tests/test_msw/test_meteo_grid.py index b5108842e..66d7fe95e 100644 --- a/imod/tests/test_msw/test_meteo_grid.py +++ b/imod/tests/test_msw/test_meteo_grid.py @@ -7,54 +7,21 @@ import numpy as np import pandas as pd import pytest -import xarray as xr -from numpy import nan from numpy.testing import assert_equal from imod.mf6.utilities.regrid import RegridderWeightsCache from imod.msw import MeteoGrid, MeteoGridCopy -def setup_meteo_grid(): - x = [1.0, 2.0, 3.0] - y = [1.0, 2.0, 3.0] - time = pd.date_range(start="2000-01-01", end="2000-01-02", freq="D") - - dx = 1.0 - dy = -1.0 - # fmt: off - precipitation = xr.DataArray( - np.array( - [ - [[1.0, 1.0, 1.0], - [nan, nan, nan], - [1.0, 1.0, 1.0]], - - [[2.0, 2.0, 1.0], - [nan, nan, nan], - [1.0, 2.0, 1.0]], - ] - ), - dims=("time", "y", "x"), - coords={"time": time, "y": y, "x": x, "dx": dx, "dy": dy} - ) - - evapotranspiration = xr.DataArray( - np.array( - [1.0, 3.0] - ), - dims=("time",), - coords={"time": time} - ) - # fmt: on - - meteo_grid = MeteoGrid(precipitation, evapotranspiration) +def test_meteo_grid_init(meteo_grids): + meteo_grids = MeteoGrid(*meteo_grids) - return meteo_grid + assert meteo_grids.dataset["precipitation"].dims == ("time", "y", "x") + assert meteo_grids.dataset["evapotranspiration"].dims == ("time",) -def test_meteo_grid(): - meteo_grid = setup_meteo_grid() +def test_meteo_grid_write(meteo_grids): + meteo_grid = MeteoGrid(*meteo_grids) with tempfile.TemporaryDirectory() as output_dir: output_dir = Path(output_dir) @@ -83,29 +50,8 @@ def test_meteo_grid(): assert gridnames == expected_filenames -def test_meteo_no_time_grid(): - x = [1.0, 2.0, 3.0] - y = [1.0, 2.0, 3.0] - time = pd.date_range(start="2000-01-01", end="2000-01-02", freq="D") - - dx = 1.0 - dy = -1.0 - # fmt: off - precipitation = xr.DataArray( - np.array( - [ - [[1.0, 1.0, 1.0], - [nan, nan, nan], - [1.0, 1.0, 1.0]], - - [[2.0, 2.0, 1.0], - [nan, nan, nan], - [1.0, 2.0, 1.0]], - ] - ), - dims=("time", "y", "x"), - coords={"time": time, "y": y, "x": x, "dx": dx, "dy": dy} - ) +def test_meteo_no_time_grid(meteo_grids): + precipitation, _ = meteo_grids evapotranspiration = 3.0 # fmt: on @@ -114,8 +60,8 @@ def test_meteo_no_time_grid(): MeteoGrid(precipitation, evapotranspiration) -def test_regrid_meteo(simple_2d_grid_with_subunits): - meteo = setup_meteo_grid() +def test_regrid_meteo(meteo_grids, simple_2d_grid_with_subunits): + meteo = MeteoGrid(*meteo_grids) new_grid = simple_2d_grid_with_subunits regrid_context = RegridderWeightsCache() @@ -126,44 +72,40 @@ def test_regrid_meteo(simple_2d_grid_with_subunits): assert np.all(regridded_ponding.dataset["y"].values == new_grid["y"].values) -def test_meteogridcopy_write(): - # Arrange - meteo_grid = setup_meteo_grid() - - with tempfile.TemporaryDirectory() as output_dir: - grid_dir = Path(output_dir) / "grid" - grid_dir.mkdir(exist_ok=True, parents=True) - meteo_grid.write(grid_dir) +def setup_written_meteo_grids(meteo_grids, tmp_path) -> Path: + meteo_grid = MeteoGrid(*meteo_grids) + grid_dir = tmp_path / "grid" + grid_dir.mkdir(exist_ok=True, parents=True) + meteo_grid.write(grid_dir) + return grid_dir - meteo_grid_copy = MeteoGridCopy(grid_dir / "mete_grid.inp") - copy_dir = Path(output_dir) / "copied" - copy_dir.mkdir(exist_ok=True, parents=True) - # Act - meteo_grid_copy.write(copy_dir) - # Assert - assert filecmp.cmp(grid_dir / "mete_grid.inp", copy_dir / "mete_grid.inp") +def test_meteogridcopy_write(meteo_grids, tmp_path): + # Arrange + grid_dir = setup_written_meteo_grids(meteo_grids, tmp_path) + meteo_grid_copy = MeteoGridCopy(grid_dir / "mete_grid.inp") + copy_dir = tmp_path / "copied" + copy_dir.mkdir(exist_ok=True, parents=True) + # Act + meteo_grid_copy.write(copy_dir) + # Assert + assert filecmp.cmp(grid_dir / "mete_grid.inp", copy_dir / "mete_grid.inp") -def test_meteogridcopy_from_imod5(): - meteo_grid = setup_meteo_grid() - with tempfile.TemporaryDirectory() as output_dir: - grid_dir = Path(output_dir) / "grid" - grid_dir.mkdir(exist_ok=True, parents=True) - meteo_grid.write(grid_dir) - - imod5_data = {} - imod5_data["extra"] = {} - imod5_data["extra"]["paths"] = [ - ["foo"], - [(grid_dir / "mete_grid.inp").resolve()], - ["bar"], - ] - - meteo_grid_copy = MeteoGridCopy.from_imod5_data(imod5_data) - copy_dir = Path(output_dir) / "copied" - copy_dir.mkdir(exist_ok=True, parents=True) - # Act - meteo_grid_copy.write(copy_dir) - # Assert - assert filecmp.cmp(grid_dir / "mete_grid.inp", copy_dir / "mete_grid.inp") +def test_meteogridcopy_from_imod5(meteo_grids, tmp_path): + # Arrange + grid_dir = setup_written_meteo_grids(meteo_grids, tmp_path) + imod5_data = {} + imod5_data["extra"] = {} + imod5_data["extra"]["paths"] = [ + ["foo"], + [(grid_dir / "mete_grid.inp").resolve()], + ["bar"], + ] + meteo_grid_copy = MeteoGridCopy.from_imod5_data(imod5_data) + copy_dir = tmp_path / "copied" + copy_dir.mkdir(exist_ok=True, parents=True) + # Act + meteo_grid_copy.write(copy_dir) + # Assert + assert filecmp.cmp(grid_dir / "mete_grid.inp", copy_dir / "mete_grid.inp") diff --git a/imod/tests/test_msw/test_model.py b/imod/tests/test_msw/test_model.py index 0e97a59c5..f0904e96f 100644 --- a/imod/tests/test_msw/test_model.py +++ b/imod/tests/test_msw/test_model.py @@ -1,11 +1,15 @@ +from copy import copy from pathlib import Path import pytest import xarray as xr from numpy.testing import assert_almost_equal, assert_equal +from pytest_cases import parametrize_with_cases from imod import mf6, msw from imod.mf6.utilities.regrid import RegridderWeightsCache +from imod.msw.model import DEFAULT_SETTINGS +from imod.typing import GridDataArray, Imod5DataDict def test_msw_model_write(msw_model, coupled_mf6_model, coupled_mf6wel, tmp_path): @@ -138,3 +142,156 @@ def test_coupled_model_regrid(msw_model, coupled_mf6_model, tmp_path): ) regridded_msw_model.write(tmp_path, mf6_discretization, regridded_mf6_wel) + + +def setup_written_meteo_grids( + meteo_grids: tuple[GridDataArray], tmp_path: Path +) -> Path: + precipitation, _ = meteo_grids + meteo_grid = msw.MeteoGrid(precipitation, precipitation) + grid_dir = tmp_path / "grid" + grid_dir.mkdir(exist_ok=True, parents=True) + meteo_grid.write(grid_dir) + return grid_dir + + +def setup_parasim_inp(directory: Path): + settings = copy(DEFAULT_SETTINGS) + settings["iybg"] = 2000 + settings["tdbg"] = 200.45 + settings["unsa_svat_path"] = str(directory) + + filename = directory / "para_sim.inp" + with open(filename, "w") as f: + rendered = msw.MetaSwapModel._template.render(settings=settings) + f.write(rendered) + return filename + + +def write_test_files(directory: Path, filenames: list[str]) -> list[Path]: + paths = [directory / filename for filename in filenames] + for p in paths: + with open(p, mode="w") as f: + f.write("test") + return paths + + +def setup_extra_files(meteo_grids: tuple[GridDataArray], directory: Path): + grid_dir = setup_written_meteo_grids(meteo_grids, directory) + setup_parasim_inp(grid_dir) + write_test_files(grid_dir, ["a.inp", "b.inp"]) + return { + "paths": [ + [str(grid_dir / fn)] + for fn in ["a.inp", "mete_grid.inp", "para_sim.inp", "b.inp"] + ] + } + + +class Imod5DataCases: + def case_grid(self, imod5_cap_data: Imod5DataDict) -> tuple[Imod5DataDict, bool]: + has_scaling_factor = True + return imod5_cap_data, has_scaling_factor + + def case_no_scaling_factors( + self, imod5_cap_data: Imod5DataDict + ) -> tuple[Imod5DataDict, bool]: + has_scaling_factor = False + cap_data = imod5_cap_data["cap"] + # open_projectfile_data adds layer kwargs to constants + layer_kwargs = {"coords": {"layer": [1]}, "dims": ("layer",)} + cap_data["perched_water_table_level"] = xr.DataArray([-9999.0], **layer_kwargs) + cap_data["soil_moisture_fraction"] = xr.DataArray([1.0], **layer_kwargs) + cap_data["conductivitiy_factor"] = xr.DataArray([1.0], **layer_kwargs) + return imod5_cap_data, has_scaling_factor + + def case_constants( + self, imod5_cap_data: Imod5DataDict + ) -> tuple[Imod5DataDict, bool]: + has_scaling_factor = False + cap_data = imod5_cap_data["cap"] + # open_projectfile_data adds layer kwargs to constants + layer_kwargs = {"coords": {"layer": [1]}, "dims": ("layer",)} + cap_data["perched_water_table_level"] = xr.DataArray([-9999.0], **layer_kwargs) + cap_data["soil_moisture_fraction"] = xr.DataArray([1.0], **layer_kwargs) + cap_data["conductivitiy_factor"] = xr.DataArray([1.0], **layer_kwargs) + cap_data["urban_ponding_depth"] = xr.DataArray([1.0], **layer_kwargs) + cap_data["rural_ponding_depth"] = xr.DataArray([1.0], **layer_kwargs) + cap_data["urban_runoff_resistance"] = xr.DataArray([1.0], **layer_kwargs) + cap_data["rural_runoff_resistance"] = xr.DataArray([1.0], **layer_kwargs) + cap_data["urban_runon_resistance"] = xr.DataArray([1.0], **layer_kwargs) + cap_data["rural_runon_resistance"] = xr.DataArray([1.0], **layer_kwargs) + cap_data["urban_infiltration_capacity"] = xr.DataArray([1.0], **layer_kwargs) + cap_data["rural_infiltration_capacity"] = xr.DataArray([1.0], **layer_kwargs) + return imod5_cap_data, has_scaling_factor + + +@parametrize_with_cases("imod5_data, has_scaling_factor", cases=Imod5DataCases) +def test_import_from_imod5( + imod5_data: Imod5DataDict, + has_scaling_factor: bool, + meteo_grids: tuple[GridDataArray], + coupled_mf6_model: mf6.Modflow6Simulation, + tmp_path: Path, +): + # Arrange + imod5_data["extra"] = setup_extra_files(meteo_grids, tmp_path) + times = coupled_mf6_model["time_discretization"].dataset.coords["time"] + dis_pkg = coupled_mf6_model["GWF_1"]["dis"] + # Act + model = msw.MetaSwapModel.from_imod5_data(imod5_data, dis_pkg, times) + # Assert + grid_packages = { + "grid", + "infiltration", + "ponding", + "sprinkling", + "idf_mapping", + } + expected_keys = grid_packages | { + "meteo_grid", + "prec_mapping", + "evt_mapping", + "coupling", + "extra_files", + "time_oc", + } + missing_keys = expected_keys - set(model.keys()) + assert not missing_keys + for pkgname in grid_packages: + # Test if all grid packages broadcasted to grid. + missing_dims = {"y", "x"} - set(model[pkgname].dataset.dims.keys()) + assert not missing_dims + assert "time" in model["time_oc"].dataset.dims.keys() + assert len(model["meteo_grid"].dataset.dims) == 0 + assert ("scaling_factor" in model.keys()) == has_scaling_factor + + +@parametrize_with_cases("imod5_data, has_scaling_factor", cases=Imod5DataCases) +def test_import_from_imod5_and_write( + imod5_data: Imod5DataDict, + has_scaling_factor: bool, + meteo_grids: tuple[GridDataArray], + coupled_mf6_model: mf6.Modflow6Simulation, + tmp_path: Path, +): + # Arrange + imod5_data["extra"] = setup_extra_files(meteo_grids, tmp_path) + times = coupled_mf6_model["time_discretization"].dataset.coords["time"] + dis_pkg = coupled_mf6_model["GWF_1"]["dis"] + npf_pkg = coupled_mf6_model["GWF_1"]["npf"] + active = dis_pkg["idomain"] == 1 + modeldir = tmp_path / "modeldir" + # Act + model = msw.MetaSwapModel.from_imod5_data(imod5_data, dis_pkg, times) + well_pkg = mf6.LayeredWell.from_imod5_cap_data(imod5_data) + mf6_wel_pkg = well_pkg.to_mf6_pkg( + active, dis_pkg["top"], dis_pkg["bottom"], npf_pkg["k"] + ) + model.write(modeldir, dis_pkg, mf6_wel_pkg, validate=False) + + # Assert + expected_n_files = 13 + if has_scaling_factor: + expected_n_files += 1 + assert len(list(modeldir.rglob(r"*.inp"))) == expected_n_files diff --git a/imod/tests/test_msw/test_utilities/test_mask.py b/imod/tests/test_msw/test_utilities/test_mask.py new file mode 100644 index 000000000..424f7dbc7 --- /dev/null +++ b/imod/tests/test_msw/test_utilities/test_mask.py @@ -0,0 +1,131 @@ +from copy import deepcopy +from typing import Callable + +import numpy as np +import pytest +import xarray as xr +from pytest_cases import parametrize_with_cases + +from imod.msw import GridData +from imod.msw.utilities.mask import ( + MetaSwapActive, + mask_and_broadcast_cap_data, + mask_and_broadcast_pkg_data, +) +from imod.typing import GridDataDict + + +@pytest.fixture(scope="function") +def coords_planar() -> dict: + x = [1.0, 2.0, 3.0] + y = [3.0, 2.0, 1.0] + dx = 1.0 + dy = 1.0 + return {"y": y, "x": x, "dx": dx, "dy": dy} + + +@pytest.fixture(scope="function") +def coords_subunit(coords_planar: dict) -> dict: + coords_subunit = deepcopy(coords_planar) + coords_subunit["subunit"] = [0, 1] + return coords_subunit + + +@pytest.fixture(scope="function") +def mask_fixture(coords_subunit: dict) -> MetaSwapActive: + mask_per_subunit = xr.DataArray( + np.array( + [ + [[0, 0, 0], [0, 1, 1], [0, 0, 0]], + [[1, 1, 1], [0, 1, 1], [0, 0, 0]], + ] + ).astype(bool), + dims=("subunit", "y", "x"), + coords=coords_subunit, + ) + mask_all = mask_per_subunit.any(dim="subunit") + + return MetaSwapActive(mask_all, mask_per_subunit) + + +@pytest.fixture(scope="function") +def data_subunit(coords_subunit: dict): + return xr.DataArray( + np.array( + [ + [[2.0, 2.0, 2.0], [2.0, 3.0, 3.0], [2.0, 2.0, 2.0]], + [[3.0, 3.0, 3.0], [2.0, 3.0, 3.0], [2.0, 2.0, 2.0]], + ] + ), + dims=("subunit", "y", "x"), + coords=coords_subunit, + ) + + +@pytest.fixture(scope="function") +def data_planar(data_subunit: xr.DataArray) -> xr.DataArray: + return data_subunit.sel(subunit=0, drop=True) + + +def integer_nan(da): + return da == 0 + + +def case_grid_float( + data_subunit: xr.DataArray, data_planar: xr.DataArray +) -> tuple[GridDataDict, Callable]: + return {"landuse": data_subunit, "elevation": data_planar}, np.isnan + + +def case_grid_integer( + data_subunit: xr.DataArray, data_planar: xr.DataArray +) -> tuple[GridDataDict, Callable]: + return { + "landuse": data_subunit.astype(int), + "elevation": data_planar.astype(int), + }, integer_nan + + +def case_constant_float(): + data = xr.DataArray(2.0) + return {"landuse": data, "elevation": data}, np.isnan + + +def case_constant_integer(): + data = xr.DataArray(2) + return {"landuse": data, "elevation": data}, integer_nan + + +@parametrize_with_cases("imod5_grid_data, isnan", cases=".") +def test_mask_and_broadcast_pkg_data( + imod5_grid_data: GridDataDict, isnan: Callable, mask_fixture: MetaSwapActive +): + # Act + masked_data = mask_and_broadcast_pkg_data(GridData, imod5_grid_data, mask_fixture) + + # Assert + assert isnan(masked_data["landuse"]).to_numpy().sum() == 11 + assert isnan(masked_data["elevation"]).to_numpy().sum() == 4 + + +@parametrize_with_cases("imod5_grid_data, isnan", cases=".") +def test_mask_and_broadcast_cap_data( + imod5_grid_data: GridDataDict, + isnan: Callable, + mask_fixture: MetaSwapActive, + coords_planar: dict, +): + # Arrange, cap data doesn't have subunit coords + imod5_grid_data["landuse"] = imod5_grid_data["landuse"].isel( + subunit=0, drop=True, missing_dims="ignore" + ) + + # Act + masked_data = mask_and_broadcast_cap_data(imod5_grid_data, mask_fixture) + + # Assert + coords_expected = xr.Coordinates(coords_planar) + xr.testing.assert_equal(masked_data["landuse"].coords, coords_expected) + xr.testing.assert_equal(masked_data["elevation"].coords, coords_expected) + assert isnan(masked_data["landuse"]).to_numpy().sum() == 4 + assert isnan(masked_data["elevation"]).to_numpy().sum() == 4 diff --git a/imod/tests/test_msw/test_utilities/test_parse.py b/imod/tests/test_msw/test_utilities/test_parse.py new file mode 100644 index 000000000..a4293252a --- /dev/null +++ b/imod/tests/test_msw/test_utilities/test_parse.py @@ -0,0 +1,28 @@ +from pytest_cases import parametrize_with_cases + +from imod.msw.utilities.parse import _try_parsing_string_to_number + + +class ParseCases: + def case_int(self): + return "1", int + + def case_float(self): + return "1.0", float + + def case_string(self): + return "a", str + + def case_aterisk(self): + return "*", str + + def case_exclamation(self): + return "!", str + + +@parametrize_with_cases(["s", "expected_type"], cases=ParseCases) +def test_try_parsing_string_to_value(s, expected_type): + # Act + parsed = _try_parsing_string_to_number(s) + # Assert + assert type(parsed) is expected_type diff --git a/imod/typing/__init__.py b/imod/typing/__init__.py index 0217f9782..da8bb6c19 100644 --- a/imod/typing/__init__.py +++ b/imod/typing/__init__.py @@ -12,7 +12,6 @@ GridDataArray: TypeAlias = Union[xr.DataArray, xu.UgridDataArray] GridDataset: TypeAlias = Union[xr.Dataset, xu.UgridDataset] GridDataDict: TypeAlias = dict[str, GridDataArray] -Imod5DataDict: TypeAlias = dict[str, GridDataDict | dict[str, list[str]]] ScalarAsDataArray: TypeAlias = Union[xr.DataArray, xu.UgridDataArray] ScalarAsDataset: TypeAlias = Union[xr.Dataset, xu.UgridDataset] UnstructuredData: TypeAlias = Union[xu.UgridDataset, xu.UgridDataArray] @@ -26,6 +25,11 @@ class SelSettingsType(TypedDict, total=False): missing_dims: Literal["raise", "warn", "ignore"] +class Imod5DataDict(TypedDict, total=False): + cap: GridDataDict + extra: dict[str, list[str]] + + # Types for optional dependencies. if TYPE_CHECKING: import geopandas as gpd