Source code for imod.msw.model
import collections
import inspect
import warnings
from copy import copy, deepcopy
from datetime import datetime
from pathlib import Path
from typing import Any, Optional, Union, cast
import cftime
import jinja2
import numpy as np
import tomli
import xarray as xr
import imod.msw
from imod.common.constants import MaskValues
from imod.common.interfaces.idict import IDict
from imod.common.serializer import EngineType
from imod.common.utilities.clip import clip_by_grid
from imod.common.utilities.dump_model import dump_model
from imod.common.utilities.partitioninfo import create_partition_info
from imod.common.utilities.regrid import regrid_imod5_cap_data
from imod.common.utilities.version import prepend_content_with_version_info
from imod.mf6.dis import StructuredDiscretization
from imod.mf6.mf6_wel_adapter import Mf6Wel
from imod.msw import GridData
from imod.msw.copy_files import FileCopier
from imod.msw.coupler_mapping import CouplerMapping
from imod.msw.idf_mapping import IdfMapping
from imod.msw.infiltration import Infiltration
from imod.msw.initial_conditions import (
InitialConditionsEquilibrium,
InitialConditionsPercolation,
InitialConditionsRootzonePressureHead,
InitialConditionsSavedState,
)
from imod.msw.landuse import LanduseOptions
from imod.msw.meteo_grid import MeteoGrid, MeteoGridCopy
from imod.msw.meteo_mapping import (
EvapotranspirationMapping,
MeteoMapping,
PrecipitationMapping,
)
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.mask import (
MetaSwapActive,
mask_and_broadcast_cap_data,
mask_and_broadcast_pkg_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.regrid import RegridderWeightsCache
from imod.util.time import to_datetime_internal
REQUIRED_PACKAGES = (
GridData,
CouplerMapping,
Infiltration,
LanduseOptions,
EvapotranspirationMapping,
PrecipitationMapping,
IdfMapping,
TimeOutputControl,
AnnualCropFactors,
)
METEO_PACKAGES = (
MeteoGrid,
MeteoGridCopy,
)
INITIAL_CONDITIONS_PACKAGES = (
InitialConditionsEquilibrium,
InitialConditionsPercolation,
InitialConditionsRootzonePressureHead,
InitialConditionsSavedState,
)
DEFAULT_SETTINGS: dict[str, Any] = {
"vegetation_mdl": 1,
"evapotranspiration_mdl": 1,
"saltstress_mdl": 0,
"surfacewater_mdl": 0,
"infilimsat_opt": 0,
"netcdf_per": 0,
"postmsw_opt": 0,
"dtgw": 1.0,
"dtsw": 1.0,
"ipstep": 2,
"nxlvage_dim": 366,
"co2": 404.32,
"fact_beta2": 1.0,
"rcsoil": 0.15,
"iterur1": 3,
"iterur2": 5,
"tdbgsm": 91.0,
"tdedsm": 270.0,
"clocktime": 0,
}
class Model(collections.UserDict[str, Any]):
def __setitem__(self, key, value):
# TODO: Add packagecheck
super().__setitem__(key, value)
def update(self, *args, **kwargs):
for k, v in dict(*args, **kwargs).items():
self[k] = v
[docs]
class MetaSwapModel(Model, IDict):
"""
Contains data and writes consistent model input files
Parameters
----------
unsaturated_database: Path-like or str
Path to the MetaSWAP soil physical database folder.
settings: dict
"""
_pkg_id = "model"
_file_name = "para_sim.inp"
_model_name = "MSW"
_template = jinja2.Template(
"{%for setting, value in settings.items()%}{{setting}} = {{value}}\n{%endfor%}"
)
[docs]
def __init__(
self,
unsaturated_database: Path | str,
settings: Optional[dict[str, Any]] = None,
starttime: Optional[str] = None,
):
super().__init__()
if settings is None:
self.simulation_settings = copy(DEFAULT_SETTINGS)
else:
self.simulation_settings = settings
self.starttime = starttime
self.simulation_settings["unsa_svat_path"] = unsaturated_database
def _render_unsaturated_database_path(self, unsaturated_database: Union[str, Path]):
# Force to Path object
unsaturated_database = Path(unsaturated_database)
# Render to string for MetaSWAP
if unsaturated_database.is_absolute():
return f'"{unsaturated_database}\\"'
else:
# TODO: Test if this is how MetaSWAP accepts relative paths
return f'"${unsaturated_database}\\"'
def _check_required_packages(self) -> None:
pkg_types_included = {type(pkg) for pkg in self.values()}
required_packages_set = cast(set[type[Any]], set(REQUIRED_PACKAGES))
missing_packages = required_packages_set - pkg_types_included
if len(missing_packages) > 0:
raise ValueError(
f"Missing the following required packages: {missing_packages}"
)
meteo_set = pkg_types_included & set(METEO_PACKAGES)
if len(meteo_set) < 1:
raise ValueError(f"Missing meteo package, assign one of {METEO_PACKAGES}")
initial_condition_set = pkg_types_included & set(INITIAL_CONDITIONS_PACKAGES)
if len(initial_condition_set) < 1:
raise ValueError(
"Missing InitialCondition package, assign one of "
f"{INITIAL_CONDITIONS_PACKAGES}"
)
elif len(initial_condition_set) > 1:
raise ValueError(
"Multiple InitialConditions assigned, choose one of "
f"{initial_condition_set}"
)
def _check_landuse_indices_in_lookup_options(self):
grid_key = self.get_pkgkey(GridData)
landuse_options_key = self.get_pkgkey(LanduseOptions)
indices_in_grid = set(self[grid_key]["landuse"].values.ravel())
indices_in_options = set(
self[landuse_options_key].dataset.coords["landuse_index"].values
)
missing_indices = indices_in_grid - indices_in_options
if len(missing_indices) > 0:
raise ValueError(
"Found the following landuse indices in GridData which "
f"were not in LanduseOptions: {missing_indices}"
)
def _check_vegetation_indices_in_annual_crop_factors(self):
landuse_options_key = self.get_pkgkey(LanduseOptions)
annual_crop_factors_key = self.get_pkgkey(AnnualCropFactors)
indices_in_options = set(
np.unique(self[landuse_options_key]["vegetation_index"])
)
indices_in_crop_factors = set(
self[annual_crop_factors_key].dataset.coords["vegetation_index"].values
)
missing_indices = indices_in_options - indices_in_crop_factors
if len(missing_indices) > 0:
raise ValueError(
"Found the following vegetation indices in LanduseOptions "
f"which were not in AnnualCropGrowth: {missing_indices}"
)
def _has_file_copier(self) -> bool:
pkg_types_included = {type(pkg) for pkg in self.values()}
return FileCopier in pkg_types_included
def _get_starttime(self):
if self.starttime is not None:
starttime = to_datetime_internal(self.starttime, use_cftime=False)
else:
starttime = self._get_minimum_package_time()
year, time_since_start_year = to_metaswap_timeformat([starttime])
year = int(year.item())
time_since_start_year = float(time_since_start_year.item())
return year, time_since_start_year
def _get_minimum_package_time(self):
"""
Loop over all packages to get the minimum time.
MetaSWAP requires a starttime in its simulation settings (para_sim.inp)
"""
starttimes = []
for pkgname in self:
ds = self[pkgname].dataset
if "time" in ds.coords:
starttimes.append(ds["time"].min().values)
if len(starttimes) == 0:
raise ValueError(
"No package with a coordinate 'time' found. Please add a MeteoGrid or TimeOutputControl package."
)
starttime = min(starttimes)
return starttime
@property
def nsubunits(self) -> int:
"""
Get the number of subunits in the model from the first dataset having a subunit dimension.
Defaults to 1 if no package has a subunit dimension.
Returns
-------
int
The number of subunits in the model.
"""
nsub = 1
for pkg in self.values():
if "subunit" in pkg.dataset.dims:
nsub = pkg.dataset.dims["subunit"]
break
return nsub
[docs]
def get_pkgkey(
self, pkg_type: type[MetaSwapPackage], optional_package: bool = False
) -> str | None:
"""
Get the package key for a package of type ``pkg_type``. Returns the
first occurrence of the package type.
Parameters
----------
pkg_type: type[MetaSwapPackage]
Type of the package to get the key for.
optional_package: bool
If True, the method will not raise an error if the package is not
found. Method returns None in this case.
Returns
-------
str
The key of the package of type ``pkg_type``.
"""
for pkg_key, pkg in self.items():
if isinstance(pkg, pkg_type):
return pkg_key
if not optional_package:
raise KeyError(f"Could not find package of type: {pkg_type}")
return None
def _get_pkg_key(
self, pkg_type: type[MetaSwapPackage], optional_package: bool = False
) -> str | None:
""" "
Preserves backwards compatibility with old code (primod) that used this.
"""
warnings.warn(
"Method '_get_pkg_key' is deprecated, use 'get_pkgkey' instead.",
DeprecationWarning,
)
return self.get_pkgkey(pkg_type, optional_package)
def _model_checks(self, validate: Optional[bool] = True):
if validate and not self._has_file_copier():
self._check_required_packages()
self._check_vegetation_indices_in_annual_crop_factors()
self._check_landuse_indices_in_lookup_options()
def _write_simulation_settings(self, directory: Path) -> None:
"""
Write simulation settings to para_sim.inp.
Parameters
----------
directory: Path or str
directory to write model in.
"""
simulation_settings = deepcopy(self.simulation_settings)
# Add time settings
year, time_since_start_year = self._get_starttime()
simulation_settings["iybg"] = year
simulation_settings["tdbg"] = time_since_start_year
# Add IdfMapping settings
idf_key = self.get_pkgkey(IdfMapping)
idf_pkg = cast(IdfMapping, self[idf_key]) # type: ignore[index]
simulation_settings.update(idf_pkg._get_output_settings())
simulation_settings["unsa_svat_path"] = self._render_unsaturated_database_path(
simulation_settings["unsa_svat_path"]
)
filename = directory / self._file_name
rendered = self._template.render(settings=simulation_settings)
# Prepend version information
rendered = prepend_content_with_version_info(rendered, comment_char="*")
with open(filename, "w") as f:
f.write(rendered)
[docs]
def write(
self,
directory: Union[str, Path],
mf6_dis: StructuredDiscretization,
mf6_wel: Optional[Mf6Wel] = None,
validate: Optional[bool] = True,
):
"""
Write packages and simulation settings (para_sim.inp).
Parameters
----------
directory: Path or str
directory to write model in.
"""
# Model checks
self._model_checks(validate)
# Force to Path
directory = Path(directory)
directory.mkdir(exist_ok=True, parents=True)
# Write simulation settings
self._write_simulation_settings(directory)
# Get index and svat
grid_key = self.get_pkgkey(GridData)
grid_pkg = cast(GridData, self[grid_key]) # type: ignore[index]
index, svat = grid_pkg.generate_isactive_svat_arrays()
# write package contents
for pkgname in self:
self[pkgname].write(directory, index, svat, mf6_dis, mf6_wel)
@classmethod
def from_file(cls, directory, modelname):
pkg_classes = {
name: pkg_cls
for name, pkg_cls in inspect.getmembers(imod.msw, inspect.isclass)
if issubclass(pkg_cls, MetaSwapPackage)
}
modeldirectory = Path(directory) / modelname
toml_path = modeldirectory / f"{modelname}.toml"
with open(toml_path, "rb") as f:
toml_content = tomli.load(f)
parentdir = toml_path.parent
simulation_settings = toml_content.get("simulation_settings", {})
unsa_svat_path = simulation_settings.get("unsa_svat_path", "")
instance = cls(unsa_svat_path, simulation_settings)
for key, entry in toml_content.items():
if key != "simulation_settings":
for pkgname, path in entry.items():
pkg_cls = pkg_classes[key]
instance[pkgname] = pkg_cls.from_file(parentdir / path)
return instance
[docs]
def dump(
self,
directory: Union[str, Path],
modelname: Optional[str] = None,
validate: bool = True,
mdal_compliant: bool = False,
crs: Optional[str] = None,
engine: EngineType = "netcdf4",
):
"""
Dump model packages to netCDF files and create a toml file with paths to these files.
The MetaSWAP model can be reloaded from the dumped files using the :func:`from_file` method.
Parameters
----------
directory: Path or str
Directory to dump model in. A subdirectory with the name of the model will be created in this directory, and the files will be dumped there.
modelname: str, optional
Name of the model. This will be used as the name of the subdirectory where the files are dumped, and in the name of the toml file. If not provided, it defaults to the value of ``self._model_name``.
validate: bool, optional
Whether to perform validation before dumping the model. If True, the model will be validated using the validation functionality in iMOD. If validation errors are found, a ValidationError is raised and the model is not dumped. Default is True.
mdal_compliant: bool, optional
Whether to write the files in a format compliant with the MDAL specification. This can be used to make the files compatible with software that supports MDAL, such as QGIS. Default is False.
crs: str, optional
Coordinate reference system to use in the dumped files. This should be a string in a format recognized by the pyproj library, for example "EPSG:28992". If not provided, no CRS information is included in the files.
engine: EngineType, optional
File engine used to write packages.
engine : str, optional
File engine used to write packages. Options are ``'netcdf4'``,
``'zarr'``, and ``'zarr.zip'``. NetCDF4 is readable by many other
softwares, for example QGIS. Zarr is optimized for big data, cloud
storage and parallel access. The ``'zarr.zip'`` option is an
experimental option which creates a zipped zarr store in a single
file, which is easier to copy and automatically compresses data as
well. Default is ``'netcdf4'``.
Returns
-------
Path
Path to the created toml file which contains the paths to the dumped package files. The package files are dumped in the same directory as the toml file.
Example
-------
>>> tmp_path = tmpdir_factory.mktemp(name)
>>> toml_path = msw_model.dump(tmp_path, name, engine=engine, validate=False)
>>> back = MetaSwapModel.from_file(tmp_path, name)
"""
toml_path = dump_model(
self,
directory=directory,
modelname=modelname,
validate=validate,
mdal_compliant=mdal_compliant,
crs=crs,
engine=engine,
)
return toml_path
[docs]
def regrid_like(
self,
mf6_regridded_dis: StructuredDiscretization,
regrid_cache: Optional[RegridderWeightsCache] = None,
) -> "MetaSwapModel":
"""
Creates a model by regridding the packages of this model to another
discretization. It regrids all the arrays in the package using the
default regridding methods. At the moment only regridding to a different
planar grid is supported, meaning ``target_grid`` has different ``"x"``
and ``"y"`` or different ``cell2d`` coords.
Parameters
----------
mf6_regridded_dis: StructuredDiscretization
Modflow6 Discretization with same discretization as the one we want
to regrid the package to.
regrid_cache: RegridderWeightsCache, optional
stores regridder weights for different regridders. Can be used to
speed up regridding, if the same regridders are used several times
for regridding different arrays.
Returns
-------
A model with similar packages to the input model, and with all the
data-arrays regridded to another discretization, similar to the one used
in input argument "mf6_regridded_dis"
"""
unsat_database = cast(str, self.simulation_settings["unsa_svat_path"])
regridded_model = MetaSwapModel(unsat_database)
target_grid = mf6_regridded_dis["idomain"]
mod2svat_name = None
for pkgname in self:
msw_package = self[pkgname]
if isinstance(msw_package, CouplerMapping):
# there can be only one couplermapping
mod2svat_name = pkgname
elif msw_package._is_regridding_supported():
regridded_package = msw_package.regrid_like(
target_grid, regrid_cache, None
)
else:
raise ValueError(f"package {pkgname} cannot be regridded")
regridded_model[pkgname] = regridded_package
if mod2svat_name is not None:
regridded_model[mod2svat_name] = CouplerMapping()
return regridded_model
[docs]
def mask_all_packages(
self,
mask: GridDataArray,
):
"""
This function applies a mask to all packages in a model. The mask must
be presented as a GridDataArray, which contains idomain-like integers.
The mask is applied to all packages in the model, and the values in the mask determine which cells are active and which are inactive. The mask is applied to all packages, regardless of whether they have a subunit dimension or not.
Parameters
----------
mask: GridDataArray
idomain-like integers. >0 sets cells to active, 0 sets cells to inactive,
mask is applied on a per-subunit basis if the mask grid has a subunit dimension.
If the package does not have a subunit dimension, the combined mask grid over all of its subunits is applied,
i.e. a cell is set to active if it is active in any of the subunits.
Example
-------
>>> mask = 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 = {
>>> "x" : [1.0, 2.0, 3.0],
>>> "y" : [3.0, 2.0, 1.0],
>>> "dx" : 1.0,
>>> "dy" : 1.0,
>>> "subunit" : [0, 1]
>>> }
>>> )
"""
if "subunit" in mask.dims:
mask_all = mask.any(dim="subunit")
msw_active = MetaSwapActive(all=mask_all, per_subunit=mask)
else:
mask_per_subunit = mask.expand_dims(dim={"subunit": self.nsubunits})
msw_active = MetaSwapActive(all=mask, per_subunit=mask_per_subunit)
for pkg in self.values():
if "x" in pkg.dataset.dims and "y" in pkg.dataset.dims:
data_dict = {
key: pkg.dataset[key] for key in pkg.dataset.data_vars.keys()
}
masked_data = mask_and_broadcast_pkg_data(pkg, data_dict, msw_active)
for key, data in masked_data.items():
pkg.dataset[key] = data
return
[docs]
def clip_box(
self,
time_min: Optional[cftime.datetime | np.datetime64 | str] = None,
time_max: Optional[cftime.datetime | np.datetime64 | str] = None,
x_min: Optional[float] = None,
x_max: Optional[float] = None,
y_min: Optional[float] = None,
y_max: Optional[float] = None,
) -> "MetaSwapModel":
"""
Clip a model by a bounding box (time, y, x). If a package of type
:class:`imod.msw.MeteoGridCopy` is present, packages of type
:class:`imod.msw.PrecipitationMapping` and
:class:`imod.msw.EvapotranspirationMapping` will not be clipped.
Otherwise incorrect mappings to meteo grids referenced in
``mete_grid.inp`` copied by :class:`imod.msw.MeteoGridCopy` would be
computed.
Parameters
----------
time_min: optional, np.datetime64
Start time to select. Data will be forward filled to this date. If
time_min is before the start time of the dataset, data is
backfilled.
time_max: optional
End time to select.
x_min: optional, float
Minimum x-coordinate to select.
x_max: optional, float
Maximum x-coordinate to select.
y_min: optional, float
Minimum y-coordinate to select.
y_max: optional, float
Maximum y-coordinate to select.
Returns
-------
clipped : MetaSwapModel
A new model that is clipped to the specified bounding box.
Examples
--------
Slicing intervals may be half-bounded, by providing None:
To select 500.0 <= x <= 1000.0:
>>> msw_model.clip_box(x_min=500.0, x_max=1000.0)
To select x <= 1000.0:
>>> msw_model.clip_box(x_max=1000.0)``
To select x >= 500.0:
>>> msw_model.clip_box(x_min=500.0)
To select a time interval, you can use datetime64:
>>> msw_model.clip_box(time_min=np.datetime64("2020-01-01"), time_max=np.datetime64("2020-12-31"))
"""
settings = deepcopy(self.simulation_settings)
unsa_svat_path = settings.pop("unsa_svat_path")
has_meteogrid_copy = MeteoGridCopy in [type(pkg) for pkg in self.values()]
clipped = type(self)(unsa_svat_path, settings)
for key, pkg in self.items():
# Skip clipping MeteoMapping if MeteoGridCopy is present, as meteo
# grid is independent of model grid and we do not want to perform
# transformations on meteo data in this case.
if has_meteogrid_copy and isinstance(pkg, MeteoMapping):
clipped[key] = deepcopy(pkg)
else:
clipped[key] = pkg.clip_box(
time_min=time_min,
time_max=time_max,
x_min=x_min,
x_max=x_max,
y_min=y_min,
y_max=y_max,
)
return clipped
def split(
self,
submodel_labels: GridDataArray,
) -> dict[str, "MetaSwapModel"]:
"""
Split a MetaSWAP model in different partitions using a submodel_labels
array. Note: for specifying meteorological grid data, splitting is only
supported when the MeteoGridCopy instance is being used.
Parameters
----------
submodel_labels: xr.DataArray or xu.UgridDataArray
A grid that defines how the simulation will be split. The array
should have the same topology as the domain being split, i.e.
similar shape as a layer in the domain. The values in the array
indicate to which partition a cell belongs. The values should be
zero or greater.
Returns
-------
dict[str, MetaSwapModel]
A mapping from generated submodel names to the corresponding clipped
MetaSWAP submodel.
Examples
--------
>>> submodel_labels = mf6_sim.create_partition_labels(n_partitions=4)
>>> msw_splitted = msw.split(submodel_labels)
"""
has_meteogrid = MeteoGrid in [type(pkg) for pkg in self.values()]
if has_meteogrid:
raise ValueError(
"Splitting packages of type MeteoGrid is not supported, use MeteoGridCopy instead."
)
settings = deepcopy(self.simulation_settings)
starttime = self.starttime
unsa_svat_path = settings.pop("unsa_svat_path")
partition_info_list = create_partition_info(submodel_labels)
# Initialize mapping from partition IDs to models
partition_id_to_models: dict[int, dict[str, MetaSwapModel]] = {}
for submodel_partition_info in partition_info_list:
partition_id_to_models[submodel_partition_info.id] = {}
partitioned_submodels = {}
submodel_to_partition = {}
# Create empty MetaSWAP model for each partition
for submodel_partition_info in partition_info_list:
submodel_name = f"{self._model_name}_{submodel_partition_info.id}"
submodel = MetaSwapModel(unsa_svat_path, settings, starttime)
partitioned_submodels[submodel_name] = submodel
submodel_to_partition[submodel_name] = submodel_partition_info
partition_id_to_models[submodel_partition_info.id][submodel_name] = submodel
# First, handle the grid data and determine the overlap.
grid_key = self.get_pkgkey(GridData)
grid_pkg = cast(GridData, self[grid_key]) # type: ignore[index]
is_in_active_domain = {}
for submodel_name, submodel in partitioned_submodels.items():
partition_info = submodel_to_partition[submodel_name]
sliced_grid_pkg = cast(
GridData, clip_by_grid(grid_pkg, partition_info.active_domain)
)
sliced_isactive = sliced_grid_pkg._generate_isactive_array().values
# Add package to model if it has data in the active domain.
if bool(sliced_isactive.any()):
is_in_active_domain[submodel_name] = True
submodel[grid_key] = sliced_grid_pkg
else:
is_in_active_domain[submodel_name] = False
# Second, add other packages to each partitioned submodel.
for submodel_name, submodel in partitioned_submodels.items():
partition_info = submodel_to_partition[submodel_name]
if not is_in_active_domain[submodel_name]:
continue
# Add packages to models
for pkg_name, pkg in self.items():
if isinstance(pkg, GridData):
continue
if isinstance(pkg, MeteoMapping):
submodel[pkg_name] = deepcopy(pkg)
else:
# Slice and add the package to the partitioned model
sliced_package = clip_by_grid(pkg, partition_info.active_domain)
submodel[pkg_name] = sliced_package
return partitioned_submodels
[docs]
@classmethod
def from_imod5_data(
cls,
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.
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, 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"])
# 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_regridded, target_dis
)
cap_data_masked = mask_and_broadcast_cap_data(
imod5_regridded["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 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)
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
# make a read function for packages from netcdf