diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst index a69de959a..e71ba80b4 100644 --- a/docs/api/changelog.rst +++ b/docs/api/changelog.rst @@ -64,6 +64,10 @@ Fixed - :meth:`imod.mf6.NodePropertyFlow.regrid_like` now regrids ``k33`` using the correct method, namely ``mean`` instead of ``harmonic_mean``. As this is the appropriate method for horizontal regridding of ``k33``. +- :meth:`imod.msw.MetaSwapModel.from_imod5_data`, + :meth:`imod.mf6.Recharge.from_imod5_cap_data`, + :meth:`imod.mf6.LayeredWell.from_imod5_cap_data` now regrids the iMOD5 CAP + data to the MODFLOW6 target discretization. Changed ~~~~~~~ diff --git a/imod/common/utilities/regrid.py b/imod/common/utilities/regrid.py index 4fbea0a31..bffb13cb4 100644 --- a/imod/common/utilities/regrid.py +++ b/imod/common/utilities/regrid.py @@ -20,13 +20,14 @@ from imod.common.utilities.dataclass_type import DataclassType, EmptyRegridMethod from imod.common.utilities.dtype import is_integer from imod.common.utilities.value_filters import is_valid +from imod.typing import Imod5DataDict from imod.typing.grid import ( GridDataArray, GridDataset, is_unstructured, ones_like, ) -from imod.util.dims import enforced_dim_order +from imod.util.dims import drop_layer_dim_cap_data, enforced_dim_order from imod.util.regrid import ( RegridderType, RegridderWeightsCache, @@ -77,7 +78,7 @@ def _regrid_array( # the dataarray might be a scalar. If it is, then it does not need regridding. if scalar_da: - return da.values[()] + return da.squeeze(drop=True) if isinstance(da, xr.DataArray): coords = da.coords @@ -459,3 +460,32 @@ def _get_regridding_domain( ).astype(int) return new_idomain + + +def regrid_imod5_cap_data( + imod5_data: Imod5DataDict, + target_dis: IRegridPackage, + regridder_types: DataclassType, + regrid_cache: RegridderWeightsCache, +) -> Imod5DataDict: + """ + Regrid iMOD5 CAP data to consistent grid. This is necessary to be able to + use iMOD5 data in MetaSWAP, as the grid of the iMOD5 CAP data is not + necessarily the same as the grid of the target MODFLOW6 discretization. The + regridding process ensures consistency between the iMOD5 CAP data and the + target MODFLOW6 grid. Used in both ``imod.msw.MetaSwapModel.from_imod5_data`` + and ``imod.mf6.Recharge.from_imod5_cap_data``. + """ + # Drop layer coords + imod5_cap_no_layer = drop_layer_dim_cap_data(imod5_data) + target_grid = target_dis.dataset["idomain"].isel(layer=0, drop=True) + # Regrid the input data + cap_data_regridded = _regrid_package_data( + imod5_cap_no_layer["cap"], target_grid, regridder_types, regrid_cache + ) + extra_paths = imod5_data["extra"]["paths"] + imod5_regridded: Imod5DataDict = { + "cap": cap_data_regridded, + "extra": {"paths": extra_paths}, + } + return imod5_regridded diff --git a/imod/formats/prj/prj.py b/imod/formats/prj/prj.py index 56e582347..dd9859a92 100644 --- a/imod/formats/prj/prj.py +++ b/imod/formats/prj/prj.py @@ -98,7 +98,7 @@ "rural_infiltration_capacity", "perched_water_table_level", "soil_moisture_fraction", - "conductivitiy_factor", + "conductivity_factor", "plot_number", "steering_location", "plot_drainage_level", diff --git a/imod/mf6/model_gwf.py b/imod/mf6/model_gwf.py index 5e7b544af..f834aa1a7 100644 --- a/imod/mf6/model_gwf.py +++ b/imod/mf6/model_gwf.py @@ -45,7 +45,7 @@ SimulationDistributingOptions, ) from imod.schemata import TypeSchema -from imod.typing import GridDataArray, StressPeriodTimesType +from imod.typing import GridDataArray, Imod5DataDict, StressPeriodTimesType from imod.typing.grid import zeros_like from imod.util.regrid import RegridderWeightsCache @@ -354,8 +354,13 @@ def from_imod5_data( ) if "cap" in imod5_keys: - result["msw-sprinkling"] = LayeredWell.from_imod5_cap_data(imod5_data) # type: ignore - result["msw-rch"] = Recharge.from_imod5_cap_data(imod5_data, dis_pkg) # type: ignore + imod5_cap_data_dict = cast(Imod5DataDict, imod5_data) + result["msw-sprinkling"] = LayeredWell.from_imod5_cap_data( + imod5_cap_data_dict, dis_pkg, regrid_cache=regrid_cache + ) + result["msw-rch"] = Recharge.from_imod5_cap_data( + imod5_cap_data_dict, dis_pkg, regrid_cache=regrid_cache + ) # import ghb's ghb_keys = [key for key in imod5_keys if key[0:3] == "ghb"] diff --git a/imod/mf6/rch.py b/imod/mf6/rch.py index 642c7b637..1933f480e 100644 --- a/imod/mf6/rch.py +++ b/imod/mf6/rch.py @@ -6,10 +6,14 @@ from imod.common.interfaces.iregridpackage import IRegridPackage from imod.common.utilities.dataclass_type import DataclassType +from imod.common.utilities.regrid import regrid_imod5_cap_data from imod.logging import init_log_decorator from imod.mf6.aggregate.aggregate_schemes import RechargeAggregationMethod from imod.mf6.dis import StructuredDiscretization, VerticesDiscretization -from imod.mf6.regrid.regrid_schemes import RechargeRegridMethod +from imod.mf6.regrid.regrid_schemes import ( + CapDataRechargeRegridMethod, + RechargeRegridMethod, +) from imod.mf6.topsystem import TopSystemBoundaryCondition from imod.mf6.utilities.imod5_converter import ( convert_unit_rch_rate, @@ -39,7 +43,6 @@ enforce_dim_order, is_planar_grid, ) -from imod.util.dims import drop_layer_dim_cap_data from imod.util.regrid import RegridderWeightsCache @@ -286,6 +289,8 @@ def from_imod5_cap_data( cls, imod5_data: Imod5DataDict, target_dis: StructuredDiscretization, + regridder_types: CapDataRechargeRegridMethod = CapDataRechargeRegridMethod(), + regrid_cache: RegridderWeightsCache = RegridderWeightsCache(), ) -> "Recharge": """ Construct an rch-package from iMOD5 data in the CAP package, loaded with @@ -293,7 +298,9 @@ def from_imod5_cap_data( used to couple MODFLOW6 to MetaSWAP models. Active cells will have a recharge rate of 0.0. """ - cap_data = drop_layer_dim_cap_data(imod5_data)["cap"] + cap_data = regrid_imod5_cap_data( + imod5_data, target_dis, regridder_types, regrid_cache + )["cap"] msw_area = get_cell_area_from_imod5_data(cap_data) msw_active = is_msw_active_cell(target_dis, cap_data, msw_area) diff --git a/imod/mf6/regrid/regrid_schemes.py b/imod/mf6/regrid/regrid_schemes.py index 0af7267af..f2a957316 100644 --- a/imod/mf6/regrid/regrid_schemes.py +++ b/imod/mf6/regrid/regrid_schemes.py @@ -450,3 +450,28 @@ class StorageCoefficientRegridMethod(DataclassType): convertible: RegridVarType = (RegridderType.OVERLAP, "mode") storage_coefficient: RegridVarType = (RegridderType.OVERLAP, "mean") specific_yield: RegridVarType = (RegridderType.OVERLAP, "mean") + + +@dataclass(config=_CONFIG) +class CapDataRechargeRegridMethod(DataclassType): + """ + Object containing regridder methods for CAP data for the + :class:`imod.mf6.Recharge.from_imod5_cap_data` method. This contains regridder + methods for only the relevant CAP variables for the recharge package. + """ + + boundary: RegridVarType = (RegridderType.OVERLAP, "mode") + wetted_area: RegridVarType = (RegridderType.RELATIVEOVERLAP, "conductance") + urban_area: RegridVarType = (RegridderType.RELATIVEOVERLAP, "conductance") + + +@dataclass(config=_CONFIG) +class CapDataWellRegridMethod(DataclassType): + """ + Object containing regridder methods for CAP data for the + :class:`imod.mf6.LayeredWell.from_imod5_cap_data` method. This contains + regridder methods for only the relevant CAP variables for the well package. + """ + + artificial_recharge: RegridVarType = (RegridderType.OVERLAP, "mean") + artificial_recharge_layer: RegridVarType = (RegridderType.OVERLAP, "mode") diff --git a/imod/mf6/riv.py b/imod/mf6/riv.py index 15f82e2e3..0957ebbad 100644 --- a/imod/mf6/riv.py +++ b/imod/mf6/riv.py @@ -141,7 +141,7 @@ def _create_drain_from_leftover_riv_imod5_data( else: drain_leftover_data = infiltration_drn_data - return Drainage(**drain_leftover_data) # type: ignore[arg-type] + return Drainage(**drain_leftover_data) # type: ignore class River(TopSystemBoundaryCondition, IRegridPackage): diff --git a/imod/mf6/utilities/imod5_converter.py b/imod/mf6/utilities/imod5_converter.py index 5044fd79f..adee28a6e 100644 --- a/imod/mf6/utilities/imod5_converter.py +++ b/imod/mf6/utilities/imod5_converter.py @@ -4,12 +4,12 @@ import pandas as pd import xarray as xr +from imod.common.interfaces.iregridpackage import IRegridPackage from imod.common.utilities.dataclass_type import DataclassType -from imod.common.utilities.regrid import _regrid_package_data +from imod.common.utilities.regrid import _regrid_package_data, regrid_imod5_cap_data from imod.mf6.package import Package from imod.typing import GridDataArray, GridDataDict, Imod5DataDict from imod.typing.grid import full_like -from imod.util.dims import drop_layer_dim_cap_data from imod.util.regrid import RegridderWeightsCache @@ -83,7 +83,12 @@ def _well_from_imod5_cap_grid_data(cap_data: GridDataDict) -> dict[str, np.ndarr return data -def well_from_imod5_cap_data(imod5_data: Imod5DataDict) -> dict[str, np.ndarray]: +def well_from_imod5_cap_data( + imod5_data: Imod5DataDict, + target_dis: IRegridPackage, + regridder_types: DataclassType, + regrid_cache: RegridderWeightsCache, +) -> 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 @@ -110,12 +115,16 @@ def well_from_imod5_cap_data(imod5_data: Imod5DataDict) -> dict[str, np.ndarray] capacity is already defined in the point data. This is an ``n:1`` mapping: multiple grid cells can map to one well. """ - cap_data = drop_layer_dim_cap_data(imod5_data)["cap"] + cap_data = imod5_data["cap"] + has_ipf_well = isinstance(cap_data["artificial_recharge_layer"], pd.DataFrame) - if isinstance(cap_data["artificial_recharge_layer"], pd.DataFrame): + if has_ipf_well: return _well_from_imod5_cap_point_data(cap_data) else: - return _well_from_imod5_cap_grid_data(cap_data) + cap_data_regridded = regrid_imod5_cap_data( + imod5_data, target_dis, regridder_types, regrid_cache + )["cap"] + return _well_from_imod5_cap_grid_data(cap_data_regridded) def regrid_imod5_pkg_data( diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py index ab4a62fa8..83eca04bd 100644 --- a/imod/mf6/wel.py +++ b/imod/mf6/wel.py @@ -28,6 +28,7 @@ ) from imod.mf6.mf6_wel_adapter import Mf6Wel, concat_indices_to_cellid from imod.mf6.package import Package +from imod.mf6.regrid.regrid_schemes import CapDataWellRegridMethod from imod.mf6.utilities.dataset import remove_inactive from imod.mf6.utilities.imod5_converter import well_from_imod5_cap_data from imod.mf6.validation_settings import ValidationSettings @@ -45,6 +46,7 @@ from imod.typing import GridDataArray, Imod5DataDict, StressPeriodTimesType from imod.typing.grid import is_spatial_grid, ones_like from imod.util.expand_repetitions import average_timeseries, resample_timeseries +from imod.util.regrid import RegridderWeightsCache from imod.util.structured import values_within_range ABSTRACT_METH_ERROR_MSG = "Method in abstract base class called" @@ -1443,7 +1445,13 @@ def _validate_imod5_depth_information( raise ValueError(log_msg) @classmethod - def from_imod5_cap_data(cls, imod5_data: Imod5DataDict): + def from_imod5_cap_data( + cls, + imod5_data: Imod5DataDict, + target_dis: StructuredDiscretization, + regridder_types: CapDataWellRegridMethod = CapDataWellRegridMethod(), + regrid_cache: RegridderWeightsCache = RegridderWeightsCache(), + ): """ Create LayeredWell from imod5_data in "cap" package. Abstraction data for sprinkling is defined in iMOD5 either with grids (IDF) or points @@ -1481,7 +1489,9 @@ def from_imod5_cap_data(cls, imod5_data: Imod5DataDict): belongs, as returned by :func:`imod.formats.prj.open_projectfile_data`. """ - data = well_from_imod5_cap_data(imod5_data) + data = well_from_imod5_cap_data( + imod5_data, target_dis, regridder_types, regrid_cache + ) return cls(**data) # type: ignore @standard_log_decorator() diff --git a/imod/msw/model.py b/imod/msw/model.py index eaf81fdbf..1c604c352 100644 --- a/imod/msw/model.py +++ b/imod/msw/model.py @@ -13,6 +13,7 @@ from imod.common.constants import MaskValues from imod.common.utilities.clip import clip_by_grid from imod.common.utilities.partitioninfo import create_partition_info +from imod.common.utilities.regrid import regrid_imod5_cap_data from imod.common.utilities.value_filters import enforce_scalar from imod.common.utilities.version import prepend_content_with_version_info from imod.mf6.dis import StructuredDiscretization @@ -38,16 +39,18 @@ from imod.msw.output_control import TimeOutputControl from imod.msw.pkgbase import MetaSwapPackage from imod.msw.ponding import Ponding +from imod.msw.regrid.regrid_schemes import CapDataRegridMethod 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.imod5_converter import ( + has_active_scaling_factor, +) from imod.msw.utilities.mask import mask_and_broadcast_cap_data from imod.msw.utilities.parse import read_para_sim from imod.msw.vegetation import AnnualCropFactors from imod.typing import GridDataArray, Imod5DataDict -from imod.util.dims import drop_layer_dim_cap_data from imod.util.regrid import RegridderWeightsCache from imod.util.time import to_datetime_internal @@ -605,7 +608,9 @@ def from_imod5_data( imod5_data: Imod5DataDict, target_dis: StructuredDiscretization, times: list[datetime], - ): + regridder_types: CapDataRegridMethod = CapDataRegridMethod(), + regrid_cache: RegridderWeightsCache = RegridderWeightsCache(), + ) -> "MetaSwapModel": """ Construct a MetaSWAP model from iMOD5 data in the CAP package, loaded with the :func:`imod.formats.prj.open_projectfile_data` function. @@ -616,29 +621,43 @@ def from_imod5_data( 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. + Target discretization, iMOD5 CAP data will be regridded to this + 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. + regridder_types: CapDataRegridMethod, default CapDataRegridMethod() + Custom regrid method for CAP data. + regrid_cache: RegridderWeightsCache, default RegridderWeightsCache() + Cache for regridder weights, can be used to speed up regridding if the + same regridders are used multiple times. Returns ------- MetaSwapModel MetaSWAP model imported from imod5 data. """ + # Path and settings management 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) + # Regrid iMOD5 CAP data to target discretization. + imod5_regridded = regrid_imod5_cap_data( + imod5_data, target_dis, regridder_types, regrid_cache + ) + # Test with regridded data instead of masked, as masking broadcasts + # scalars to grids, which causes the is_scalar check in + # has_active_scaling_factor to always return False. + active_scaling_factor = has_active_scaling_factor(imod5_regridded["cap"]) + # Setup model model = cls(unsa_svat_path, parasim_settings) model["grid"], msw_active = GridData.from_imod5_data( - imod5_cap_no_layer, target_dis + imod5_regridded, target_dis ) cap_data_masked = mask_and_broadcast_cap_data( - imod5_cap_no_layer["cap"], msw_active + imod5_regridded["cap"], msw_active ) imod5_masked: Imod5DataDict = { "cap": cap_data_masked, @@ -650,7 +669,7 @@ def from_imod5_data( 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"]): + if active_scaling_factor: model["scaling_factor"] = ScalingFactors.from_imod5_data(imod5_masked) area = model["grid"]["area"].isel(subunit=0, drop=True) model["idf_mapping"] = IdfMapping(area, MaskValues.msw_default) diff --git a/imod/msw/regrid/regrid_schemes.py b/imod/msw/regrid/regrid_schemes.py index 9caedc6dd..cf82018b1 100644 --- a/imod/msw/regrid/regrid_schemes.py +++ b/imod/msw/regrid/regrid_schemes.py @@ -10,6 +10,71 @@ ) +@dataclass(config=_CONFIG) +class CapDataRegridMethod(DataclassType): + """ + Object containing regridder methods for the iMOD5 CAP data. This can be + provided to the ``regrid_like`` method to regrid with custom settings. This + is used when constructing a MetaSWAP model from iMOD5 data. + + Parameters + ---------- + boundary: tuple, default (RegridderType.OVERLAP, "mode") + landuse: tuple, default (RegridderType.OVERLAP, "mode") + rootzone_thickness: tuple, default (RegridderType.OVERLAP, "mean") + soil_physical_unit: tuple, default (RegridderType.OVERLAP, "mode") + meteo_station_number: tuple, default (RegridderType.OVERLAP, "mode") + surface_elevation: tuple, default (RegridderType.OVERLAP, "mean") + artificial_recharge: tuple, default (RegridderType.OVERLAP, "mean") + artificial_recharge_layer: tuple, default (RegridderType.OVERLAP, "mode") + artificial_recharge_capacity: tuple, default (RegridderType.OVERLAP, "mean") + wetted_area: tuple, default (RegridderType.RELATIVEOVERLAP, "conductance") + urban_area: tuple, default (RegridderType.RELATIVEOVERLAP, "conductance") + urban_ponding_depth: tuple, default (RegridderType.OVERLAP, "mean") + rural_ponding_depth: tuple, default (RegridderType.OVERLAP, "mean") + urban_runoff_resistance: tuple, default (RegridderType.OVERLAP, "mean") + rural_runoff_resistance: tuple, default (RegridderType.OVERLAP, "mean") + urban_runon_resistance: tuple, default (RegridderType.OVERLAP, "mean") + rural_runon_resistance: tuple, default (RegridderType.OVERLAP, "mean") + urban_infiltration_capacity: tuple, default (RegridderType.OVERLAP, "mean") + rural_infiltration_capacity: tuple, default (RegridderType.OVERLAP, "mean") + perched_water_table_level: tuple, default (RegridderType.OVERLAP, "mean") + soil_moisture_fraction: tuple, default (RegridderType.OVERLAP, "mean") + conductivity_factor: tuple, default (RegridderType.OVERLAP, "mean") + plot_number: tuple, default (RegridderType.OVERLAP, "mode") + steering_location: tuple, default (RegridderType.OVERLAP, "mode") + plot_drainage_level: tuple, default (RegridderType.OVERLAP, "mean") + plot_drainage_resistance: tuple, default (RegridderType.OVERLAP, "mean") + """ + + boundary: RegridVarType = (RegridderType.OVERLAP, "mode") + landuse: RegridVarType = (RegridderType.OVERLAP, "mode") + rootzone_thickness: RegridVarType = (RegridderType.OVERLAP, "mean") + soil_physical_unit: RegridVarType = (RegridderType.OVERLAP, "mode") + meteo_station_number: RegridVarType = (RegridderType.OVERLAP, "mode") + surface_elevation: RegridVarType = (RegridderType.OVERLAP, "mean") + artificial_recharge: RegridVarType = (RegridderType.OVERLAP, "mean") + artificial_recharge_layer: RegridVarType = (RegridderType.OVERLAP, "mode") + artificial_recharge_capacity: RegridVarType = (RegridderType.OVERLAP, "mean") + wetted_area: RegridVarType = (RegridderType.RELATIVEOVERLAP, "conductance") + urban_area: RegridVarType = (RegridderType.RELATIVEOVERLAP, "conductance") + urban_ponding_depth: RegridVarType = (RegridderType.OVERLAP, "mean") + rural_ponding_depth: RegridVarType = (RegridderType.OVERLAP, "mean") + urban_runoff_resistance: RegridVarType = (RegridderType.OVERLAP, "mean") + rural_runoff_resistance: RegridVarType = (RegridderType.OVERLAP, "mean") + urban_runon_resistance: RegridVarType = (RegridderType.OVERLAP, "mean") + rural_runon_resistance: RegridVarType = (RegridderType.OVERLAP, "mean") + urban_infiltration_capacity: RegridVarType = (RegridderType.OVERLAP, "mean") + rural_infiltration_capacity: RegridVarType = (RegridderType.OVERLAP, "mean") + perched_water_table_level: RegridVarType = (RegridderType.OVERLAP, "mean") + soil_moisture_fraction: RegridVarType = (RegridderType.OVERLAP, "mean") + conductivity_factor: RegridVarType = (RegridderType.OVERLAP, "mean") + plot_number: RegridVarType = (RegridderType.OVERLAP, "mode") + steering_location: RegridVarType = (RegridderType.OVERLAP, "mode") + plot_drainage_level: RegridVarType = (RegridderType.OVERLAP, "mean") + plot_drainage_resistance: RegridVarType = (RegridderType.OVERLAP, "mean") + + @dataclass(config=_CONFIG) class SprinklingRegridMethod(DataclassType): """ diff --git a/imod/msw/scaling_factors.py b/imod/msw/scaling_factors.py index 85eb37324..e4dd10dae 100644 --- a/imod/msw/scaling_factors.py +++ b/imod/msw/scaling_factors.py @@ -109,7 +109,7 @@ def from_imod5_data(cls, imod5_data: Imod5DataDict) -> "ScalingFactors": cap_data["soil_moisture_fraction"], grid_ones ) data["scale_hydraulic_conductivity"] = concat_imod5( - cap_data["conductivitiy_factor"], grid_ones + cap_data["conductivity_factor"], grid_ones ) data["scale_pressure_head"] = concat_imod5(grid_ones, grid_ones) data["depth_perched_water_table"] = cap_data["perched_water_table_level"] diff --git a/imod/msw/utilities/imod5_converter.py b/imod/msw/utilities/imod5_converter.py index 58fcce04c..c6a3cd32d 100644 --- a/imod/msw/utilities/imod5_converter.py +++ b/imod/msw/utilities/imod5_converter.py @@ -106,7 +106,7 @@ def has_active_scaling_factor(imod5_cap: GridDataDict): variable_inactive_mapping = { "perched_water_table_level": MaskValues.msw_default, "soil_moisture_fraction": 1.0, - "conductivitiy_factor": 1.0, + "conductivity_factor": 1.0, } scaling_factor_inactive = True for var, inactive_value in variable_inactive_mapping.items(): diff --git a/imod/tests/conftest.py b/imod/tests/conftest.py index 9bc606476..d1f66a799 100644 --- a/imod/tests/conftest.py +++ b/imod/tests/conftest.py @@ -23,6 +23,8 @@ ) from .fixtures.flow_transport_simulation_fixture import flow_transport_simulation from .fixtures.imod5_cap_data import ( + cap_coupled_dis_grid, + cap_coupled_dis_grid__big, cap_data_sprinkling_grid, cap_data_sprinkling_grid__big, cap_data_sprinkling_points, diff --git a/imod/tests/fixtures/imod5_cap_data.py b/imod/tests/fixtures/imod5_cap_data.py index e74f1664c..88a2a5d7d 100644 --- a/imod/tests/fixtures/imod5_cap_data.py +++ b/imod/tests/fixtures/imod5_cap_data.py @@ -3,6 +3,7 @@ import pytest import xarray as xr +from imod.mf6 import StructuredDiscretization from imod.typing import Imod5DataDict @@ -24,6 +25,42 @@ def zeros_dask_grid(n): return da.chunk({"x": n, "y": n}) +@pytest.fixture(scope="function") +def cap_coupled_dis_grid() -> StructuredDiscretization: + n = 3 + layer_coord = np.arange(n) + 1 + layer_template = xr.DataArray( + layer_coord, dims="layer", coords={"layer": layer_coord} + ) + template_grid = layer_template * zeros_grid(n) + + idomain = template_grid.astype(int) + 1 + top = 10.0 + bottom = xr.DataArray( + [0.0, -10.0, -20.0], dims="layer", coords={"layer": layer_coord} + ) + + return StructuredDiscretization(idomain=idomain, top=top, bottom=bottom) + + +@pytest.fixture(scope="function") +def cap_coupled_dis_grid__big() -> StructuredDiscretization: + n = 1000 + layer_coord = np.arange(3) + 1 + layer_template = xr.DataArray( + layer_coord, dims="layer", coords={"layer": layer_coord} + ) + template_grid = layer_template * zeros_dask_grid(n) + + idomain = template_grid.astype(int) + 1 + top = 10.0 + bottom = xr.DataArray( + [0.0, -10.0, -20.0], dims="layer", coords={"layer": layer_coord} + ) + + return StructuredDiscretization(idomain=idomain, top=top, bottom=bottom) + + @pytest.fixture(scope="function") def cap_data_sprinkling_grid() -> Imod5DataDict: n = 3 @@ -46,7 +83,7 @@ def cap_data_sprinkling_grid() -> Imod5DataDict: "artificial_recharge_capacity": xr.DataArray(25.0), } - return {"cap": cap_data} + return {"cap": cap_data, "extra": {"paths": ["path1", "path2"]}} @pytest.fixture(scope="function") @@ -69,7 +106,7 @@ def cap_data_sprinkling_grid__big() -> Imod5DataDict: "artificial_recharge_capacity": xr.DataArray(25.0), } - return {"cap": cap_data} + return {"cap": cap_data, "extra": {"paths": ["path1", "path2"]}} @pytest.fixture(scope="function") @@ -94,7 +131,7 @@ def cap_data_sprinkling_grid_with_layer() -> Imod5DataDict: "artificial_recharge_capacity": xr.DataArray(25.0), } - return {"cap": cap_data} + return {"cap": cap_data, "extra": {"paths": ["path1", "path2"]}} @pytest.fixture(scope="function") @@ -126,4 +163,4 @@ def cap_data_sprinkling_points() -> Imod5DataDict: "artificial_recharge_capacity": xr.DataArray(25.0), } - return {"cap": cap_data} + return {"cap": cap_data, "extra": {"paths": ["path1", "path2"]}} diff --git a/imod/tests/fixtures/msw_imod5_cap_fixture.py b/imod/tests/fixtures/msw_imod5_cap_fixture.py index c57925015..28c53335b 100644 --- a/imod/tests/fixtures/msw_imod5_cap_fixture.py +++ b/imod/tests/fixtures/msw_imod5_cap_fixture.py @@ -222,7 +222,7 @@ def imod5_cap_data() -> GridDataDict: ), **da_kwargs ) - d["conductivitiy_factor"]= xr.DataArray( + d["conductivity_factor"]= xr.DataArray( np.array([ [ [2.5, 2.5, 2.5], diff --git a/imod/tests/test_mf6/test_mf6_rch.py b/imod/tests/test_mf6/test_mf6_rch.py index 16fef6bb9..5867234f0 100644 --- a/imod/tests/test_mf6/test_mf6_rch.py +++ b/imod/tests/test_mf6/test_mf6_rch.py @@ -578,6 +578,7 @@ def test_from_imod5_cap_data(imod5_dataset): # Arrange data = deepcopy(imod5_dataset[0]) target_discretization = StructuredDiscretization.from_imod5_data(data) + data["extra"] = {"paths": ["path1", "path2"]} data["cap"] = {} msw_bound = data["bnd"]["ibound"].isel(layer=0, drop=False) data["cap"]["boundary"] = msw_bound @@ -613,6 +614,7 @@ def test_from_imod5_cap_data__regrid(imod5_dataset): # Arrange data = deepcopy(imod5_dataset[0]) target_discretization = StructuredDiscretization.from_imod5_data(data) + data["extra"] = {"paths": ["path1", "path2"]} data["cap"] = {} msw_bound = data["bnd"]["ibound"].isel(layer=0, drop=False) data["cap"]["boundary"] = msw_bound @@ -637,6 +639,7 @@ def test_from_imod5_cap_data__clip_box(imod5_dataset): # Arrange data = deepcopy(imod5_dataset[0]) target_discretization = StructuredDiscretization.from_imod5_data(data) + data["extra"] = {"paths": ["path1", "path2"]} data["cap"] = {} msw_bound = data["bnd"]["ibound"].isel(layer=0, drop=False) data["cap"]["boundary"] = msw_bound diff --git a/imod/tests/test_mf6/test_mf6_simulation.py b/imod/tests/test_mf6/test_mf6_simulation.py index ef7386a01..35406a52c 100644 --- a/imod/tests/test_mf6/test_mf6_simulation.py +++ b/imod/tests/test_mf6/test_mf6_simulation.py @@ -615,6 +615,7 @@ def test_from_imod5__has_cap_data(imod5_dataset): imod5_data = deepcopy(imod5_dataset[0]) period_data = imod5_dataset[1] + imod5_data["extra"] = {"paths": ["path1", "path2"]} imod5_data["cap"] = {} msw_bound = imod5_data["bnd"]["ibound"].isel(layer=0, drop=True) imod5_data["cap"]["boundary"] = msw_bound diff --git a/imod/tests/test_mf6/test_mf6_wel.py b/imod/tests/test_mf6/test_mf6_wel.py index c37f71e7d..94d240f87 100644 --- a/imod/tests/test_mf6/test_mf6_wel.py +++ b/imod/tests/test_mf6/test_mf6_wel.py @@ -1124,14 +1124,16 @@ def test_logmessage_for_missing_filter_settings( @pytest.mark.unittest_jit -def test_from_imod5_cap_data__grid(cap_data_sprinkling_grid): +def test_from_imod5_cap_data__grid(cap_data_sprinkling_grid, cap_coupled_dis_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) + well = LayeredWell.from_imod5_cap_data( + cap_data_sprinkling_grid, cap_coupled_dis_grid + ) # Assert ds = well.dataset @@ -1143,20 +1145,26 @@ def test_from_imod5_cap_data__grid(cap_data_sprinkling_grid): @pytest.mark.unittest_jit @pytest.mark.timeout(300) -def test_from_imod5_cap_data__big_grid(cap_data_sprinkling_grid__big): +def test_from_imod5_cap_data__big_grid( + cap_data_sprinkling_grid__big, cap_coupled_dis_grid__big +): """Test if performance is acceptable for large grids.""" # Arrange bnd_2d = cap_data_sprinkling_grid__big["cap"]["boundary"] layer = xr.DataArray([1, 1], coords={"layer": [1, 2]}, dims=["layer"]) bnd = layer * bnd_2d # Act - well = LayeredWell.from_imod5_cap_data(cap_data_sprinkling_grid__big) + well = LayeredWell.from_imod5_cap_data( + cap_data_sprinkling_grid__big, cap_coupled_dis_grid__big + ) mf6_wel = well.to_mf6_pkg(bnd.astype(bool), bnd, bnd, bnd) # Assert assert mf6_wel.dataset.sizes["ncellid"] == bnd_2d.size @pytest.mark.unittest_jit -def test_from_imod5_cap_data__points(cap_data_sprinkling_points): +def test_from_imod5_cap_data__points(cap_data_sprinkling_points, cap_coupled_dis_grid): with pytest.raises(NotImplementedError): - LayeredWell.from_imod5_cap_data(cap_data_sprinkling_points) + LayeredWell.from_imod5_cap_data( + cap_data_sprinkling_points, cap_coupled_dis_grid + ) diff --git a/imod/tests/test_msw/test_model.py b/imod/tests/test_msw/test_model.py index 80f136466..aeaaf887c 100644 --- a/imod/tests/test_msw/test_model.py +++ b/imod/tests/test_msw/test_model.py @@ -1,5 +1,6 @@ from copy import copy from pathlib import Path +from typing import cast import pytest import xarray as xr @@ -15,6 +16,7 @@ from imod.typing import GridDataArray, Imod5DataDict from imod.typing.grid import zeros_like from imod.util.regrid import RegridderWeightsCache +from imod.util.spatial import empty_2d def test_msw_model_write(msw_model, coupled_mf6_model, coupled_mf6wel, tmp_path): @@ -320,7 +322,7 @@ def case_no_scaling_factors( 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["conductivity_factor"] = xr.DataArray([1.0], **layer_kwargs) return imod5_cap_data, has_scaling_factor def case_constants( @@ -332,7 +334,7 @@ def case_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["conductivity_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) @@ -344,19 +346,46 @@ def case_constants( return imod5_cap_data, has_scaling_factor +class TargetDisCases: + def case_dis_same_grid( + self, coupled_mf6_model: mf6.Modflow6Simulation + ) -> tuple[mf6.StructuredDiscretization :, xr.DataArray]: + """ + Modflow6 discretization on the same grid as the CAP data, thus no regridding needed. + """ + dis_pkg = cast(mf6.StructuredDiscretization, coupled_mf6_model["GWF_1"]["dis"]) + times = coupled_mf6_model["time_discretization"].dataset.coords["time"] + return dis_pkg, times + + def case_dis_different_grid( + self, coupled_mf6_model: mf6.Modflow6Simulation + ) -> tuple[mf6.StructuredDiscretization :, xr.DataArray]: + """ + Modflow6 discretization on different grid as the CAP data, thus regridding needed. + """ + dis_pkg = cast(mf6.StructuredDiscretization, coupled_mf6_model["GWF_1"]["dis"]) + target_grid = empty_2d( + dx=0.5, xmin=0.75, xmax=3.25, dy=-0.5, ymin=0.75, ymax=3.25 + ) + times = coupled_mf6_model["time_discretization"].dataset.coords["time"] + regridder_types = RegridderWeightsCache() + dis_pkg_regridded = dis_pkg.regrid_like(target_grid, regridder_types) + return dis_pkg_regridded, times + + @pytest.mark.unittest_jit @parametrize_with_cases("imod5_data, has_scaling_factor", cases=Imod5DataCases) +@parametrize_with_cases("dis_pkg, times", cases=TargetDisCases) def test_import_from_imod5( imod5_data: Imod5DataDict, has_scaling_factor: bool, meteo_grids: tuple[GridDataArray], - coupled_mf6_model: mf6.Modflow6Simulation, + dis_pkg: mf6.StructuredDiscretization, + times: xr.DataArray, 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 @@ -384,6 +413,8 @@ def test_import_from_imod5( assert "time" in model["time_oc"].dataset.sizes.keys() assert len(model["meteo_grid"].dataset.dims) == 0 assert ("scaling_factor" in model.keys()) == has_scaling_factor + assert dis_pkg.dataset["x"].equals(model["grid"].dataset["x"]) + assert dis_pkg.dataset["y"].equals(model["grid"].dataset["y"]) @pytest.mark.unittest_jit @@ -404,7 +435,7 @@ def test_import_from_imod5_and_write( 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) + well_pkg = mf6.LayeredWell.from_imod5_cap_data(imod5_data, dis_pkg) mf6_wel_pkg = well_pkg.to_mf6_pkg( active, dis_pkg["top"], dis_pkg["bottom"], npf_pkg["k"] ) diff --git a/imod/tests/test_msw/test_scaling_factors.py b/imod/tests/test_msw/test_scaling_factors.py index d2610577e..e61e5b1f9 100644 --- a/imod/tests/test_msw/test_scaling_factors.py +++ b/imod/tests/test_msw/test_scaling_factors.py @@ -139,7 +139,7 @@ def test_from_imod5_data(fixed_format_parser): 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"]["conductivity_factor"] = scale_rural imod5_data["cap"]["perched_water_table_level"] = depth_perched_water_table scaling_factors = ScalingFactors.from_imod5_data(imod5_data)