API Documentation

nlmod

read

nlmod.read.ahn._assert_as_data_array_is_none(as_data_array: bool | None) None[source]
nlmod.read.ahn._delete_tiles_from_config_in_data(config: dict = None)[source]

Delete the location of the map tiles of AHN from the data-directory of nlmod.

This method deletes the files downloaded with _save_tiles_from_config_in_data from the data directory again.

Parameters:

config (dict, optional) – A dictionary with properties of the data sources of the different AHN-versions. When None, the configuration is retreived from the method get_configuration(). The default is None.

Return type:

None.

nlmod.read.ahn._download_ahn_ellipsis(extent: list[float], identifier: str, merge_tiles: bool = True, cut_extent: bool = True, cachedir=None, cachename=None, **kwargs) DataArray | list[DataArray][source]

Download and merge AHN tiles from the ellipsis.

Parameters:
  • extent (list, tuple or np.array) – extent

  • identifier (str) – The identifier determines the AHN-version, the resolution and the type of height data.

  • merge_tiles (bool, optional) – If True, the function returns a merged DataArray. If False, the function returns a list of DataArrays with the original tiles. The default is True.

  • cut_extent (bool, optional) – If True, only keep the requested extent from the data. The defualts is True.

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

DataArray of the AHN

Return type:

xr.DataArray

nlmod.read.ahn._download_ahn_hwh(extent: list[float], version: str = 'AHN4', data_kind: str = 'DTM', tile_size: str = '5x6.25km', resolution: float = 5.0, merge_tiles: bool = True, cut_extent: bool = True, config: dict = None, cachedir=None, cachename=None)[source]

Download AHN-data.

Parameters:
  • extent (list, tuple or np.array) – The extent to be downloaded, consisting of 4 floats: xmin, xmax, ymin, ymax.

  • version (str, optional) – The AHN_version, which can be ‘AHN2’, ‘AHN3’, ‘AHN4’, ‘AHN5’ and ‘AHN6’. The default is “AHN4”, as this is available in the whole of the Netherlands.

  • data_kind (str, optional) – The kind of data. This can be ‘DTM’ (terrain elevation) or “DSM” (surface elevation). The default is “DTM”.

  • tile_size (str, optional) – The size of the map tiles that the data is offered by the webservice. This can be ‘5x6.25km’ or ‘1x1km’. The default is ‘5x6.25km’.

  • resolution (float, optional) – The resolution of the AHN-data, which can be 0.5 and 5.0 m. The default is 5.0.

  • merge_tiles (bool, optional) – If True, the function returns a merged DataArray. If False, the function returns a list of DataArrays with the original tiles. The default is True.

  • cut_extent (bool, optional) – If True, only keep the requested extent from the data. The defualts is True.

  • config (dict, optional) – A dictionary with properties of the data sources of the different AHN-versions. When None, the configuration is retreived from the method get_configuration(). The default is None.

Raises:
  • AssertionError – If the specified value for resolution or data_kind is not supported

  • KeyError – When there is no data-source that matches the requested parameters.

:raises cachedir : str or None, optional: directory to save cache. If None no cache is used. Default is None. :raises cachename : str or None, optional: filename of netcdf cache. If None no cache is used. Default is None.

Returns:

A DataArray with the AHN-data, or a list of DataArrays with AHN-data if merge_tiles=False.

Return type:

xr.DataArray or list of xr.DataArray

nlmod.read.ahn._download_and_save_json_file(url, fname, timeout=120.0)[source]
nlmod.read.ahn._download_tiles_ellipsis(extent: list[float] | None = None, crs: int = 28992, timeout: float = 120.0, base_url: str = 'https://api.ellipsis-drive.com/v3', path_id: str = 'a9d410ad-a2f6-404c-948a-fdf6b43e77a6', timestamp_id: str = '05931403-2510-43af-9cc3-f60a066d4482') GeoDataFrame[source]
nlmod.read.ahn._download_tiles_hwh(extent: list[float], version: str = 'AHN4', data_kind: str = 'DTM', tile_size: str = '5x6.25km', resolution: float = 5.0, config: dict = None, timeout: float = 120.0, cachedir=None, cachename=None)[source]

Download tiles that contain AHN-data.

Parameters:
  • extent (list, tuple or np.array) – The extent to be downloaded, consisting of 4 floats: xmin, xmax, ymin, ymax.

  • version (str, optional) – The AHN_version, which can be ‘AHN2’, ‘AHN3’, ‘AHN4’, ‘AHN5’ and ‘AHN6’. The default is “AHN4”, as this is available in the whole of the Netherlands.

  • data_kind (str, optional) – The kind of data. This can be ‘DTM’ (terrain elevation) or “DSM” (surface elevation). The default is “DTM”.

  • tile_size (str, optional) – The size of the map tiles that the data is offered by the webservice. This can be ‘5x6.25km’ or ‘1x1km’. The default is ‘5x6.25km’.

  • resolution (float, optional) – The resolution of the AHN-data, which can be 0.5 and 5.0 m. The default is 5.0.

  • config (dict, optional) – A dictionary with properties of the data sources of the different AHN-versions. When None, the configuration is retreived from the method get_configuration(). The default is None.

  • timeout (float, optional) – The amount of seconds to wait for the server to send data before giving up. The default is 120.

Raises:
  • AssertionError – If the specified value for resolution or data_kind is not supported

  • KeyError – When there is no data-source that matches the requested parameters.

:raises cachedir : str or None, optional: directory to save cache. If None no cache is used. Default is None. :raises cachename : str or None, optional: filename of netcdf cache. If None no cache is used. Default is None.

Returns:

A GeoDataFrame containing the location of the map tiles, and url’s with the location of the AHN-data.

Return type:

gpd.GeoDataFrame

nlmod.read.ahn._get_tiles_from_file(fname: str, extent: list[float] | None = None, crs: int = 28992) GeoDataFrame[source]
nlmod.read.ahn._rename_identifier(identifier: str) str[source]
nlmod.read.ahn._round_coordinates(geom: Geometry, ndigits: int = 2) Geometry[source]
nlmod.read.ahn._save_tiles_from_config_in_data(config: dict = None)[source]

Saves the location of the map tiles of AHN to the data-directory of nlmod.

For each available combination of version, data_kind, tile_size and resolution, a json-file is downloaded and saved to the folder ahn in the data-directory of nlmod. When requesting ahn-data after running this method, the location of the map-tiles is then not downloaded anymore, but taken from the previously downloaded json-files.

Parameters:

config (dict, optional) – A dictionary with properties of the data sources of the different AHN-versions. When None, the configuration is retreived from the method get_configuration(). The default is None.

Return type:

None.

nlmod.read.ahn._update_ellipsis_tiles_in_data() None[source]
nlmod.read.ahn.discretize_ahn(ds: Dataset, ahn_da: DataArray, method: str = 'average', cachedir=None, cachename=None) Dataset[source]

Discretize ahn data to model the model grid.

Parameters:
  • ds (xr.Dataset) – dataset with the model information.

  • ahn_da (xr.DataArray) – ahn data within model extent.

  • method (str, optional) – Method used to resample ahn to grid of ds. See documentation of nlmod.resample.structured_da_to_ds for possible values. The default is ‘average’.

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

ds_out – Dataset with the ahn variable.

Return type:

xr.Dataset

nlmod.read.ahn.download_ahn(extent: list[float], identifier: str = None, merge_tiles=True, cachedir=None, cachename=None, **kwargs) DataArray[source]

Download ahn data within an extent.

Parameters:
  • extent (list, tuple or np.array) – The extent to be downloaded, consisting of 4 floats: xmin, xmax, ymin, ymax.

  • version (str, optional) – The AHN_version, which can be ‘AHN2’, ‘AHN3’, ‘AHN4’, ‘AHN5’ and ‘AHN6’. version is ignored if identifier is specified. The default is “AHN4”, as this is available in the whole of the Netherlands.

  • data_kind (str, optional) – The kind of data. This can be ‘DTM’ (terrain elevation) or “DSM” (surface elevation). data_kind is ignored if identifier is specified. The default is “DTM”.

  • tile_size (str, optional) – The size of the map tiles that the data is offered by the webservice. This can be ‘5x6.25km’ or ‘1x1km’. tile_size is ignored if identifier is specified. The default is ‘5x6.25km’.

  • resolution (float, optional) – The resolution of the AHN-data, which can be 0.5 and 5.0 m. resolution is ignored if identifier is specified. The default is 5.0.

  • merge_tiles (bool, optional) – If True, the function returns a merged DataArray. If False, the function returns a list of DataArrays with the original tiles. The default is True.

  • cut_extent (bool, optional) – If True, only keep the requested extent from the data. The defualts is True.

  • config (dict, optional) – A dictionary with properties of the data sources of the different AHN-versions. When None, the configuration is retreived from the method get_configuration(). The default is None.

  • identifier (str, optional) –

    The identifier determines the AHN-version, the resolution and the type of height data. Possible values are (casing is important):

    AHN1: ‘AHN1 maaiveldmodel (DTM) 5m’ AHN2: ‘AHN2 maaiveldmodel (DTM) ½m, geïnterpoleerd’,

    ’AHN2 maaiveldmodel (DTM) ½m’, ‘AHN2 DSM ½m’, ‘AHN2 maaiveldmodel (DTM) 5m’

    AHN3: ‘AHN3 maaiveldmodel (DTM) ½m’,

    ’AHN3 DSM ½m’, ‘AHN3 maaiveldmodel (DTM) 5m’, ‘AHN3 DSM 5m’

    AHN4: ‘AHN4 maaiveldmodel (DTM) ½m’,

    ’AHN4 DSM ½m’, ‘AHN4 maaiveldmodel (DTM) 5m’, ‘AHN4 DSM 5m’

    AHN5: ‘AHN5 maaiveldmodel (DTM) 5m’,

    ’AHN5 DSM 5m’, ‘AHN5 maaiveldmodel (DTM) ½m’, ‘AHN5 DSM ½m’

    When no identifier is specified (the default), the url to download data from is taken from config, using the arguments version, data_kind, tile_size and resolution of this method. The default is None.

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

ahn_da – DataArray with the ahn variable.

Return type:

xr.DataArray

nlmod.read.ahn.download_ahn1(extent: list[float], identifier: Literal['AHN1 maaiveldmodel (DTM) 5m'] = 'AHN1 maaiveldmodel (DTM) 5m', as_data_array: bool | None = None, cachedir=None, cachename=None, **kwargs) DataArray[source]

Download AHN1.

Parameters:
  • extent (list, tuple or np.array) – extent

  • identifier (str, optional) – The identifier determines the resolution and the type of height data. The only allowed value is ‘AHN1 maaiveldmodel (DTM) 5m’. The default is “AHN1 maaiveldmodel (DTM) 5m”.

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

DataArray of the AHN

Return type:

xr.DataArray

nlmod.read.ahn.download_ahn2(extent: list[float], identifier: Literal['AHN2 maaiveldmodel (DTM) ½m, geïnterpoleerd', 'AHN2 maaiveldmodel (DTM) ½m', 'AHN2 DSM ½m', 'AHN2 maaiveldmodel (DTM) 5m'] = None, as_data_array: bool | None = None, cachedir=None, cachename=None, **kwargs) DataArray[source]

Download AHN2.

Parameters:
  • extent (list, tuple or np.array) – extent

  • identifier (str, optional) – The identifier determines the resolution and the type of height data. Possible values are ‘AHN2 maaiveldmodel (DTM) ½m, geïnterpoleerd’, ‘AHN2 maaiveldmodel (DTM) ½m’, ‘AHN2 DSM ½m’, and ‘AHN2 maaiveldmodel (DTM) 5m’. The default is None. When no identifier is specified (the default), the url to download data from is taken from config, using the keyword-arguments data_kind, tile_size and resolution, which all can be specified as kwargs.

  • **kwargs (dict) – See docstring of download_ahn for a descption of extra arguments.

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

DataArray of the AHN

Return type:

xr.DataArray

nlmod.read.ahn.download_ahn3(extent: list[float], identifier: Literal['AHN3 maaiveldmodel (DTM) ½m', 'AHN3 DSM ½m', 'AHN3 maaiveldmodel (DTM) 5m', 'AHN3 DSM 5m'] = None, as_data_array: bool | None = None, cachedir=None, cachename=None, **kwargs) DataArray[source]

Download AHN3.

Parameters:
  • extent (list, tuple or np.array) – extent

  • identifier (str, optional) – The identifier determines the resolution and the type of height data. Possible values are ‘AHN3 maaiveldmodel (DTM) ½m’, ‘AHN3 DSM ½m’, ‘AHN3 maaiveldmodel (DTM) 5m’, and ‘AHN3 DSM 5m’. The default is None. When no identifier is specified (the default), the url to download data from is taken from config, using the keyword-arguments data_kind, tile_size and resolution, which all can be specified as kwargs.

  • **kwargs (dict) – See docstring of download_ahn for a descption of extra arguments.

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

DataArray of the AHN

Return type:

xr.DataArray

nlmod.read.ahn.download_ahn4(extent: list[float], identifier: Literal['AHN4 maaiveldmodel (DTM) ½m', 'AHN4 DSM ½m', 'AHN4 maaiveldmodel (DTM) 5m', 'AHN4 DSM 5m'] = None, as_data_array: bool | None = None, cachedir=None, cachename=None, **kwargs) DataArray[source]

Download AHN4.

Parameters:
  • extent (list, tuple or np.array) – extent

  • identifier (str, optional) – The identifier determines the resolution and the type of height data. Possible values are ‘AHN4 maaiveldmodel (DTM) ½m’, ‘AHN4 DSM ½m’, ‘AHN4 maaiveldmodel (DTM) 5m’, and ‘AHN4 DSM 5m’. The default is None. When no identifier is specified (the default), the url to download data from is taken from config, using the keyword-arguments data_kind, tile_size and resolution, which all can be specified as kwargs.

  • **kwargs (dict) – See docstring of download_ahn for a descption of extra arguments.

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

DataArray of the AHN

Return type:

xr.DataArray

nlmod.read.ahn.download_ahn5(extent: list[float], identifier: Literal['AHN5 maaiveldmodel (DTM) 5m', 'AHN5 DSM 5m', 'AHN5 maaiveldmodel (DTM) ½m', 'AHN5 DSM ½m'] = None, cachedir=None, cachename=None, **kwargs) DataArray[source]

Download AHN5.

Parameters:
  • extent (list, tuple or np.array) – extent

  • identifier (str, optional) – The identifier determines the resolution and the type of height data. Possible values are ‘AHN5 maaiveldmodel (DTM) 5m’, ‘AHN5 DSM 5m’, ‘AHN5 maaiveldmodel (DTM) ½m’, and ‘AHN5 DSM ½m’. The default is None. When no identifier is specified (the default), the url to download data from is taken from config, using the keyword-arguments data_kind, tile_size and resolution, which all can be specified as kwargs.

  • **kwargs (dict) – See docstring of download_ahn for a descption of extra arguments.

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

DataArray of the AHN

Return type:

xr.DataArray

nlmod.read.ahn.download_ahn6(extent: list[float], tile_size: str = '1x1km', **kwargs)[source]

Download AHN6.

Parameters:
  • extent (list, tuple or np.array) – extent

  • tile_size (str, optional) – The size of the map tiles that the data is offered by the webservice. For AHN6, this can only be ‘1x1km’. The default is “1x1km”.

  • **kwargs (dict) – See docstring of download_ahn for a descption of extra arguments.

Returns:

DataArray of the AHN

Return type:

xr.DataArray

nlmod.read.ahn.download_ahn_along_line(line: LineString, ahn: DataArray | None = None, dx: float | None = None, num: int | None = None, method: str = 'linear', plot: bool = False) DataArray[source]

Get the height of the surface level along a line.

Parameters:
  • line (shapely.LineString) – The line along which the surface level is calculated.

  • ahn (xr.DataArray, optional) – The 2d DataArray containing surface level values. If None, ahn4-values are downloaded from the web. The default is None.

  • dx (float, optional) – The distance between the points along the line at which the surface level is calculated. Only used when num is None. When dx is None, it is set to the resolution of ahn. The default is None.

  • num (int, optional) – If not None, the surface level is calculated at num equally spaced points along the line. The default is None.

  • method (string, optional) – The method to interpolate the 2d surface level values to the points along the line. The default is “linear”.

  • plot (bool, optional) – if True, plot the 2d surface level, the line and the calculated heights. The default is False.

Returns:

z – A DataArray with dimension s, containing surface level values along the line.

Return type:

xr.DataArray

nlmod.read.ahn.download_ahn_at_point(x: float, y: float, buffer: float = 0.75, return_da: bool = False, return_mean: bool = False, identifier: str = 'dsm_05m', res: float = 0.5, **kwargs) float[source]

Get the height of the surface level at a certain point, defined by x and y.

Parameters:
  • x (float) – The x-coordinate fo the point.

  • y (float) – The y-coordinate fo the point..

  • buffer (float, optional) – The buffer around x and y that is downloaded. The default is 0.75.

  • return_da (bool, optional) – Return the downloaded DataArray when True. The default is False.

  • return_mean (bool, optional) – Resturn the mean of all non-nan pixels within buffer. Return the center pixel when False. The default is False.

  • identifier (str, optional) – The identifier passed onto download_latest_ahn_from_wcs. The default is “dsm_05m”.

  • res (float, optional) – The resolution that is passed onto download_latest_ahn_from_wcs. The default is 0.5.

  • **kwargs (dict) – kwargs are passed onto the method download_latest_ahn_from_wcs.

Returns:

The surface level value at the requested point.

Return type:

float

nlmod.read.ahn.download_latest_ahn_from_wcs(extent: list[float] = None, identifier: Literal['dsm_05m', 'dtm_05m'] = 'dsm_05m', res: float | None = None, version: str = '1.0.0', fmt: str = 'image/tiff', crs: str = 'EPSG:28992', maxsize: int = 2000, cachedir=None, cachename=None) DataArray[source]

Get the latest AHN from the wcs service.

Parameters:
  • extent (list, tuple or np.array, optional) – extent. The default is None.

  • identifier (str, optional) –

    Possible values for identifier are:

    ’dsm_05m’ ‘dtm_05m’

    The default is ‘dsm_05m’. the identifier contains resolution and type info: - ‘dtm’ is only surface level (maaiveld), ‘dsm’ has other surfaces such as buildings. - 5m or 05m is a resolution of 5x5 or 0.5x0.5 meter.

  • res (float, optional) – resolution of ahn raster. If None the resolution is inferred from the identifier. The default is None.

  • version (str, optional) – version of wcs service, options are ‘1.0.0’ and ‘2.0.1’. The default is ‘1.0.0’.

  • fmt (str, optional) – geotif format . The default is ‘GEOTIFF_FLOAT32’.

  • crs (str, optional) – coördinate reference system. The default is ‘EPSG:28992’.

  • maxsize (float, optional) – maximum number of cells in x or y direction. The default is 2000.

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Return type:

xr.DataArray

nlmod.read.ahn.get_ahn(ds: Dataset | None = None, identifier: str = None, method: str = 'average', extent: list[float] | None = None, merge_tiles: bool = True, cachedir=None, cachename=None, **kwargs) Dataset[source]

Get a model dataset with ahn variable.

Parameters:
  • ds (xr.Dataset, optional) – dataset with the model information. Using ds=None is deprecated, instead use nlmod.read.ahn.download_ahn().

  • method (str, optional) – Method used to resample ahn to grid of ds. See documentation of nlmod.resample.structured_da_to_ds for possible values. The default is ‘average’.

  • extent (list, tuple or np.array) – The extent to be downloaded, consisting of 4 floats: xmin, xmax, ymin, ymax. Only used if ds is None. The default is None.

  • version (str, optional) – The AHN_version, which can be ‘AHN2’, ‘AHN3’, ‘AHN4’, ‘AHN5’ and ‘AHN6’. version is ignored if identifier is specified. The default is “AHN4”, as this is available in the whole of the Netherlands.

  • data_kind (str, optional) – The kind of data. This can be ‘DTM’ (terrain elevation) or “DSM” (surface elevation). data_kind is ignored if identifier is specified. The default is “DTM”.

  • tile_size (str, optional) – The size of the map tiles that the data is offered by the webservice. This can be ‘5x6.25km’ or ‘1x1km’. tile_size is ignored if identifier is specified. The default is ‘5x6.25km’.

  • resolution (float, optional) – The resolution of the AHN-data, which can be 0.5 and 5.0 m. resolution is ignored if identifier is specified. The default is 5.0.

  • merge_tiles (bool, optional) – If True, the function returns a merged DataArray. If False, the function returns a list of DataArrays with the original tiles. The default is True.

  • cut_extent (bool, optional) – If True, only keep the requested extent from the data. The defualts is True.

  • config (dict, optional) – A dictionary with properties of the data sources of the different AHN-versions. When None, the configuration is retreived from the method get_configuration(). The default is None.

  • identifier (str, optional) –

    The identifier determines the AHN-version, the resolution and the type of height data. Possible values are (casing is important):

    AHN1: ‘AHN1 maaiveldmodel (DTM) 5m’ AHN2: ‘AHN2 maaiveldmodel (DTM) ½m, geïnterpoleerd’,

    ’AHN2 maaiveldmodel (DTM) ½m’, ‘AHN2 DSM ½m’, ‘AHN2 maaiveldmodel (DTM) 5m’

    AHN3: ‘AHN3 maaiveldmodel (DTM) ½m’,

    ’AHN3 DSM ½m’, ‘AHN3 maaiveldmodel (DTM) 5m’, ‘AHN3 DSM 5m’

    AHN4: ‘AHN4 maaiveldmodel (DTM) ½m’,

    ’AHN4 DSM ½m’, ‘AHN4 maaiveldmodel (DTM) 5m’, ‘AHN4 DSM 5m’

    AHN5: ‘AHN5 maaiveldmodel (DTM) 5m’,

    ’AHN5 DSM 5m’, ‘AHN5 maaiveldmodel (DTM) ½m’, ‘AHN5 DSM ½m’

    When no identifier is specified (the default), the url to download data from is taken from config, using the arguments version, data_kind, tile_size and resolution of this method. The default is None.

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

ds_out – Dataset with the ahn variable.

Return type:

xr.Dataset

nlmod.read.ahn.get_ahn1(extent: list[float], identifier: Literal['AHN1 maaiveldmodel (DTM) 5m'] = 'AHN1 maaiveldmodel (DTM) 5m', as_data_array: bool | None = None, cachedir=None, cachename=None, **kwargs) DataArray[source]

Download AHN1.

Deprecated since version 0.10.0: get_ahn1 will be removed in nlmod 1.0.0, it is replaced by download_ahn1 because of new naming convention https://github.com/gwmod/nlmod/issues/47

Parameters:
  • extent (list, tuple or np.array) – extent

  • identifier (str, optional) – The identifier determines the resolution and the type of height data. The only allowed value is ‘AHN1 maaiveldmodel (DTM) 5m’. The default is “AHN1 maaiveldmodel (DTM) 5m”.

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

DataArray of the AHN

Return type:

xr.DataArray

nlmod.read.ahn.get_ahn2(extent: list[float], identifier: Literal['AHN2 maaiveldmodel (DTM) ½m, geïnterpoleerd', 'AHN2 maaiveldmodel (DTM) ½m', 'AHN2 DSM ½m', 'AHN2 maaiveldmodel (DTM) 5m'] = 'AHN2 maaiveldmodel (DTM) 5m', as_data_array: bool | None = None, cachedir=None, cachename=None, **kwargs) DataArray[source]

Download AHN2.

Deprecated since version 0.10.0: get_ahn2 will be removed in nlmod 1.0.0, it is replaced by download_ahn2 because of new naming convention https://github.com/gwmod/nlmod/issues/47

Parameters:
  • extent (list, tuple or np.array) – extent

  • identifier (str, optional) – The identifier determines the resolution and the type of height data. Possible values are ‘AHN2 maaiveldmodel (DTM) ½m, geïnterpoleerd’, ‘AHN2 maaiveldmodel (DTM) ½m’, ‘AHN2 DSM ½m’, and ‘AHN2 maaiveldmodel (DTM) 5m’. The default is “AHN2 maaiveldmodel (DTM) 5m”.

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

DataArray of the AHN

Return type:

xr.DataArray

nlmod.read.ahn.get_ahn3(extent: list[float], identifier: Literal['AHN3 maaiveldmodel (DTM) ½m', 'AHN3 DSM ½m', 'AHN3 maaiveldmodel (DTM) 5m', 'AHN3 DSM 5m'] = 'AHN3 maaiveldmodel (DTM) 5m', as_data_array: bool | None = None, cachedir=None, cachename=None, **kwargs) DataArray[source]

Download AHN3.

Deprecated since version 0.10.0: get_ahn3 will be removed in nlmod 1.0.0, it is replaced by download_ahn3 because of new naming convention https://github.com/gwmod/nlmod/issues/47

Parameters:
  • extent (list, tuple or np.array) – extent

  • identifier (str, optional) – The identifier determines the resolution and the type of height data. Possible values are ‘AHN3 maaiveldmodel (DTM) ½m’, ‘AHN3 DSM ½m’, ‘AHN3 maaiveldmodel (DTM) 5m’, and ‘AHN3 DSM 5m’. The default is “AHN3 maaiveldmodel (DTM) 5m”.

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

DataArray of the AHN

Return type:

xr.DataArray

nlmod.read.ahn.get_ahn4(extent: list[float], identifier: Literal['AHN4 maaiveldmodel (DTM) ½m', 'AHN4 DSM ½m', 'AHN4 maaiveldmodel (DTM) 5m', 'AHN4 DSM 5m'] = 'AHN4 maaiveldmodel (DTM) 5m', as_data_array: bool | None = None, cachedir=None, cachename=None, **kwargs) DataArray[source]

Download AHN4.

Deprecated since version 0.10.0: get_ahn4 will be removed in nlmod 1.0.0, it is replaced by download_ahn4 because of new naming convention https://github.com/gwmod/nlmod/issues/47

Parameters:
  • extent (list, tuple or np.array) – extent

  • identifier (str, optional) – The identifier determines the resolution and the type of height data. Possible values are ‘AHN4 maaiveldmodel (DTM) ½m’, ‘AHN4 DSM ½m’, ‘AHN4 maaiveldmodel (DTM) 5m’, and ‘AHN4 DSM 5m’. The default is “AHN4 maaiveldmodel (DTM) 5m”.

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

DataArray of the AHN

Return type:

xr.DataArray

nlmod.read.ahn.get_ahn5(extent: list[float], identifier: Literal['AHN5 maaiveldmodel (DTM) 5m', 'AHN5 DSM 5m', 'AHN5 maaiveldmodel (DTM) ½m', 'AHN5 DSM ½m'] = 'AHN5 maaiveldmodel (DTM) 5m', cachedir=None, cachename=None, **kwargs) DataArray[source]

Download AHN5.

Deprecated since version 0.10.0: get_ahn5 will be removed in nlmod 1.0.0, it is replaced by download_ahn5 because of new naming convention https://github.com/gwmod/nlmod/issues/47

Parameters:
  • extent (list, tuple or np.array) – extent

  • identifier (str, optional) – The identifier determines the resolution and the type of height data. Possible values are ‘AHN5 maaiveldmodel (DTM) 5m’, ‘AHN5 DSM 5m’, ‘AHN5 maaiveldmodel (DTM) ½m’, and ‘AHN5 DSM ½m’. The default is “AHN5 maaiveldmodel (DTM) 5m”.

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

DataArray of the AHN

Return type:

xr.DataArray

nlmod.read.ahn.get_ahn_along_line(line: LineString, ahn: DataArray | None = None, dx: float | None = None, num: int | None = None, method: str = 'linear', plot: bool = False) DataArray[source]

Get the height of the surface level along a line.

Deprecated since version 0.10.0: get_ahn_along_line will be removed in nlmod 1.0.0, it is replaced by download_ahn_along_line because of new naming convention https://github.com/gwmod/nlmod/issues/47

Parameters:
  • line (shapely.LineString) – The line along which the surface level is calculated.

  • ahn (xr.DataArray, optional) – The 2d DataArray containing surface level values. If None, ahn4-values are downloaded from the web. The default is None.

  • dx (float, optional) – The distance between the points along the line at which the surface level is calculated. Only used when num is None. When dx is None, it is set to the resolution of ahn. The default is None.

  • num (int, optional) – If not None, the surface level is calculated at num equally spaced points along the line. The default is None.

  • method (string, optional) – The method to interpolate the 2d surface level values to the points along the line. The default is “linear”.

  • plot (bool, optional) – if True, plot the 2d surface level, the line and the calculated heights. The default is False.

Returns:

z – A DataArray with dimension s, containing surface level values along the line.

Return type:

xr.DataArray

nlmod.read.ahn.get_ahn_at_point(x: float, y: float, buffer: float = 0.75, return_da: bool = False, return_mean: bool = False, identifier: str = 'dsm_05m', res: float = 0.5, **kwargs) float[source]

Get the height of the surface level at a certain point, defined by x and y.

Deprecated since version 0.10.0: get_ahn_at_point will be removed in nlmod 1.0.0, it is replaced by download_ahn_at_point because of new naming convention https://github.com/gwmod/nlmod/issues/47

Parameters:
  • x (float) – The x-coordinate fo the point.

  • y (float) – The y-coordinate fo the point..

  • buffer (float, optional) – The buffer around x and y that is downloaded. The default is 0.75.

  • return_da (bool, optional) – Return the downloaded DataArray when True. The default is False.

  • return_mean (bool, optional) – Resturn the mean of all non-nan pixels within buffer. Return the center pixel when False. The default is False.

  • identifier (str, optional) – The identifier passed onto download_latest_ahn_from_wcs. The default is “dsm_05m”.

  • res (float, optional) – The resolution that is passed onto download_latest_ahn_from_wcs. The default is 0.5.

  • **kwargs (dict) – kwargs are passed onto the method download_latest_ahn_from_wcs.

Returns:

The surface level value at the requested point.

Return type:

float

nlmod.read.ahn.get_configuration()[source]

Get the location of json-files with positions of map tiles with AHN-data

Returns:

A dictionary with the locations of the json-files that contain the position of the map-tiles of AHN-data, for several combinations of version, data_kind, tile_size and resolution.

Return type:

dict

nlmod.read.ahn.get_latest_ahn_from_wcs(extent: list[float] = None, identifier: Literal['dsm_05m', 'dtm_05m'] = 'dsm_05m', res: float | None = None, version: str = '1.0.0', fmt: str = 'image/tiff', crs: str = 'EPSG:28992', maxsize: int = 2000, cachedir=None, cachename=None) DataArray[source]

Get the latest AHN from the wcs service.

Deprecated since version 0.10.0: get_latest_ahn_from_wcs will be removed in nlmod 1.0.0, it is replaced by download_latest_ahn_from_wcs because of new naming convention https://github.com/gwmod/nlmod/issues/47

Parameters:
  • extent (list, tuple or np.array, optional) – extent. The default is None.

  • identifier (str, optional) –

    Possible values for identifier are:

    ’dsm_05m’ ‘dtm_05m’

    The default is ‘dsm_05m’. the identifier contains resolution and type info: - ‘dtm’ is only surface level (maaiveld), ‘dsm’ has other surfaces such as buildings. - 5m or 05m is a resolution of 5x5 or 0.5x0.5 meter.

  • res (float, optional) – resolution of ahn raster. If None the resolution is inferred from the identifier. The default is None.

  • version (str, optional) – version of wcs service, options are ‘1.0.0’ and ‘2.0.1’. The default is ‘1.0.0’.

  • fmt (str, optional) – geotif format . The default is ‘GEOTIFF_FLOAT32’.

  • crs (str, optional) – coördinate reference system. The default is ‘EPSG:28992’.

  • maxsize (float, optional) – maximum number of cells in x or y direction. The default is 2000.

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Return type:

xr.DataArray

nlmod.read.bgt._update_bronhouder_json(fname)[source]

“Update the json-file within nlmod with bronhouder-names

Parameters:

fname (str) – The path and name of a file with bronhouder-names, downloaded from https://www.kadaster.nl/-/bgt-bronhoudercodes.

nlmod.read.bgt.download_bgt(extent, layer='waterdeel', cut_by_extent=True, make_valid=False, fname=None, geometry=None, remove_expired=True, add_bronhouder_names=True, timeout=1200, cachedir=None, cachename=None)[source]

Get geometries within an extent or polygon from the Basis Registratie Grootschalige Topografie (BGT)

Parameters:
  • extent (list or tuple of length 4 or shapely Polygon) – The extent (xmin, xmax, ymin, ymax) or polygon for which shapes are requested.

  • layer (string or list of strings, optional) – The layer(s) for which shapes are requested. When layer is “all”, all layers are requested. The default is “waterdeel”.

  • cut_by_extent (bool, optional) – Only return the intersection with the extent if True. The default is True

  • make_valid (bool, optional) – Make geometries valid by appying a buffer of 0 m when True. The default is False.

  • fname (string, optional) – Save the zipfile that is received by the request to file. The default is None, which does not save anything to file.

  • geometry (string, optional) – When geometry is specified, this attribute is used as the geometry of the resulting GeoDataFrame. Some layers have multiple geometry-attributes. An example is the layer ‘pand’, where each buidling (polygon) also contains a Point-geometry for the label. When geometry is None, the last attribute starting with the word “geometrie” is used as the geometry. The default is None.

  • remove_expired (bool, optional) – Remove expired items (that contain a value for ‘eindRegistratie’) when True. The default is True.

  • add_bronhouder_names (bool, optional) – Add bronhouder in a column called ‘bronhouder_name. names when True. The default is True.

  • timeout (int optional) – The amount of time in seconds to wait for the server to send data before giving up. The default is 1200 (20 minutes).

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

gdf – A GeoDataFrame (when only one layer is requested) or a dict of GeoDataFrames containing all geometries and properties.

Return type:

GeoPandas GeoDataFrame or dict of GeoPandas GeoDataFrame

nlmod.read.bgt.get_bgt(*args, **kwargs)[source]

Get geometries within an extent or polygon from the Basis Registratie Grootschalige Topografie (BGT)

Deprecated since version 0.10.0: get_bgt will be removed in nlmod 1.0.0, it is replaced by download_bgt because of new naming convention https://github.com/gwmod/nlmod/issues/47

Parameters:
  • extent (list or tuple of length 4 or shapely Polygon) – The extent (xmin, xmax, ymin, ymax) or polygon for which shapes are requested.

  • layer (string or list of strings, optional) – The layer(s) for which shapes are requested. When layer is “all”, all layers are requested. The default is “waterdeel”.

  • cut_by_extent (bool, optional) – Only return the intersection with the extent if True. The default is True

  • make_valid (bool, optional) – Make geometries valid by appying a buffer of 0 m when True. The default is False.

  • fname (string, optional) – Save the zipfile that is received by the request to file. The default is None, which does not save anything to file.

  • geometry (string, optional) – When geometry is specified, this attribute is used as the geometry of the resulting GeoDataFrame. Some layers have multiple geometry-attributes. An example is the layer ‘pand’, where each buidling (polygon) also contains a Point-geometry for the label. When geometry is None, the last attribute starting with the word “geometrie” is used as the geometry. The default is None.

  • remove_expired (bool, optional) – Remove expired items (that contain a value for ‘eindRegistratie’) when True. The default is True.

  • add_bronhouder_names (bool, optional) – Add bronhouder in a column called ‘bronhouder_name. names when True. The default is True.

  • timeout (int optional) – The amount of time in seconds to wait for the server to send data before giving up. The default is 1200 (20 minutes).

Returns:

gdf – A GeoDataFrame (when only one layer is requested) or a dict of GeoDataFrames containing all geometries and properties.

Return type:

GeoPandas GeoDataFrame or dict of GeoPandas GeoDataFrame

nlmod.read.bgt.get_bgt_layers(timeout=1200)[source]

Get the layers in the Basis Registratie Grootschalige Topografie (BGT)

Parameters:

timeout (int optional) – The amount of time in seconds to wait for the server to send data before giving up. The default is 1200 (20 minutes).

Returns:

A list with the layer names.

Return type:

list

nlmod.read.bgt.get_bronhouder_names() Dict[str, str][source]

Get the names of the bronhouders of the BGT data.

Returns:

A dictionary with the bronhouder codes as keys and the names as values.

Return type:

dict

Notes

The bronhouder names are retrieved from https://www.kadaster.nl/-/bgt-bronhoudercodes. The Toelichting sheet in the .ods file gives the changes compared to the old file. If changes are made, please add them manually to the bgt_bronhouder_names dictionary. A test is added for to check if the dictionary is up to date with the latest file from the Kadaster up to date with the .ods file from 2024-01-01.

nlmod.read.bgt.read_bgt_gml(fname, geometry='geometrie2dGrondvlak', crs='epsg:28992')[source]
nlmod.read.bgt.read_bgt_zipfile(fname, geometry=None, files=None, cut_by_extent=True, make_valid=False, extent=None, remove_expired=True, add_bronhouder_names=True)[source]

Read data from a zipfile that was downloaded using get_bgt().

Parameters:
  • fname (string) – The filename of the zip-file containing the BGT-data.

  • geometry (str, optional) – DESCRIPTION. The default is None.

  • files (string of list of strings, optional) – The files to read from the zipfile. Read all files when files is None. The default is None.

  • cut_by_extent (bool, optional) – Cut the geoemetries by the supplied extent. When no extent is supplied, cut_by_extent is set to False. The default is True.

  • make_valid (bool, optional) – Make geometries valid by appying a buffer of 0 m when True. THe defaults is False.

  • extent (list or tuple of length 4 or shapely Polygon) – The extent (xmin, xmax, ymin, ymax) or polygon by which the geometries are clipped. Only used when cut_by_extent is True. The defult is None.

  • remove_expired (bool, optional) – Remove expired items (that contain a value for ‘eindRegistratie’) when True. The default is True.

  • add_bronhouder_names (bool, optional) – Add bronhouder in a column called ‘bronhouder_name. names when True. The default is True.

Returns:

gdf – A dict of GeoDataFrames containing all geometries and properties.

Return type:

dict of GeoPandas GeoDataFrame

nlmod.read.geotop._get_kh_kv_from_df(df, ilithok, istrat=None, anisotropy=1.0, mask=None)[source]
nlmod.read.geotop._handle_nans_in_stochastic_approach(kh, kv, kh_method, kv_method)[source]
nlmod.read.geotop._save_excel_files_as_csv()[source]

This method takes the files REF_GTP_STR_UNIT.xlsx and REF_GTP_LITHO_CLASS.xlsx that are taken from the GeoTOP 1.6 zipfile downloaded from DINOloket, and saves them as csv-files. In this way version-control can better process the changes in future versions of GeoTOP.

Return type:

None.

nlmod.read.geotop.add_kh_and_kv(gt, df, stochastic=None, kh_method='arithmetic_mean', kv_method='harmonic_mean', anisotropy=1.0, kh='kh', kv='kv', kh_df='kh', kv_df='kv')[source]

Add kh and kv variables to the voxels of the GeoTOP Dataset.

Parameters:
  • gt (xr.Dataset) – The geotop dataset, at least with variable lithok.

  • df (pd.DataFrame) – A DataFrame with information about the kh and optionally kv, for different lithoclasses or stratigraphic units. The DataFrame must contain the columns ‘lithok’ and ‘kh’, and optionally ‘strat’ and ‘kv’. As an example see nlmod.read.geotop.get_kh_kv_table().

  • stochastic (bool, str or None, optional) – When stochastic is True or a string, use the stochastic data of GeoTOP. The only supported method right now is “linear”, which means kh and kv are determined from a linear weighted mean of the voxels. For kh the method from kh_method is used to calculate the mean. For kv the method from kv_method is used. When stochastic is False or None, the stochastic data of GeoTOP is not used. The default is None.

  • kh_method (str, optional) – The method to calculate the weighted mean of kh values when stochastic is True or “linear”. Allowed values are “arithmetic_mean” and “harmonic_mean”. The default is “arithmetic_mean”.

  • kv_method (str, optional) – The method to calculate the weighted mean of kv values when stochastic is True or “linear”. Allowed values are “arithmetic_mean” and “harmonic_mean”. The default is “arithmetic_mean”.

  • anisotropy (float, optional) – The anisotropy value used when there are no kv values in df. The default is 1.0.

  • kh (str, optional) – THe name of the new variable with kh values in gt. The default is “kh”.

  • kv (str, optional) – The name of the new variable with kv values in gt. The default is “kv”.

  • kh_df (str, optional) – The name of the column with kh values in df. The default is “kh”.

  • kv_df (str, optional) – THe name of the column with kv values in df. The default is “kv”.

Raises:

DESCRIPTION.

Returns:

gt – Datset with voxel-data, with the added variables ‘kh’ and ‘kv’.

Return type:

xr.Dataset

nlmod.read.geotop.add_top_and_botm(ds)[source]

Add the top and bottom of the voxels to the GeoTOP Dataset.

This makes sure the structure of the geotop dataset is more like regis, and we can use the cross-section class (DatasetCrossSection from nlmod.

Parameters:

ds (xr.Dataset) – The geotop-dataset.

Returns:

ds – The geotop-dataset, with added variables “top” and “botm”.

Return type:

xr.Dataset

nlmod.read.geotop.aggregate_to_ds(gt, ds, kh='kh', kv='kv', kd='kD', c='c', kh_gt='kh', kv_gt='kv', add_kd_and_c=False)[source]

Aggregate voxels from GeoTOP to layers in a model DataSet with top and botm, to calculate kh and kv.

Parameters:
  • gt (xr.Dataset) – A Dataset containing the Geotop voxel data.

  • ds (xr.Dataset) – A Dataset containing the top and botm of the layers.

  • kh (str, optional) – The name of the new variable for the horizontal conductivity in ds. The default is “kh”.

  • kv (str, optional) – The name of the new variable for the vertical conductivity in ds. The default is “kv”.

  • kd (str, optional) – The name of the variable for the horizontal transmissivity. Only used when add_kd_and_c is True. The default is “kD”.

  • c (str, optional) – The name of the variable for the vertical reistance. Only used when add_kd_and_c is True The default is “c”.

  • kh_gt (str, optional) – The name of the variable for the horizontal conductivity in gt. The default is “kh”.

  • kv_gt (str, optional) – The name of the variable for the vertical conductivity in gt. The default is “kv”.

  • add_kd_and_c (bool, optional) – Add the variables kd and c to ds. The default is False.

Returns:

ds – The Dataset ds, with added variables kh and kv (and optionally kd and c).

Return type:

xr.Dataset

nlmod.read.geotop.download_geotop(extent, url=None, probabilities=False, chunks='auto', cachedir=None, cachename=None)[source]

Get a slice of the geotop netcdf url within the extent, set the x and y coordinates to match the cell centers and keep only the strat and lithok data variables.

Parameters:
  • extent (list, tuple or np.array) – desired model extent (xmin, xmax, ymin, ymax)

  • url (str, optional) – url of geotop netcdf file. The default is nlmod.read.geotop.GEOTOP_URL: http://www.dinodata.nl/opendap/GeoTOP/geotop.nc

  • probabilities (bool, optional) – if True, also download probability data. The default is False.

  • chunks (int, dict, 'auto' or None. The default is 'auto'.) –

    If provided, used to load the data into dask arrays. - chunks="auto" will use dask auto chunking. - chunks=None skips using dask. This uses xarray’s internally private lazy

    indexing classes, but data is eagerly loaded into memory as numpy arrays when accessed. This can be more efficient for smaller arrays or when large arrays are sliced before computation.

    See dask chunking for more details.

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

gt – slices geotop netcdf.

Return type:

xarray Dataset

nlmod.read.geotop.get_geotop(*args, **kwargs)[source]

Get a slice of the geotop netcdf url within the extent, set the x and y coordinates to match the cell centers and keep only the strat and lithok data variables.

Deprecated since version 0.10.0: get_geotop will be removed in nlmod 1.0.0, it is replaced by download_geotop because of new naming convention https://github.com/gwmod/nlmod/issues/47

Parameters:
  • extent (list, tuple or np.array) – desired model extent (xmin, xmax, ymin, ymax)

  • url (str, optional) – url of geotop netcdf file. The default is nlmod.read.geotop.GEOTOP_URL: http://www.dinodata.nl/opendap/GeoTOP/geotop.nc

  • probabilities (bool, optional) – if True, also download probability data. The default is False.

Returns:

gt – slices geotop netcdf.

Return type:

xarray Dataset

nlmod.read.geotop.get_kh_kv_table(kind='Brabant')[source]
nlmod.read.geotop.get_lithok_colors()[source]
nlmod.read.geotop.get_lithok_props(rgb_colors=True)[source]
nlmod.read.geotop.get_strat_props()[source]
nlmod.read.geotop.split_layers_on_geul(strat, units_no_geul, geulen)[source]

Modifies the stratigraphic data from geotop in such a way that every stratigraphic unit is completely above or completely below any other unit (and never both above and below the same unit). This is useful for creating a layer model from the stratigraphic units.

This function splits the geulen over multiple layers and splits other layers locally when a geul is partly above and partly below that layer.

Some extra logic is added to minimize the number of geul layers by finding the most efficient layer to add the geul to.

Parameters:
  • strat (xarray DataArray) – with dimensions z, y and x, with values corresponding to the stratigraphic unit in each voxel.

  • units_no_geul (list) – Ordered list of stratigraphic units without geulen.

  • geulen (list) – Ordered list of geulen units.

Returns:

  • strat (xarray DataArray) – Modified stratigraphic data, including all the new geul units and the non-geul units that are split because the geul was in between them.

  • new_unit_order (list) – Ordered list with new stratigraphic units, including all the new geul units and the non-geul units that are split because the geul was in between them.

Notes

the new stratigraphic data contains more unit numbers than the original stratigraphic data. When a layer is split by a geul, the part of the layer below the geul gets a new unit number. The same goes for the geul itself when it occurs across multiple layers.

nlmod.read.geotop.to_model_layers(geotop_ds, strat_props=None, method_geulen='add_to_layer_below', drop_layer_dim_from_top=True, cachedir=None, cachename=None, **kwargs)[source]

Convert geotop voxel dataset to layered dataset.

Converts geotop data to dataset with layer elevations and hydraulic conductivities. Optionally uses hydraulic conductivities provided present in geotop_ds.

Parameters:
  • geotop_ds (xr.DataSet) – geotop voxel dataset (download using get_geotop(extent))

  • strat_props (pd.DataFrame, optional) – The properties (code and name) of the stratigraphic units. Load with get_strat_props() when None. The default is None.

  • method_geulen (str, optional) – strat-units >=6000 are so-called ‘geulen’ (paleochannels, gullies). These are difficult to add to the layer model, because they can occur above and/or below any other unit. Multiple methods are available to handle these ‘geulen’. The method “add_to_layer_below” adds the thickness of the ‘geul’ to the layer with a positive thickness below the ‘geul’. The method “add_to_layer_above” adds the thickness of the ‘geul’ to the layer with a positive thickness above the ‘geul’. The method “add_as_layer” tries to add the ‘geulen’ as one or more layers, which can fail if a ‘geul’ is locally both below the top and above the bottom of another layer (splitting the layer in two, which is not supported). The method “split_layers” splits layers when a ‘geul’ is both below the top and above the bottom of another layer. The default is “add_to_layer_below”.

  • drop_layer_dim_from_top (bool, optional) – When True, fill NaN values in top and botm and drop the layer dimension from top. This will transform top and botm to the data model in MODFLOW. An advantage of this data model is that the layer model is consistent by definition, with no possibilities of gaps between layers. The default is True.

  • kwargs (dict) – Kwargs are passed to aggregate_to_ds()

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

ds – dataset with top and botm (and optionally kh and kv) per geotop layer

Return type:

xr.DataSet

Module with functions to deal with the northsea.

  • identify model cells with the north sea

  • add bathymetry of the northsea to the layer model

  • extrapolate the layer model below the northsea bed.

Note: if you like jazz please check this out: https://www.northseajazz.com

nlmod.read.jarkus.discretize_bathymetry(ds, bathymetry_da, da_name='bathymetry', datavar_sea='northsea', method='average', cachedir=None, cachename=None)[source]

Discretize bathymetry data to model the model grid.

Parameters:
  • ds (xarray.Dataset) – dataset with model data where bathymetry is added to.

  • bathymetry_da (xarray.DataArray) – bathymetry data

  • da_name (str, optional) – name of the datavar that is used to store the bathymetry data. The default is ‘bathymetry’.

  • datavar_sea (str, optional) – datavariable in the dataset that is used to identify cells with sea in them. The default is ‘northsea’.

  • method (str, optional) – Method used to resample bathymetry data to the modelgrid. See the documentation of nlmod.resample.structured_da_to_ds for possible values. The default is ‘average’.

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

ds_out – dataset with bathymetry

Return type:

xarray.Dataset

Notes

The nan values in the original bathymetry are filled and then the data is resampled to the modelgrid. Maybe we can speed up things by changing the order in which operations are executed.

nlmod.read.jarkus.download_bathymetry(extent, kind='jarkus', cachedir=None, cachename=None)[source]

Download bathymetry data from the jarkus dataset.

Parameters:
  • extent (list, tuple or np.array) – extent xmin, xmax, ymin, ymax.

  • kind (str, optional) – The kind of data. Can be “jarkus”, “kusthoogte” or “vaklodingen”. The default is “jarkus”.

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

da_bathymetry – bathymetry data

Return type:

xarray.DataArray

nlmod.read.jarkus.get_bathymetry(ds=None, extent=None, da_name='bathymetry', datavar_sea='northsea', kind='jarkus', method='average', cachedir=None, cachename=None)[source]

Get bathymetry of the Northsea from the jarkus dataset.

Parameters:
  • ds (xarray.Dataset) – dataset with model data where bathymetry is added to. Using ds=None is deprecated, instead use nlmod.read.jarkus.download_bathymetry().

  • extent (list, tuple or np.array, optional) – extent xmin, xmax, ymin, ymax. Only used if ds is None. The default is None.

  • da_name (str, optional) – name of the datavar that is used to store the bathymetry data. The default is ‘bathymetry’.

  • datavar_sea (str, optional) – datavariable in the dataset that is used to identify cells with sea in them. The default is ‘northsea’.

  • method (str, optional) – Method used to resample bathymetry data to the modelgrid. See the documentation of nlmod.resample.structured_da_to_ds for possible values. The default is ‘average’.

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

ds_out – dataset with bathymetry

Return type:

xarray.Dataset

Notes

The nan values in the original bathymetry are filled and then the data is resampled to the modelgrid. Maybe we can speed up things by changing the order in which operations are executed.

nlmod.read.jarkus.get_dataset_jarkus(extent, kind='jarkus', return_tiles=False, time=-1, cachedir=None, cachename=None)[source]

Get bathymetry from Jarkus within a certain extent. If return_tiles is False, the following actions are performed: 1. find Jarkus tiles within the extent 2. download netcdf files of Jarkus tiles 3. read Jarkus tiles and combine the ‘z’ parameter of the last time step of each tile (when time=1), to a dataarray.

Parameters:
  • extent (list, tuple or np.array) – extent (xmin, xmax, ymin, ymax) of the desired grid. Should be RD-new coordinates (EPSG:28992)

  • kind (str, optional) – The kind of data. Can be “jarkus”, “kusthoogte” or “vaklodingen”. The default is “jarkus”.

  • return_tiles (bool, optional) – Return the individual tiles when True. The default is False.

  • time (str, int or pd.TimeStamp, optional) – The time to return data for. When time=”last_non_nan”, this returns the last non-NaN-value for each pixel. This can take a while, as all tiles need to be checked. When time is an integer, it is used as the time index. When set to -1, this then downloads the last time available in each tile (which can contain large areas with NaN-values). When time is a string (other than “last_non_nan”) or a pandas Timestamp, only data on this exact time are downloaded. The default is -1.

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

z – dataset containing bathymetry data

Return type:

xr.DataSet

nlmod.read.jarkus.get_jarkus_tilenames(extent, kind='jarkus')[source]

Find all Jarkus tilenames within a certain extent.

Parameters:

extent (list, tuple or np.array) – extent (xmin, xmax, ymin, ymax) of the desired grid. Should be RD-new coordinates (EPSG:28992)

Returns:

netcdf_urls – list of the urls of all netcdf files of the tiles with Jarkus data.

Return type:

list of str

nlmod.read.jarkus.get_netcdf_tiles(kind='jarkus')[source]

Find all Jarkus netcdf tile names.

Returns:

netcdf_urls – list of the urls of all netcdf files of the tiles with Jarkus data.

Return type:

list of str

Notes

This function would be redundant if the jarkus catalog (http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/jarkus/grids/catalog.nc) had a proper way of displaying the url’s of each tile. It seems like an attempt was made to do this because there is a data variable named ‘urlPath’ in the catalog. However the dataarray of ‘urlPath’ has the same string for each tile.

nlmod.read.knmi._add_ts_to_ds(timeseries, loc_sel, variable, ds)[source]

Add a timeseries to a variable at location loc_sel in model DataSet.

nlmod.read.knmi._download_knmi_at_locations(locations, start=None, end=None)[source]

Get precipitation (RD) and evaporation (EV24) data from the knmi at the locations

Parameters:
  • locations (pd.DataFrame) – each row contains a location (x and y) and the relevant precipitation (stn_RD) and evaporation (stn_EV24) stations.

  • start (str or datetime, optional) – start date of measurements that you want, The default is ‘2010’.

  • end (str or datetime, optional) – end date of measurements that you want, The default is None.

Returns:

hpd.ObsCollection

Return type:

oc_knmi

nlmod.read.knmi._get_locations_structured(ds)[source]

Get dataframe with the locations of the grid cells of a structured grid.

Parameters:

ds (xr.DataSet) – dataset containing relevant model grid information

Returns:

locations – DataFrame with the locations of all active grid cells. includes the columns: x, y, row, col and layer

Return type:

pandas DataFrame

nlmod.read.knmi._get_locations_vertex(ds)[source]

Get dataframe with the locations of the grid cells of a vertex grid.

Parameters:

ds (xr.DataSet) – dataset containing relevant model grid information

Returns:

locations – DataFrame with the locations of all active grid cells. includes the columns: x, y and layer

Return type:

pandas DataFrame

nlmod.read.knmi._resample_df_to_model_time(df, ds)[source]
nlmod.read.knmi.discretize_knmi(ds, oc_knmi, method='linear', most_common_station=True, add_stn_dimensions=False, to_model_time=True, hourly_precision=None, cachedir=None, cachename=None)[source]

discretize knmi data to model grid

Create a dataset with recharge (and evaporation) data by following these steps:
  1. Check for each cell (structured or vertex) which knmi measurement stations (prec and evap) are the closest.

  2. Download precipitation and evaporation data for all knmi stations that were found at 1

  3. Create a recharge array in which each cell has a reference to a timeseries. Timeseries are created for each unique combination of precipitation and evaporation. The following packages are created:

    1. the rch package itself in which cells with the same precipitation and evaporation stations are defined. This package also refers to all the time series package (see b).

    2. the time series packages in which the recharge flux is defined for the time steps of the model. Each package contains the time series for one or more cels (defined in a).

Parameters:
  • ds (xr.DataSet) – dataset containing relevant model grid information

  • oc_knmi (hpd.ObsCollection) – ObsCollection with precipitation (RD) and evaporation (EV24) data from the knmi.

  • method (str, optional) – If ‘linear’, calculate recharge by subtracting evaporation from precipitation. If ‘separate’, add precipitation as ‘recharge’ and evaporation as ‘evaporation’. Method is only used when add_stn_dimensions=False. The default is ‘linear’.

  • most_common_station (bool, optional) – When True, only use data from the station that is most common in the model area. The default is True

  • add_stn_dimensions (bool, optional) – When True, add the dimension time to the variable recharge (and evaporation when method=’seperate’). When True, add dimension stn_RD and stn_EV24 to the variable recharge and evaporation, and add variables “recharge_stn” and “evaporation_stn” that specify for every grid cell which KNMI-stations are used. When add_stn_dimensions is False, specify recharge (and evaporation when method=’seperate’) for every gridcell. The default is False.

  • to_model_time (bool, optional) – When True, resample the recharge and evaporation to the dimension time in ds. When False, save the original times of the KNMI-data in variables time_RD and time_EV24. to_model_time=False is only supported for add_stn_dimensions=True. The default is True.

  • hourly_precision (bool, optional) – When True, take into account the hour of the daily knmi-measurements (1:00 for evaporation and 9:00 for precipitation). The default is None, which will raise a warning, mentioning that hourly_precision is set to False, but the default will change to True in the future.

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

ds – dataset with ‘recharge’ (and ‘evaporation’) variables.

Return type:

xr.DataSet

Raises:
  • AttributeError – if the dataset has no time dimension

  • TypeError – if time dimension does not have a datetime64[ns] type

  • ValueError – if grid is not vertex

nlmod.read.knmi.download_knmi(ds=None, extent=None, delr=None, delc=None, start=None, end=None, most_common_station=False, cachedir=None, cachename=None)[source]

Get precipitation (RD) and evaporation (EV24) data from the knmi at the grid cells.

Parameters:
  • ds (xr.DataSet) – dataset containing relevant model grid information

  • most_common_station (bool, optional) – When True, only download data from the station that is most common in the model area. The default is False

  • start (str or datetime, optional) – start date of measurements that you want, The default is ‘2010’.

  • end (str or datetime, optional) – end date of measurements that you want, The default is None.

Raises:

ValueError – If data is not available for the entire model period

:raises cachedir : str or None, optional: directory to save cache. If None no cache is used. Default is None. :raises cachename : str or None, optional: filename of netcdf cache. If None no cache is used. Default is None.

Returns:

hpd.ObsCollection

Return type:

oc_knmi

nlmod.read.knmi.get_knmi(ds, most_common_station=False, start=None, end=None)[source]

Get precipitation (RD) and evaporation (EV24) data from the knmi at the grid cells.

Deprecated since version 0.10.0: get_knmi will be removed in nlmod 1.0.0, it is replaced by download_knmi because of new naming convention https://github.com/gwmod/nlmod/issues/47

Parameters:
  • ds (xr.DataSet) – dataset containing relevant model grid information

  • most_common_station (bool, optional) – When True, only download data from the station that is most common in the model area. The default is False

  • start (str or datetime, optional) – start date of measurements that you want, The default is ‘2010’.

  • end (str or datetime, optional) – end date of measurements that you want, The default is None.

Returns:

hpd.ObsCollection

Return type:

oc_knmi

nlmod.read.knmi.get_knmi_at_locations(locations, ds=None, start=None, end=None)[source]

Get precipitation (RD) and evaporation (EV24) data from the knmi at the locations

Deprecated since version 0.10.0: get_knmi_at_locations will be removed in nlmod 1.0.0, it is replaced by _download_knmi_at_locations because of new naming convention https://github.com/gwmod/nlmod/issues/47

Parameters:
  • locations (pd.DataFrame) – each row contains a location (x and y) and the relevant precipitation (stn_RD) and evaporation (stn_EV24) stations.

  • ds (xr.DataSet or None, optional) – dataset containing relevant time information. If None provide start and end.

  • start (str or datetime, optional) – start date of measurements that you want, The default is ‘2010’.

  • end (str or datetime, optional) – end date of measurements that you want, The default is None.

Returns:

hpd.ObsCollection

Return type:

oc_knmi

nlmod.read.knmi.get_locations(ds, oc_knmi=None, most_common_station=False)[source]

Get the locations of the active grid cells in ds and the nearest (or most common) precipitation and evaporation station.

Parameters:
  • ds (xr.DataSet) – dataset containing relevant model grid information

  • oc_knmi (hpd.ObsCollection or None, optional) – ObsCollection with knmi station data. If None the nearest of all knmi stations is used.

  • most_common_station (bool, optional) – When True, only download data from the station that is most common in the model area. The default is False

Raises:

ValueError – wrong grid type specified.

Returns:

locations – each row contains a location (x and y) and the relevant precipitation (stn_RD) and evaporation (stn_EV24) stations.

Return type:

pd.DataFrame

nlmod.read.knmi.get_recharge(ds, oc_knmi=None, method='linear', most_common_station=False, add_stn_dimensions=False, to_model_time=True, hourly_precision=None, cachedir=None, cachename=None)[source]

Add recharge to model dataset from KNMI data.

Add recharge to the model dataset with knmi data by following these steps:
  1. check for each cell (structured or vertex) which knmi measurement stations (prec and evap) are the closest.

  2. download precipitation and evaporation data for all knmi stations that were found at 1

  3. create a recharge array in which each cell has a reference to a timeseries. Timeseries are created for each unique combination of precipitation and evaporation. The following packages are created:

    1. the rch package itself in which cells with the same precipitation and evaporation stations are defined. This package also refers to all the time series package (see b).

    2. the time series packages in which the recharge flux is defined for the time steps of the model. Each package contains the time series for one or more cels (defined in a).

Supports structured and unstructred datasets.

Parameters:
  • ds (xr.DataSet) – dataset containing relevant model grid information

  • oc_knmi (hpd.ObsCollection) – ObsCollection with precipitation (RD) and evaporation (EV24) data from the knmi.

  • method (str, optional) – If ‘linear’, calculate recharge by subtracting evaporation from precipitation. If ‘separate’, add precipitation as ‘recharge’ and evaporation as ‘evaporation’. Method is only used when add_stn_dimensions=False. The default is ‘linear’.

  • most_common_station (bool, optional) – When True, only download data from the station that is most common in the model area. The default is False

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

ds – dataset with spatial model data including the rch raster

Return type:

xr.DataSet

nlmod.read.regis.add_geotop_to_regis_layers(rg, gt, layers='HLc', geotop_k=None, remove_nan_layers=True, anisotropy=1.0, gt_layered=None)[source]

Combine geotop and regis in such a way that the one or more layers in Regis are replaced by the geo_eenheden of geotop.

Parameters:
  • rg (xarray.DataSet) – regis dataset

  • gt (xarray.DataSet) – geotop dataset

  • layers (str or list of strings) – The regis layers to be replaced by geotop-layers

  • geotop_k (pd.DataFrame, optional) – The DataFrame with information about kh and kv of the GeoTOP-data. This DataFrame must at least contain columns ‘lithok’ and ‘kh’.

  • remove_nan_layers (bool, optional) – When True, layers with only 0 or NaN thickness are removed. The default is True.

  • anisotropy (float, optional) – The anisotropy value (kh/kv) used when there are no kv values in df. The default is 1.0.

  • gt_layered (xarray.Dataset) – A layered representation of the geotop-dataset. By supplying this parameter, the user can change the GeoTOP-layering, which is usueally defined by nlmod.read.geotop.to_model_layers(gt).

Returns:

gt – combined dataset

Return type:

xr.DataSet

nlmod.read.regis.download_regis(extent, botm_layer='AKc', variables=('top', 'botm', 'kh', 'kv'), remove_nan_layers=True, drop_layer_dim_from_top=True, probabilities=False, nodata=-9999, rename_layers_to_version_2_2_2=True, cachedir=None, cachename=None)[source]

Download a regis dataset within an extent.

Parameters:
  • extent (list, tuple or np.array) – desired model extent (xmin, xmax, ymin, ymax)

  • botm_layer (str, optional) – regis layer that is used as the bottom of the model. This layer is included in the model. the Default is “AKc” which is the bottom layer of regis. call nlmod.read.regis.get_layer_names() to get a list of regis names.

  • variables (tuple or list, optional) – The variables to keep from the regis Dataset. Possible entries in the list are ‘top’, ‘botm’, ‘kD’, ‘c’, ‘kh’, ‘kv’, ‘sdh’ and ‘sdv’. The default is (“top”, “botm”, “kh”, “kv”).

  • remove_nan_layers (bool, optional) – When True, layers that do not occur in the requested extent (layers that contain only NaN values for the botm array) are removed. The default is True.

  • drop_layer_dim_from_top (bool, optional) – When True, fill NaN values in top and botm and drop the layer dimension from top. This will transform top and botm to the data model in MODFLOW. An advantage of this data model is that the layer model is consistent by definition, with no possibilities of gaps between layers. The default is True.

  • probabilities (bool, optional) – if True, also download probability data. The default is False.

  • nodata (int or float, optional) – When nodata is not None, set values equal to nodata to nan. The default is -9999.

  • rename_layers_to_version_2_2_2 (bool, toptional) – From version 2.2.3 of regis, the names of stratigraphic layers change, compared to previous versions. If rename_layers_to_version_2_2_2 is True, the layer-names are renamed back to their original names. The default is True.

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

regis_ds – dataset with regis data projected on the modelgrid.

Return type:

xarray dataset

nlmod.read.regis.extract_version_from_title(version_string)[source]

Extract version number in format X.Y.Z from a string like “REGIS vXXrYsZ”.

Parameters:

version_string (str) – The input string containing version information in format “REGIS vXXrYsZ”.

Returns:

Extracted version in format “X.Y.Z”.

Return type:

packaging.version.Version

Examples

>>> extract_version("REGIS v02r2s3")
<Version('2.2.3')>
nlmod.read.regis.get_combined_layer_models(extent, regis_botm_layer='AKc', use_regis=True, use_geotop=True, remove_nan_layers=True, geotop_layers='HLc', geotop_k=None, gt_layered=None, cachedir=None, cachename=None)[source]

Combine layer models into a single layer model.

Possibilities so far include:
  • use_regis -> full model based on regis

  • use_regis and use_geotop -> holoceen of REGIS is filled with geotop

Parameters:
  • extent (list, tuple or np.array) – desired model extent (xmin, xmax, ymin, ymax)

  • regis_botm_layer (binary str, optional) – regis layer that is used as the bottom of the model. This layer is included in the model. the Default is ‘AKc’ which is the bottom layer of regis. call nlmod.read.regis.get_layer_names() to get a list of regis names.

  • use_regis (bool, optional) – True if part of the layer model should be REGIS. The default is True.

  • use_geotop (bool, optional) – True if part of the layer model should be geotop. The default is True.

  • remove_nan_layers (bool, optional) – When True, layers which contain only NaNs for the botm array are removed. The default is True.

  • geotop_layers (str or list of strings) – The regis layers to be replaced by geotop-layers

  • geotop_k (pd.DataFrame, optional) – The DataFrame with information about kh and kv of the GeoTOP-data. This DataFrame must at least contain columns ‘lithok’ and ‘kh’.

  • gt_layered (xarray.Dataset) – A layered representation of the geotop-dataset. By supplying this parameter, the user can change the GeoTOP-layering, which is usually defined by nlmod.read.geotop.to_model_layers(gt).

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

combined_ds – combination of layer models.

Return type:

xarray dataset

Raises:

ValueError – if an invalid combination of layers is used.

nlmod.read.regis.get_layer_names(ds=None, rename_layers_to_version_2_2_2=True)[source]

Get all the available regis layer names.

Parameters:
  • ds (xarray.Dataset, optional) – The regis dataset. If None, a connection is made to the REGIS server. The default is None.

  • rename_layers_to_version_2_2_2 (bool, optional) – If True, the layer names are renamed to their pre-v2.2.3 names. The default is True.

Returns:

layer_names – array with names of all the regis layers.

Return type:

np.array

nlmod.read.regis.get_legend(kind='REGIS')[source]

Get a legend (DataFrame) with the colors of REGIS and/or GeoTOP layers.

These colors can be used when plotting cross-sections.

nlmod.read.regis.get_legend_lithoclass()[source]
nlmod.read.regis.get_legend_lithostratigraphy()[source]
nlmod.read.regis.get_regis(extent, botm_layer='AKc', variables=('top', 'botm', 'kh', 'kv'), remove_nan_layers=True, drop_layer_dim_from_top=True, probabilities=False, nodata=-9999, rename_layers_to_version_2_2_2=True, cachedir=None, cachename=None)[source]

Get a regis dataset projected on the modelgrid.

Deprecated since version 0.10.0: get_regis will be removed in nlmod 1.0.0, it is replaced by download_regis because of new naming convention https://github.com/gwmod/nlmod/issues/47

Parameters:
  • extent (list, tuple or np.array) – desired model extent (xmin, xmax, ymin, ymax)

  • botm_layer (str, optional) – regis layer that is used as the bottom of the model. This layer is included in the model. the Default is “AKc” which is the bottom layer of regis. call nlmod.read.regis.get_layer_names() to get a list of regis names.

  • variables (tuple or list, optional) – The variables to keep from the regis Dataset. Possible entries in the list are ‘top’, ‘botm’, ‘kD’, ‘c’, ‘kh’, ‘kv’, ‘sdh’ and ‘sdv’. The default is (“top”, “botm”, “kh”, “kv”).

  • remove_nan_layers (bool, optional) – When True, layers that do not occur in the requested extent (layers that contain only NaN values for the botm array) are removed. The default is True.

  • drop_layer_dim_from_top (bool, optional) – When True, fill NaN values in top and botm and drop the layer dimension from top. This will transform top and botm to the data model in MODFLOW. An advantage of this data model is that the layer model is consistent by definition, with no possibilities of gaps between layers. The default is True.

  • probabilities (bool, optional) – if True, also download probability data. The default is False.

  • nodata (int or float, optional) – When nodata is not None, set values equal to nodata to nan. The default is -9999.

  • rename_layers_to_version_2_2_2 (bool, toptional) – From version 2.2.3 of regis, the names of stratigraphic layers change, compared to previous versions. If rename_layers_to_version_2_2_2 is True, the layer-names are renamed back to their original names. The default is True.

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

regis_ds – dataset with regis data projected on the modelgrid.

Return type:

xarray dataset

nlmod.read.regis.get_table_name_changes()[source]

Get the table with name changes of REGIS.

Returns:

df – A DataFrame containsing old and new names.

Return type:

pd.DataFrame

nlmod.read.regis.read_gleg(fname)[source]
nlmod.read.regis.read_voleg(fname)[source]
nlmod.read.rws.calculate_sea_coverage(dtm, ds=None, zmax=0.0, xy_sea=None, diagonal=False, method='mode', nodata=-1, return_filled_dtm=False)[source]

Determine where the sea is by interpreting the digital terrain model.

This method assumes the pixel defined in xy_sea (by default top-left) of the DTM-DataArray is sea. It then determines the height of the sea that is required for other pixels to become sea as well, taking into account the pixels in between.

Parameters:
  • dtm (xr.DataArray) – The digital terrain data, which can be of higher resolution than ds, Nans are filled by the minial value of dtm.

  • ds (xr.Dataset, optional) – Dataset with model information. When ds is not None, the sea DataArray is transformed to the model grid. The default is None.

  • zmax (float, optional) – Locations thet become sea when the sea level reaches a level of zmax will get a value of 1 in the resulting DataArray. The default is 0.0.

  • xy_sea (tuble of 2 floats) – The x- and y-coordinate of a location within the dtm that is sea. From this point, calculate_sea determines at what level each cell becomes wet. When xy_cell is None, the most northwest grid cell is sea, which is appropriate for the Netherlands. The default is None.

  • diagonal (bool, optional) – When true, dtm-values are connected diagonally as well (to determine the level the sea will reach). The default is False.

  • method (str, optional) – The method used to scale the dtm to ds. The default is “mode” (mode means that if more than half of the (not-nan) cells are wet, the cell is classified as sea).

  • nodata (int or float, optional) – The value for model cells outside the coverage of the dtm. Only used internally. The default is -1.

  • return_filled_dtm (bool, optional) – When True, return the filled dtm. The default is False.

Returns:

sea – A DataArray with value of 1 where the sea is and 0 where it is not.

Return type:

xr.DataArray

nlmod.read.rws.discretize_northsea(ds, gdf=None, da_name='northsea', cachedir=None, cachename=None)[source]

Get Dataset which is 1 at the northsea and 0 everywhere else.

Sea is defined by rws surface water shapefile.

Parameters:
  • ds (xr.DataSet, None, optional) – dataset containing relevant model information

  • gdf (gpd.GeoDataFrame or None, optional) – geometries of the surface water, if None the geometries are obtained using get_gdf_surface_water. The default is None.

  • da_name (str, optional) – name of the datavar that identifies sea cells

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

ds_out – Dataset with a single DataArray, this DataArray is 1 at sea and 0 everywhere else. Grid dimensions according to ds.

Return type:

xr.DataSet

nlmod.read.rws.discretize_surface_water(ds, gdf, da_basename='rws_oppwater', cachedir=None, cachename=None)[source]

Create 3 data-arrays from the shapefile with surface water.

These data arrays are: - area: area of the shape in the cell - cond: conductance based on the area and “bweerstand” column in shapefile - stage: surface water level based on the “peil” column in the shapefile

Parameters:
  • ds (xr.DataSet) – xarray with model data

  • gdf (gpd.GeoDataFrame or None, optional) – geometries of the surface water and the columns ‘bweerstand’ and ‘peil’.

  • da_basename (str) – name of the polygon shapes, name is used as a prefix to store data arrays in ds

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

ds – dataset with modelgrid data.

Return type:

xarray.Dataset

nlmod.read.rws.download_bathymetry(extent: list[float], resolution: Literal['20m', '1m'] = '20m', res: float | None = None, method: str | Callable | None = None, chunks: str | dict[str, int] | None = 'auto', config: dict | None = None, cachedir=None, cachename=None) DataArray[source]

Download bathymetry data from RWS.

Bathymetry is available at 20m resolution and at 1m resolution. The 20m resolution is available for large water bodies, but not in the major rivers. The 1m dataset covers the major waterbodies across all of the Netherlands.

Parameters:
  • extent (tuple) – extent of the model domain

  • resolution (str, optional) – resolution of the bathymetry data, “1m” or “20m”. The default is “20m”.

  • res (float, optional) – resolution of the output data array. The default is None, which uses resolution of the input datasets. Resampling method is provided by the method kwarg.

  • method (str, optional) – Rasterio resampling method. The default is None. Pre-defined method are “first”, “last”, “min”, “nearest”, “sum” or “count”. But custom callables are also supported. See the rasterio documentation for more information.

  • chunks (dict, optional) – chunks for the output data array. The default is “auto”, which lets xarray/dask pick the chunksize. Set to None to avoid chunking.

  • config (dict, optional) – configuration dictionary containing urls and layer numbers for GDR data. The default is None, which uses the default configuration provided by the function get_gdr_configuration().

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

bathymetry – bathymetry data

Return type:

xr.DataArray

nlmod.read.rws.download_bathymetry_gdf(resolution: Literal['20m', '1m'] = '20m', extent: list[float] | None = None, config: dict | None = None) GeoDataFrame[source]

Get bathymetry dataframe from RWS.

Note that the 20m resolution does not contain bathymetry data for the major rivers. If you need the bathymetry of the major rivers, use the 1m resolution.

Parameters:
  • resolution (str, optional) – resolution of the bathymetry data, “1m” or “20m”. The default is “20m”.

  • extent (tuple, optional) – extent of the model domain. The default is None.

  • config (dict, optional) – configuration dictionary containing urls and layer numbers for GDR data. The default is None, which uses the default configuration provided by the function get_gdr_configuration().

nlmod.read.rws.get_bathymetry(extent: list[float], resolution: Literal['20m', '1m'] = '20m', res: float | None = None, method: str | Callable | None = None, chunks: str | dict[str, int] | None = 'auto', config: dict | None = None, cachedir=None, cachename=None) DataArray[source]

Get bathymetry data from RWS.

Deprecated since version 0.10.0: get_bathymetry will be removed in nlmod 1.0.0, it is replaced by download_bathymetry because of new naming convention https://github.com/gwmod/nlmod/issues/47

Bathymetry is available at 20m resolution and at 1m resolution. The 20m resolution is available for large water bodies, but not in the major rivers. The 1m dataset covers the major waterbodies across all of the Netherlands.

Parameters:
  • extent (tuple) – extent of the model domain

  • resolution (str, optional) – resolution of the bathymetry data, “1m” or “20m”. The default is “20m”.

  • res (float, optional) – resolution of the output data array. The default is None, which uses resolution of the input datasets. Resampling method is provided by the method kwarg.

  • method (str, optional) – Rasterio resampling method. The default is None. Pre-defined method are “first”, “last”, “min”, “nearest”, “sum” or “count”. But custom callables are also supported. See the rasterio documentation for more information.

  • chunks (dict, optional) – chunks for the output data array. The default is “auto”, which lets xarray/dask pick the chunksize. Set to None to avoid chunking.

  • config (dict, optional) – configuration dictionary containing urls and layer numbers for GDR data. The default is None, which uses the default configuration provided by the function get_gdr_configuration().

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

bathymetry – bathymetry data

Return type:

xr.DataArray

nlmod.read.rws.get_bathymetry_gdf(resolution: Literal['20m', '1m'] = '20m', extent: list[float] | None = None, config: dict | None = None) GeoDataFrame[source]

Get bathymetry dataframe from RWS.

Deprecated since version 0.10.0: get_bathymetry_gdf will be removed in nlmod 1.0.0, it is replaced by download_bathymetry_gdf because of new naming convention https://github.com/gwmod/nlmod/issues/47

Note that the 20m resolution does not contain bathymetry data for the major rivers. If you need the bathymetry of the major rivers, use the 1m resolution.

Parameters:
  • resolution (str, optional) – resolution of the bathymetry data, “1m” or “20m”. The default is “20m”.

  • extent (tuple, optional) – extent of the model domain. The default is None.

  • config (dict, optional) – configuration dictionary containing urls and layer numbers for GDR data. The default is None, which uses the default configuration provided by the function get_gdr_configuration().

nlmod.read.rws.get_gdf_surface_water(ds=None, extent=None, cachedir=None, cachename=None)[source]

Read a shapefile with surface water as a geodataframe, cut by the extent of the model.

Parameters:
  • ds (xr.DataSet, None, optional) – dataset containing relevant model information

  • extent (list, tuple or np.array, optional) – desired model extent (xmin, xmax, ymin, ymax)

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

gdf_opp_water – surface water geodataframe.

Return type:

GeoDataframe

nlmod.read.rws.get_gdr_configuration() dict[source]

Get configuration for GDR data.

Note: Currently only includes configuration for bathymetry data. Other datasets can be added in the future. See https://geo.rijkswaterstaat.nl/arcgis/rest/services/GDR/ for available data.

Returns:

config – configuration dictionary containing urls and layer numbers for GDR data.

Return type:

dict

nlmod.read.rws.get_northsea(ds, gdf=None, da_name='northsea', cachedir=None, cachename=None)[source]

Get Dataset which is 1 at the northsea and 0 everywhere else.

Sea is defined by rws surface water shapefile.

Deprecated since version 0.10.0: get_northsea will be removed in nlmod 1.0.0, it is replaced by discretize_northsea because of new naming convention https://github.com/gwmod/nlmod/issues/47

Parameters:
  • ds (xr.DataSet, None, optional) – dataset containing relevant model information

  • gdf (gpd.GeoDataFrame or None, optional) – geometries of the surface water, if None the geometries are obtained using get_gdf_surface_water. The default is None.

  • da_name (str, optional) – name of the datavar that identifies sea cells

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

ds_out – Dataset with a single DataArray, this DataArray is 1 at sea and 0 everywhere else. Grid dimensions according to ds.

Return type:

xr.DataSet

nlmod.read.rws.get_surface_water(ds, gdf=None, da_basename='rws_oppwater', cachedir=None, cachename=None)[source]

Create 3 data-arrays from the shapefile with surface water:

Deprecated since version 0.10.0: get_surface_water will be removed in nlmod 1.0.0, it is replaced by discretize_surface_water because of new naming convention https://github.com/gwmod/nlmod/issues/47

  • area: area of the shape in the cell

  • cond: conductance based on the area and “bweerstand” column in shapefile

  • stage: surface water level based on the “peil” column in the shapefile

Parameters:
  • ds (xr.DataSet) – xarray with model data

  • gdf (gpd.GeoDataFrame or None, optional) – geometries of the surface water, if None the geometries are obtained using get_gdf_surface_water. The default is None.

  • da_basename (str) – name of the polygon shapes, name is used as a prefix to store data arrays in ds

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

ds – dataset with modelgrid data.

Return type:

xarray.Dataset

nlmod.read.waterboard._set_column_from_columns(gdf, set_column, from_columns, nan_values=None)[source]

Populate a target column from one or more source columns.

The function creates set_column and fills missing values by iterating over from_columns in order. Each source can be either a single column name or a list of column names. For a list input, the row-wise mean of available values is used. Existing values in set_column are preserved once filled.

Parameters:
  • gdf (pandas.DataFrame or geopandas.GeoDataFrame) – Input table that contains the source columns.

  • set_column (str) – Name of the column that will be created and filled.

  • from_columns (str or list, optional) –

    Source definition used to populate set_column. Supported entries are: - str: copy values from this column. - list[str]: use row-wise mean of these columns. - list containing a mix of the above: processed in order until all

    possible values are filled.

  • nan_values (float, int or list, optional) – Value or values that should be interpreted as missing and replaced by np.nan in the resulting set_column.

Returns:

gdf with the newly created set_column.

Return type:

pandas.DataFrame or geopandas.GeoDataFrame

Raises:

ValueError – If set_column already exists in gdf.

nlmod.read.waterboard.download_data(wb, data_kind, extent=None, max_record_count=None, config=None, cachedir=None, cachename=None, **kwargs)[source]

Get the data for a Waterboard and a specific data_kind.

Parameters:
  • ws (str) – The name of the waterboard.

  • data_kind (str) – The kind of data you like to download. Possible values are ‘watercourses’, ‘level_areas’ and ‘level_deviations’

  • extent (tuple or list of length 4, optional) – THe extent of the data you like to donload: (xmin, xmax, ymin, ymax). Download everything when extent is None. The default is None.

  • max_record_count (int, optional) – THe maximum number of records that are downloaded in each call to the webservice. When max_record_count is None, the maximum is set equal to the maximum of the server. The default is None.

  • config (dict, optional) – A dictionary with properties of the data sources of the Waterboards. When None, the configuration is retreived from the method get_configuration(). The default is None.

  • **kwargs (dict) – Optional arguments which are passed onto arcrest() or wfs().

Raises:

DESCRIPTION.

cachedirstr or None, optional

directory to save cache. If None no cache is used. Default is None.

cachenamestr or None, optional

filename of netcdf cache. If None no cache is used. Default is None.

Returns:

gdf – A GeoDataFrame containing data from the waterboard (polygons for level_areas/level_deviations and lines for watercourses).

Return type:

GeoDataFrame

nlmod.read.waterboard.download_polygons(**kwargs)[source]

Get the location of the Dutch Waterboards as a Polygon GeoDataFrame.

nlmod.read.waterboard.get_configuration()[source]

Get the configuration of of the data sources of the Waterboards.

nlmod.read.waterboard.get_data(wb, data_kind, extent=None, max_record_count=None, config=None, cachedir=None, cachename=None, **kwargs)[source]

Get the data for a Waterboard and a specific data_kind.

Deprecated since version 0.10.0: get_data will be removed in nlmod 1.0.0, it is replaced by download_data because of new naming convention https://github.com/gwmod/nlmod/issues/47

Parameters:
  • ws (str) – The name of the waterboard.

  • data_kind (str) – The kind of data you like to download. Possible values are ‘watercourses’, ‘level_areas’ and ‘level_deviations’

  • extent (tuple or list of length 4, optional) – THe extent of the data you like to donload: (xmin, xmax, ymin, ymax). Download everything when extent is None. The default is None.

  • max_record_count (int, optional) – THe maximum number of records that are downloaded in each call to the webservice. When max_record_count is None, the maximum is set equal to the maximum of the server. The default is None.

  • config (dict, optional) – A dictionary with properties of the data sources of the Waterboards. When None, the configuration is retreived from the method get_configuration(). The default is None.

  • **kwargs (dict) – Optional arguments which are passed onto arcrest() or wfs().

Raises:

DESCRIPTION.

cachedirstr or None, optional

directory to save cache. If None no cache is used. Default is None.

cachenamestr or None, optional

filename of netcdf cache. If None no cache is used. Default is None.

Returns:

gdf – A GeoDataFrame containing data from the waterboard (polygons for level_areas/level_deviations and lines for watercourses).

Return type:

GeoDataFrame

nlmod.read.waterboard.get_polygons(**kwargs)[source]

Get the location of the Dutch Waterboards as a Polygon GeoDataFrame.

Deprecated since version 0.10.0: get_polygons will be removed in nlmod 1.0.0, it is replaced by download_polygons because of new naming convention https://github.com/gwmod/nlmod/issues/47

dims

nlmod.dims.base._get_structured_grid_ds(xedges, yedges, nlay=1, top=nan, botm=nan, xorigin=0.0, yorigin=0.0, angrot=0, attrs=None, crs=None)[source]

Create an xarray dataset with structured grid geometry.

Parameters:
  • xedges (array_like) – A 1D array of the x coordinates of the grid edges.

  • yedges (array_like) – A 1D array of the y coordinates of the grid edges.

  • nlay (int, optional) – The number of layers in the grid. Default is 1.

  • top (array_like, optional) – A 2D array of the top elevation of the grid cells. Default is NaN.

  • botm (array_like, optional) – A 3D array of the bottom elevation of the grid cells. Default is NaN.

  • xorigin (int or float, optional) – lower left x coordinate of the model grid. When angrot == 0, xorigin is added to the first two values of extent. Otherwise it is the x-coordinate of the point the grid is rotated around, and xorigin is added to the Dataset-attributes. The default is 0.0.

  • yorigin (int or float, optional) – lower left y coordinate of the model grid. When angrot == 0, yorigin is added to the last two values of extent. Otherwise it is the y-coordinate of the point the grid is rotated around, and yorigin is added to the Dataset-attributes. The default is 0.0.

  • angrot (int or float, optinal) – the rotation of the grid in counter clockwise degrees. When angrot != 0 the grid is rotated, and all coordinates of the model are in model coordinates. See https://nlmod.readthedocs.io/en/stable/examples/11_grid_rotation.html for more infomation. The default is 0.0.

  • attrs (dict, optional) – A dictionary of attributes to add to the xarray dataset. Default is an empty dictionary.

  • crs (dict or str, optional) – A dictionary or string describing the coordinate reference system of the grid. Default is None.

Returns:

ds – An xarray dataset with the following data variables and coordinates:

  • top : a 2D array of the top elevation of the grid cells

  • botm : a 3D array of the bottom elevation of the grid cells

  • x : a 1D array of the x coordinates of the grid cell centers

  • y : a 1D array of the y coordinates of the grid cell centers

  • layer : a 1D array of the layer indices

  • xc : a 2D array of the x coordinates of the grid cell centers, after rotation if angrot is not 0.0 (optional)

  • yc : a 2D array of the y coordinates of the grid cell centers, after rotation if angrot is not 0.0 (optional)

The dataset also includes the attributes specified in the attrs dictionary, and a coordinate reference system specified by crs, if provided.

Return type:

xarray.Dataset

nlmod.dims.base._get_vertex_grid_ds(x, y, xv, yv, cell2d, extent, nlay=1, top=nan, botm=nan, xorigin=0.0, yorigin=0.0, angrot=0.0, attrs=None, crs=None)[source]

Create an xarray dataset with vertex-based grid geometry.

Parameters:
  • x (array_like) – A 1D array of the x coordinates of the grid cell centers.

  • y (array_like) – A 1D array of the y coordinates of the grid cell centers.

  • xv (array_like) – A 1D array of the x coordinates of the grid vertices.

  • yv (array_like) – A 1D array of the y coordinates of the grid vertices.

  • cell2d (array-like) – array-like with vertex grid cell2d info.

  • extent (list) – A list of [xmin, xmax, ymin, ymax] defining the extent of the model grid.

  • nlay (int or sequence of ints, optional) – The number of layers in the grid, or a sequence of layer indices. Default is 1.

  • top (array_like, optional) – A 2D array of the top elevation of the grid cells. Default is NaN.

  • botm (array_like, optional) – A 3D array of the bottom elevation of the grid cells. Default is NaN.

  • xorigin (int or float, optional) – lower left x coordinate of the model grid. When angrot == 0, xorigin is added to the first two values of extent. Otherwise it is the x-coordinate of the point the grid is rotated around, and xorigin is added to the Dataset-attributes. The default is 0.0.

  • yorigin (int or float, optional) – lower left y coordinate of the model grid. When angrot == 0, yorigin is added to the last two values of extent. Otherwise it is the y-coordinate of the point the grid is rotated around, and yorigin is added to the Dataset-attributes. The default is 0.0.

  • angrot (int or float, optinal) – the rotation of the grid in counter clockwise degrees. When angrot != 0 the grid is rotated, and all coordinates of the model are in model coordinates. See https://nlmod.readthedocs.io/en/stable/examples/11_grid_rotation.html for more infomation. The default is 0.0.

  • attrs (dict, optional) – A dictionary of attributes to add to the xarray dataset. Default is an empty dictionary.

  • crs (dict or str, optional) – A dictionary or string describing the coordinate reference system of the grid. Default is None.

Returns:

ds – An xarray dataset with the following data variables and coordinates:

  • top : a 2D array of the top elevation of the grid cells

  • botm : a 3D array of the bottom elevation of the grid cells

  • x : a 1D array of the x coordinates of the grid cell centers

  • y : a 1D array of the y coordinates of the grid cell centers

  • layer : a 1D array of the layer indices

  • xv : a 1D array of the x coordinates of the grid vertices

  • yv : a 1D array of the y coordinates of the grid vertices

The dataset also includes the attributes specified in the attrs dictionary, and a coordinate reference system specified by crs, if provided.

Return type:

xarray.Dataset

nlmod.dims.base.extrapolate_ds(ds, mask=None, layer='layer', mask_values=None)[source]

Fill missing data in layermodel based on nearest interpolation.

Used for ensuring layer model contains data everywhere. Useful for filling in data beneath the sea for coastal groundwater models, or models near the border of the Netherlands.

Parameters:
  • ds (xarray.DataSet) – Model layer DataSet

  • mask (np.ndarray, optional) – Boolean mask for each cell, with a value of True if its value needs to be determined. When mask is None, it is determined from the botm- variable. The default is None.

  • layer (str, optional) – The name of the layer dimension. The default is ‘layer’.

  • mask_values (np.ndarray, optional) – Boolean mask for each cell, with a value of True if its value is used to fill data in mask. When mask_values is None, it is determined from mask. The default is None.

Returns:

ds – filled layermodel

Return type:

xarray.DataSet

nlmod.dims.base.get_ds(extent, delr=100.0, delc=None, model_name=None, model_ws=None, layer=None, top=0.0, botm=None, kh=10.0, kv=1.0, crs='+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725 +units=m +no_defs', xorigin=0.0, yorigin=0.0, angrot=0.0, attrs=None, extrapolate=True, fill_nan=True, transport=False, version_tag=None, download_exe=True, **kwargs)[source]

Create a model dataset from scratch.

Parameters:
  • extent (list, tuple or np.array) – desired model extent (xmin, xmax, ymin, ymax)

  • delr (int, float, list, tuple or array, optional) – The gridsize along columns (dx). The default is 100. meter.

  • delc (None, int, float, list, tuple or array, optional) – The gridsize along rows (dy). Set to delr when None. If None delc=delr. The default is None.

  • model_name (str, optional) – name of the model. The default is None

  • model_ws (str, optional) – workspace of the model. This is where modeldata is saved to. The default is None.

  • layer (int, list, tuple or ndarray, optional) – The names or index of the layers of the model. When layer is an integer it is the number of layers. When layer is None, the number of layers is caluclated from botm. When botm is None as well, the number of layers is set to 10. The default is None.

  • top (float, list or ndarray, optional) – The top of the model. It has to be of shape (len(y), len(x)) or it is transformed into that shape if top is a float. The default is 0.0.

  • botm (list or ndarray, optional) – The botm of the model layers. It has to be of shape (len(layer), len(y), len(x)) or it is transformed to that shape if botm is or a list/array of len(layer). When botm is None, a botm is generated with a constant layer thickness of 10 meter. The default is None.

  • kh (float, list or ndarray, optional) – The horizontal conductivity of the model layers. It has to be of shape (len(layer), len(y), len(x)) or it is transformed to that shape if kh is a float or a list/array of len(layer). The default is 10.0.

  • kv (float, list or ndarray, optional) – The vertical conductivity of the model layers. It has to be of shape (len(layer), len(y), len(x)) or it is transformed to that shape if kv is a float or a list/array of len(layer). The default is 1.0.

  • crs (int, optional) – The coordinate reference system of the model. The default is 28992.

  • xorigin (int or float, optional) – lower left x coordinate of the model grid. When angrot == 0, xorigin is added to the first two values of extent. Otherwise it is the x-coordinate of the point the grid is rotated around, and xorigin is added to the Dataset-attributes. The default is 0.0.

  • yorigin (int or float, optional) – lower left y coordinate of the model grid. When angrot == 0, yorigin is added to the last two values of extent. Otherwise it is the y-coordinate of the point the grid is rotated around, and yorigin is added to the Dataset-attributes. The default is 0.0.

  • angrot (int or float, optinal) – the rotation of the grid in counter clockwise degrees. When angrot != 0 the grid is rotated, and all coordinates of the model are in model coordinates. See https://nlmod.readthedocs.io/en/stable/examples/11_grid_rotation.html for more infomation. The default is 0.0.

  • attrs (dict, optional) – Attributes of the model dataset. The default is None.

  • extrapolate (bool, optional) – When true, extrapolate data-variables, into the sea or other areas with only nans. The default is True

  • fill_nan (bool, optional) – if True nan values in the top, botm, kh and kv are filled using the fill_nan_top_botm_kh_kv function. Layers with only nan values in the botm are removed.

  • transport (bool, optional) – flag indicating whether dataset includes data for a groundwater transport model (GWT). Default is False, no transport.

  • version_tag (str, default None) – GitHub release ID: for example “18.0” or “latest”. If version_tag is provided, the most recent installation location of MODFLOW is found in flopy metadata that respects version_tag. If not found, the executables are downloaded if download_exe is True (see below). Not compatible with exe_name.

  • download_exe (bool, optional) – If True, download the executable if it is not found locally. The default is True.

  • **kwargs (dict) – Kwargs are passed into nlmod.to_model_ds(). See nlmod.to_model_ds() for more information.

Returns:

The model dataset.

Return type:

xr.Dataset

nlmod.dims.base.set_ds_attrs(ds, model_name, model_ws, mfversion='mf6', exe_name=None, version_tag=None, download_exe=True)[source]

Set the attribute of a model dataset.

Parameters:
  • ds (xarray dataset) – An existing model dataset

  • model_name (str) – name of the model.

  • model_ws (str or None) – workspace of the model. This is where modeldata is saved to.

  • mfversion (str, optional) – modflow version. The default is “mf6”.

  • exe_name (str, optional) – path to modflow executable, default is None, which assumes binaries are available in nlmod/bin directory. Binaries can be downloaded using nlmod.util.download_mfbinaries().

  • version_tag (str, default None) – GitHub release ID: for example “18.0” or “latest”. If version_tag is provided, the most recent installation location of MODFLOW is found in flopy metadata that respects version_tag. If not found, the executables are downloaded. Not compatible with exe_name.

  • download_exe (bool, optional) – If True, download the executable if it is not found locally. The default is True.

Returns:

ds – model dataset.

Return type:

xarray dataset

nlmod.dims.base.to_model_ds(ds, model_name=None, model_ws=None, extent=None, delr=100.0, delc=None, fill_nan=True, extrapolate=True, anisotropy=10, fill_value_kh=1.0, fill_value_kv=0.1, xorigin=0.0, yorigin=0.0, angrot=0.0, drop_attributes=True, transport=False, remove_nan_layers=True, version_tag=None, download_exe=True)[source]

Transform an input dataset to a groundwater model dataset.

Optionally select a different grid size.

Parameters:
  • ds (xarray.dataset) – A layer model dataset.

  • model_name (str, optional) – name of the model. The default is None

  • model_ws (str, optional) – workspace of the model. This is where modeldata is saved to. The default is None

  • extent (list or tuple of length 4, optional) – The extent of the new grid. Get from ds when None. The default is None.

  • delr (int, float, list, tuple or array, optional) – The gridsize along columns (dx). The default is 100. meter.

  • delc (None, int, float, list, tuple or array, optional) – The gridsize along rows (dy). Set to delr when None. If None delc=delr The default is None.

  • fill_nan (bool, optional) – if True nan values in the top, botm, kh and kv are filled using the fill_nan_top_botm_kh_kv function. Layers with only nan values in the botm are removed.

  • extrapolate (bool, optional) – When true, extrapolate data-variables, into the sea or other areas with only nans. The default is True

  • anisotropy (int or float) – factor to calculate kv from kh or the other way around

  • fill_value_kh (int or float, optional) – use this value for kh if there is no data. The default is 1.0.

  • fill_value_kv (int or float, optional) – use this value for kv if there is no data. The default is 0.1.

  • xorigin (int or float, optional) – lower left x coordinate of the model grid. When angrot == 0, xorigin is added to the first two values of extent. Otherwise it is the x-coordinate of the point the grid is rotated around, and xorigin is added to the Dataset-attributes. The default is 0.0.

  • yorigin (int or float, optional) – lower left y coordinate of the model grid. When angrot == 0, yorigin is added to the last two values of extent. Otherwise it is the y-coordinate of the point the grid is rotated around, and yorigin is added to the Dataset-attributes. The default is 0.0.

  • angrot (int or float, optinal) – the rotation of the grid in counter clockwise degrees. When angrot != 0 the grid is rotated, and all coordinates of the model are in model coordinates. See https://nlmod.readthedocs.io/en/stable/examples/11_grid_rotation.html for more infomation. The default is 0.0.

  • drop_attributes (bool, optional) – if True drop the attributes from the layer model dataset. Otherwise keep the attributes. Default is True.

  • transport (bool, optional) – flag indicating whether dataset includes data for a groundwater transport model (GWT). Default is False, no transport.

  • remove_nan_layers (bool, optional) – if True remove layers with only nan values in the botm. Default is True.

  • version_tag (str, default None) – GitHub release ID: for example “18.0” or “latest”. If version_tag is provided, the most recent installation location of MODFLOW is found in flopy metadata that respects version_tag. If not found, the executables are downloaded. Not compatible with exe_name.

  • download_exe (bool, optional) – If True, download the executable if it is not found locally. The default is True.

Returns:

ds – the model Dataset.

Return type:

xarray.dataset

Module containing model grid functions.

  • project data on different grid types

  • obtain various types of reclists from a grid that can be used as input for a MODFLOW package

  • fill, interpolate and resample grid data

nlmod.dims.grid._agg_area_weighted(gdf, col)[source]
nlmod.dims.grid._agg_length_weighted(gdf, col)[source]
nlmod.dims.grid._agg_max_area(gdf, col)[source]
nlmod.dims.grid._agg_max_length(gdf, col)[source]
nlmod.dims.grid._agg_nearest(gdf, col, modelgrid)[source]
nlmod.dims.grid._get_aggregates_values(group, fields_methods, modelgrid=None)[source]
nlmod.dims.grid._get_attrs(ds)[source]
nlmod.dims.grid.affine_transform_gdf(gdf, affine)[source]

Apply an affine transformation to a geopandas GeoDataFrame.

nlmod.dims.grid.aggregate_vector_per_cell(gdf, fields_methods, modelgrid=None)[source]

Aggregate vector features per cell.

Parameters:
  • gdf (geopandas.GeoDataFrame) – GeoDataFrame containing points, lines or polygons per grid cell.

  • fields_methods (dict) – fields (keys) in the Geodataframe with their aggregation method (items) aggregation methods can be: max, min, mean, sum, length_weighted (lines), max_length (lines), area_weighted (polygon), max_area (polygon), and functions supported by pandas.*GroupBy.agg().

  • modelgrid (flopy Groundwater flow modelgrid) – only necesary if one of the field methods is ‘nearest’

Returns:

celldata – DataFrame with aggregated surface water parameters per grid cell

Return type:

pd.DataFrame

nlmod.dims.grid.col_to_list(col_in, ds, cellids)[source]

Convert array data in ds to a list of values for specific cells.

This function is typically used to create a rec_array with stress period data for the modflow packages. Can be used for structured and vertex grids.

Parameters:
  • col_in (xr.DataArray, str, int or float) – if col_in is a str type it is the name of the variable in ds (if it exists). if col_in is an int or a float it is a value that will be used for all cells in cellids.

  • ds (xr.Dataset) – dataset with model data. Can have dimension (layer, y, x) or (layer, icell2d).

  • cellids (tuple of numpy arrays) –

    tuple with indices of the cells that will be used to create the list with values. There are 3 options:

    1. cellids contains (layers, rows, columns)

    2. cellids contains (rows, columns) or (layers, icell2ds)

    3. cellids contains (icell2ds)

Raises:

ValueError – raised if the cellids are in the wrong format.

Returns:

col_lst – raster values from ds presented in a list per cell.

Return type:

list

nlmod.dims.grid.cols_to_reclist(ds, cellids, *args, cellid_column=0)[source]

Create a reclist for stress period data from a set of cellids.

Parameters:
  • ds (xr.Dataset) – dataset with model data. Should have dimensions (layer, icell2d).

  • cellids (tuple of length 2 or 3) – tuple with indices of the cells that will be used to create the list. For a structured grid, cellids represents (layer, row, column). For a vertex grid cellid reprsents (layer, icell2d).

  • args (xr.DataArray, str, int or float) – the args parameter represents the data to be used as the columns in the reclist. See col_to_list of the allowed values.

  • cellid_column (int, optional) – Adds the cellid ((layer, row, col) or (layer, icell2d)) to the reclist in this column number. Do not add cellid when cellid_column is None. The default is 0.

nlmod.dims.grid.da_to_reclist(ds, mask, col1=None, col2=None, col3=None, layer=0, aux=None, first_active_layer=False, only_active_cells=True)[source]

Create a reclist for stress period data from a model dataset.

Used for vertex grids.

Parameters:
  • ds (xr.Dataset) – dataset with model data and dimensions (layer, icell2d)

  • mask (xr.DataArray for booleans) – True for the cells that will be used in the rec list.

  • col1 (str, int or float, optional) –

    1st column of the reclist, if None the reclist will be a list with (cellid,) for each row.

    col1 should be the following value for each package (can also be the

    name of a timeseries): rch: recharge [L/T] ghb: head [L] drn: drain level [L] chd: head [L]

  • col2 (str, int or float, optional) –

    2nd column of the reclist, if None the reclist will be a list with (cellid, col1) for each row.

    col2 should be the following value for each package (can also be the

    name of a timeseries): ghb: conductance [L^2/T] drn: conductance [L^2/T]

  • col3 (str, int or float, optional) –

    3th column of the reclist, if None the reclist will be a list with (cellid, col1, col2) for each row.

    col3 should be the following value for each package (can also be the

    name of a timeseries): riv: bottom [L]

  • aux (str or list of str, optional) – list of auxiliary variables to include in reclist

  • layer (int, optional) – layer used in the reclist. Not used if layer is in the dimensions of mask or if first_active_layer is True. The default is 0

  • first_active_layer (bool, optional) – If True an extra mask is applied to use the first active layer of each cell in the grid. Not used if layer is in the dimensions of mask. The default is False.

  • only_active_cells (bool, optional) – If True an extra mask is used to only include cells with an idomain of 1. The default is True.

Returns:

reclist – every row consists of ((layer, icell2d), col1, col2, col3).

Return type:

list of tuples

nlmod.dims.grid.ds_to_gridprops(ds_in, gridprops, method='nearest', icvert_nodata=-1)[source]

Resample a xarray dataset of a structured grid to a new dataset with a vertex grid.

Returns a dataset with resampled variables and the untouched variables.

Parameters:
  • ds_in (xr.Dataset) – dataset with dimensions (layer, y, x). y and x are from the original structured grid

  • gridprops (dictionary) – dictionary with grid properties output from gridgen. Used as the definition of the vertex grid.

  • method (str, optional) – type of interpolation used to resample. The default is ‘nearest’.

  • icvert_nodata (int, optional) – integer to represent nodata-values in cell2d array. Defaults to -1.

Returns:

ds_out – dataset with resampled variables and the untouched variables.

Return type:

xr.Dataset

nlmod.dims.grid.gdf_area_to_da(gdf, ds, ix=None, index_name=None, desc='Intersecting with grid', silent=False, sparse=True, **kwargs)[source]

Calculate the area of overlap between polygons in a GeoDataFrame and a model grid.

This function computes the area of each geometry in gdf that falls within each model grid cell defined in the xarray ds Dataset. The result is returned as an xarray DataArray, with dimensions based on the model grid and the index of the input GeoDataFrame.

Parameters:
  • gdf (geopandas.GeoDataFrame) – GeoDataFrame containing polygon geometries to be intersected with the model grid.

  • ds (xr.Dataset) – The xarray dataset that defines the grid. When the grid is rotated, the geodataframe is transformed in model coordinates.

  • ix (flopy.utils.GridIntersect, optional) – GridIntersect, if not provided the modelgrid is determined from ds. The default is None.

  • index_name (str, optional) – Name to use for the index dimension of the resulting DataArray. If None, uses the name of gdf.index, or defaults to ‘index’.

  • desc (string, optional) – The description of the progressbar. The default is ‘Intersecting with grid’.

  • silent (bool, optional) – Do not show a progressbar when silent is True. The default is False.

  • sparse (bool, default True) – If True, attempts to return a sparse array. Falls back to dense if sparse is not available.

  • **kwargs (dict) – Additional keyword arguments passed to ix.intersect.

Returns:

A DataArray with shape (y, x, index) for structured grids or (icell2d, index) for unstructured grids. Each value represents the area of a given geometry within a specific grid cell.

Return type:

xarray.DataArray

Notes

  • Uses flopy.utils.GridIntersect under the hood for spatial intersection.

  • For rotated grids, geometries are transformed into model space using an affine transformation.

  • If the sparse package is installed and sparse=True, a sparse.COO array is returned.

  • Suitable for use in spatial weighting or disaggregation tasks.

Raises:

ImportWarning – If sparse=True but the sparse library is not installed, a warning is issued and a dense array is used instead.

nlmod.dims.grid.gdf_to_bool_da(gdf, ds, ix=None, buffer=0.0, contains_centroid=False, min_area_fraction=None, **kwargs)[source]

Return True if grid cell is partly in polygons, False otherwise.

This function returns True for grid cells that are partly in the polygons. If contains_centroid is True, only cells are returned where the centroid is in the polygon. If min_area_fraction is set, only cells are returned where the area of the intersection is larger than min_area_fraction * cell_area.

ix can be provided to speed up the process. If not provided it is computed from ds.

Parameters:
  • gdf (geopandas.GeoDataFrame or shapely.geometry) – shapes that will be rasterised.

  • ds (xr.Dataset) – Dataset with model data.

  • ix (flopy.utils.GridIntersect, optional) – If not provided it is computed from ds. Speeds up the the function

  • buffer (float, optional) – The distance to buffer around the geometries in meters. A positive distance produces a dilation, a negative distance an erosion. The default is 0.

  • contains_centroid (bool, optional) – if True, only store intersection result if cell centroid is contained within intersection shape, only used if shape type is “polygon”

  • min_area_fraction (float, optional) – float defining minimum intersection area threshold, if intersection area is smaller than min_frac_area * cell_area, do not store intersection result.

  • **kwargs (keyword arguments) – keyword arguments are passed to the intersect-method.

Returns:

da – True if polygon is in cell, False otherwise.

Return type:

xr.DataArray

nlmod.dims.grid.gdf_to_bool_ds(gdf, ds, da_name, keep_coords=None, ix=None, buffer=0.0, contains_centroid=False, min_area_fraction=None, **kwargs)[source]

Return True if grid cell is partly in polygons, False otherwise.

This function returns True for grid cells that are partly in the polygons. If contains_centroid is True, only cells are returned where the centroid is in the polygon. If min_area_fraction is set, only cells are returned where the area of the intersection is larger than min_area_fraction * cell_area.

ix can be provided to speed up the process. If not provided it is computed from ds.

Parameters:
  • gdf (geopandas.GeoDataFrame or shapely.geometry) – shapes that will be rasterised.

  • ds (xr.Dataset) – Dataset with model data.

  • da_name (str) – The name of the variable with boolean data in the ds_out

  • keep_coords (tuple or None, optional) – the coordinates in ds the you want keep in your empty ds. If None all coordinates are kept from original ds. The default is None.

  • ix (flopy.utils.GridIntersect, optional) – If not provided it is computed from ds. Speeds up the the function

  • buffer (float, optional) – The distance to buffer around the geometries in meters. A positive distance produces a dilation, a negative distance an erosion. The default is 0.

  • contains_centroid (bool, optional) – if True, only store intersection result if cell centroid is contained within intersection shape, only used if shape type is “polygon”

  • min_area_fraction (float, optional) – float defining minimum intersection area threshold, if intersection area is smaller than min_frac_area * cell_area, do not store intersection result.

  • **kwargs (keyword arguments) – keyword arguments are passed to the intersect_*-methods.

Returns:

ds_out – True if polygon is in cell, False otherwise.

Return type:

xr.Dataset

nlmod.dims.grid.gdf_to_count_da(gdf, ds, ix=None, buffer=0.0, **kwargs)[source]

Counts in how many polygons a coordinate of ds appears.

Parameters:
  • gdf (geopandas.GeoDataFrame or shapely.geometry) – shapes that will be rasterised.

  • ds (xr.Dataset) – Dataset with model data.

  • ix (flopy.utils.GridIntersect, optional) – If not provided it is computed from ds.

  • buffer (float, optional) – buffer around the geometries. The default is 0.

  • **kwargs (keyword arguments) – keyword arguments are passed to the intersect_*-methods.

Returns:

da – 1 if polygon is in cell, 0 otherwise. Grid dimensions according to ds.

Return type:

xr.DataArray

nlmod.dims.grid.gdf_to_count_ds(gdf, ds, da_name, keep_coords=None, ix=None, buffer=0.0, **kwargs)[source]

Counts in how many polygons a coordinate of ds appears.

Parameters:
  • gdf (geopandas.GeoDataFrame) – polygon shapes with surface water.

  • ds (xr.Dataset) – Dataset with model data.

  • da_name (str) – The name of the variable with boolean data in the ds_out

  • keep_coords (tuple or None, optional) – the coordinates in ds the you want keep in your empty ds. If None all coordinates are kept from original ds. The default is None.

  • ix (flopy.utils.GridIntersect, optional) – If not provided it is computed from ds.

  • buffer (float, optional) – buffer around the geometries. The default is 0.

  • **kwargs (keyword arguments) – keyword arguments are passed to the intersect_*-methods.

Returns:

ds_out – Dataset with a single DataArray, this DataArray is 1 if polygon is in cell, 0 otherwise. Grid dimensions according to ds and mfgrid.

Return type:

xr.Dataset

nlmod.dims.grid.gdf_to_da(gdf, ds, column, agg_method=None, fill_value=nan, min_total_overlap=0.0, ix=None)[source]

Project vector data on a grid. Aggregate data if multiple geometries are in a single cell. Supports structured and vertex grids. This method replaces gdf_to_data_array_struc.

Parameters:
  • gdf (geopandas.GeoDataframe) – vector data can only contain a single geometry type.

  • ds (xr.Dataset) – model Datset

  • column (str) – column name in the geodataframe.

  • agg_method (str, optional) – aggregation method to handle multiple geometries in one cell, options are: - max, min, mean, sum - length_weighted (lines), max_length (lines), - area_weighted (polygon), max_area (polygon). The default is ‘max’.

  • fill_value (float or int, optional) – The value to fill in da outside gdf. The default is np.NaN

  • min_total_overlap (float, optional) – Only assign cells with a gdf-area larger than min_total_overlap * cell-area. The default is 0.0

  • ix (GridIntersect, optional) – If not provided it is computed from ds.

Returns:

da – The DataArray with the projected vector data.

Return type:

xr DataArray

nlmod.dims.grid.gdf_to_data_array_struc(gdf, gwf, field='VALUE', agg_method=None, interp_method=None)[source]

Project vector data on a structured grid. Aggregate data if multiple geometries are in a single cell.

Parameters:
  • gdf (geopandas.GeoDataframe) – vector data can only contain a single geometry type.

  • gwf (flopy groundwater flow model) – model with a structured grid.

  • field (str, optional) – column name in the geodataframe. The default is ‘VALUE’.

  • interp_method (str or None, optional) – method to obtain values in cells without geometry by interpolating between cells with values. Options are ‘nearest’ and ‘linear’.

  • agg_method (str, optional) – aggregation method to handle multiple geometries in one cell, options are: - max, min, mean, - length_weighted (lines), max_length (lines), - area_weighted (polygon), max_area (polygon). The default is ‘max’.

Returns:

da – The DataArray with the projected vector data.

Return type:

xr DataArray

nlmod.dims.grid.gdf_to_grid(gdf, ml=None, ix=None, desc='Intersecting with grid', silent=False, cachedir=None, cachename=None, **kwargs)[source]

Intersect a geodataframe with the grid of a MODFLOW model.

Note: This method is a wrapper around the GridIntersect method in flopy.

Parameters:
  • gdf (geopandas.GeoDataFrame) – A GeoDataFrame that needs to be cut by the grid. The GeoDataFrame can consist of multiple types (Point, LineString, Polygon and the Multi-variants).

  • ml (flopy.modflow.Modflow or flopy.mf6.ModflowGwf or xr.Dataset, optional) – The flopy model or xarray dataset that defines the grid. When a Dataset is supplied, and the grid is rotated, the geodataframe is transformed in model coordinates. The default is None.

  • ix (flopy.utils.GridIntersect, optional) – GridIntersect, if not provided the modelgrid in ml is used.

  • desc (string, optional) – The description of the progressbar. The default is ‘Intersecting with grid’.

  • silent (bool, optional) – Do not show a progressbar when silent is True. The default is False.

  • **kwargs (keyword arguments) – keyword arguments are passed to the intersect_*-methods.

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

gdfg – The GeoDataFrame with the geometries per grid-cell.

Return type:

geopandas.GeoDataFrame

nlmod.dims.grid.get_affine(ds, sx=None, sy=None)[source]

Get the affine-transformation, from pixel to real-world coordinates.

nlmod.dims.grid.get_affine_mod_to_world(ds)[source]

Get the affine-transformation from model to real-world coordinates.

nlmod.dims.grid.get_affine_world_to_mod(ds)[source]

Get the affine-transformation from real-world to model coordinates.

nlmod.dims.grid.get_cell2d_from_ds(ds)[source]

Get the cell2d-list from a model dataset.

Flopy needs this list to build a disv-package

nlmod.dims.grid.get_cellids_from_xy(x, y, ds=None, modelgrid=None, gi=None)[source]

Map points to grid cells.

Note this is faster and more convenient than GridIntersect, because it provides the corresponding cellid for each point (instead of aggregating points per cell).

Parameters:
  • x (np.ndarray) – x-coordinates of points.

  • y (np.ndarray) – y-coordinates of points.

  • ds (xr.Dataset, optional) – Model dataset containing grid information. Must provide one of ds or modelgrid.

  • modelgrid (StructuredGrid or VertexGrid, optional) – Model grid object. Must provide one of ds or modelgrid.

  • gi (GridIntersect, optional) – GridIntersect instance, when passed saves some time building grid polygons.

Returns:

cellids – series with mapping between points and cellid of grid cell in which points are located.

Return type:

pd.Series

nlmod.dims.grid.get_dims_coords_from_modelgrid(mg)[source]

Get dimensions and coordinates from modelgrid.

Used to build new xarray DataArrays with appropriate dimensions and coordinates.

Parameters:

mg (flopy.discretization.Grid) – flopy modelgrid object

Returns:

  • dims (tuple of str) – tuple containing dimensions

  • coords (dict) – dictionary containing spatial coordinates derived from modelgrid

Raises:

ValueError – for unsupported grid types

nlmod.dims.grid.get_extent(ds, rotated=True, xmargin=0.0, ymargin=0.0)[source]

Get the model extent, corrected for angrot if necessary.

Parameters:
  • ds (xr.Dataset) – model dataset.

  • rotated (bool, optional) – if True, the extent is corrected for angrot. The default is True.

  • xmargin (float, optional) – margin to add to the x-extent. The default is 0.0.

  • ymargin (float, optional) – margin to add to the y-extent. The default is 0.0.

Returns:

extent – [xmin, xmax, ymin, ymax]

Return type:

list

nlmod.dims.grid.get_extent_gdf(ds, rotated=True, crs='EPSG:28992')[source]

Get the model extent as a GeoDataFrame with a polygon.

Parameters:
  • ds (xr.Dataset) – model dataset.

  • rotated (bool, optional) – if True, the extent is corrected for angrot. The default is True.

  • crs (str, optional) – Coordinate reference system. The default is “EPSG:28992”.

Returns:

gdf – GeoDataFrame containing the model extent as a polygon geometry.

Return type:

geopandas.GeoDataFrame

nlmod.dims.grid.get_extent_polygon(ds, rotated=True)[source]

Get the model extent, as a shapely Polygon.

nlmod.dims.grid.get_icell2d_from_xy(x, y, ds, gi=None, rotated=True)[source]

Get the icell2d value of a point defined by its x and y coordinates.

Parameters:
  • x (float) – The x-coordinate of the point.

  • y (float) – The y-coordinate of the point.

  • ds (xr.Dataset) – The model dataset.

  • gi (flopy.utils.GridIntersect, optional) – Can be supplied to speed up the calculation, as the generation of the GridIntersect-instance can take some time. The default is None.

  • rotated (bool, optional) – If the model grid has a rotation, and rotated is False, x and y are in model coordinates. Otherwise x and y are in real world coordinates. The default is True.

Raises:

ValueError – Raises a ValueError if the point is outside of the model grid.

Returns:

icell2d – The icell2d-number of the model cell containing the point (zero-based)

Return type:

int

nlmod.dims.grid.get_node_structured(lay, row, col, shape)[source]

Get the node number of a structured grid.

Parameters:
  • lay (int) – layer number.

  • row (int) – row number.

  • col (int) – column number.

  • shape (tuple) – shape of the model grid (nlay, nrow, ncol).

Returns:

node – node number.

Return type:

int

nlmod.dims.grid.get_node_vertex(lay, icell2d, shape)[source]

Get the node number of a vertex grid.

Parameters:
  • lay (int) – layer number.

  • icell2d (int) – icell2d number

  • shape (tuple) – shape of the model grid (nlay, ncpl).

Returns:

node – node number

Return type:

int

nlmod.dims.grid.get_row_col_from_xy(x, y, ds, rotated=True, gi=None)[source]

Get the row and column of a point defined by a x and y coordinate.

Parameters:
  • x (float) – The x-coordinate of the point.

  • y (float) – The y-coordinate of the point.

  • ds (xr.Dataset) – The model dataset.

  • rotated (bool, optional) – If the model grid has a rotation, and rotated is False, x and y are in model coordinates. Otherwise x and y are in real world coordinates. The default is True.

  • gi (flopy.utils.GridIntersect, optional) – Can be supplied to use a GridIntersect-class instead of our own calculation

Raises:

ValueError – Raises a ValueError if the point is outside of the model grid.

Returns:

  • row (int) – The row-number of the model cell containing the point (zero-based)

  • col (int) – The column-number of the model cell containing the point (zero-based)

nlmod.dims.grid.get_thickness_from_topbot(top, bot)[source]

Get thickness from data arrays with top and bots.

Parameters:
  • top (xr.DataArray) – raster with top of each cell. dimensions should be (y,x) or (icell2d).

  • bot (xr.DataArray) – raster with bottom of each cell. dimensions should be (layer, y,x) or (layer, icell2d).

Returns:

thickness – raster with thickness of each cell. dimensions should be (layer, y,x) or (layer, icell2d).

Return type:

xr.DataArray

nlmod.dims.grid.get_vertices(ds, vert_per_cid=4, epsilon=0, rotated=False)[source]

Get vertices of a vertex modelgrid from a ds or the modelgrid. Only return the 4 corners of each cell and not the corners of adjacent cells thus limiting the vertices per cell to 4 points.

This method uses the xvertices and yvertices attributes of the modelgrid. When no modelgrid is supplied, a modelgrid-object is created from ds.

Parameters:
  • ds (xr.Dataset) – model dataset, attribute grid_type should be ‘vertex’

  • modelgrid (flopy.discretization.vertexgrid.VertexGrid) – vertex grid with attributes xvertices and yvertices.

  • vert_per_cid (int or None:) –

    number of vertices per cell: - 4 return the 4 vertices of each cell - 5 return the 4 vertices of each cell + one duplicate vertex (sometimes useful if you want to create polygons) - anything else, the maximum number of vertices. For locally refined cells this includes all the vertices adjacent to the cell.

    if vert_per_cid is 4 or 5 vertices are removed using the Ramer-Douglas-Peucker Algorithm -> https://github.com/fhirschmann/rdp.

  • epsilon (int or float, optional) – epsilon in the rdp algorithm. I (Onno) think this is: the maximum distance between a line and a point for which the point is considered to be on the line. The default is 0.

Returns:

vertices_da – Vertex coördinates per cell with dimensions(cid, no_vert, 2).

Return type:

xr.DataArray

nlmod.dims.grid.get_vertices_arr(ds, modelgrid=None, vert_per_cid=4, epsilon=0, rotated=False)[source]

Get vertices of a vertex modelgrid from a ds or the modelgrid. Only return the 4 corners of each cell and not the corners of adjacent cells thus limiting the vertices per cell to 4 points.

This method uses the xvertices and yvertices attributes of the modelgrid. When no modelgrid is supplied, a modelgrid-object is created from ds.

Parameters:
  • ds (xr.Dataset) – model dataset, attribute grid_type should be ‘vertex’

  • modelgrid (flopy.discretization.vertexgrid.VertexGrid) – vertex grid with attributes xvertices and yvertices.

  • vert_per_cid (int or None:) –

    number of vertices per cell: - 4 return the 4 vertices of each cell - 5 return the 4 vertices of each cell + one duplicate vertex (sometimes useful if you want to create polygons) - anything else, the maximum number of vertices. For locally refined cells this includes all the vertices adjacent to the cell.

    if vert_per_cid is 4 or 5 vertices are removed using the Ramer-Douglas-Peucker Algorithm -> https://github.com/fhirschmann/rdp.

  • epsilon (int or float, optional) – epsilon in the rdp algorithm. I (Onno) think this is: the maximum distance between a line and a point for which the point is considered to be on the line. The default is 0.

Returns:

vertices_arr – Vertex coördinates per cell with dimensions(cid, no_vert, 2).

Return type:

numpy array

nlmod.dims.grid.get_vertices_from_ds(ds)[source]

Get the vertices-list from a model dataset.

Flopy needs needs this list to build a disv-package

nlmod.dims.grid.get_xyi_icell2d(gridprops=None, ds=None)[source]

Get x and y coordinates of the cell mids from the cellids in the grid properties.

Parameters:
  • gridprops (dictionary, optional) – dictionary with grid properties output from gridgen. If gridprops is None xyi and icell2d will be obtained from ds.

  • ds (xr.Dataset) – dataset with model data. Should have dimension (layer, icell2d).

Returns:

  • xyi (numpy.ndarray) – array with x and y coördinates of cell centers, shape(len(icell2d), 2).

  • icell2d (numpy.ndarray) – array with cellids, shape(len(icell2d))

nlmod.dims.grid.gridprops_to_vertex_ds(gridprops, ds, nodata=-1)[source]

Gridprops is a dictionary containing keyword arguments needed to generate a flopy modelgrid instance.

nlmod.dims.grid.interpolate_gdf_to_array(gdf, gwf, field='values', method='nearest')[source]

Interpolate data from a point gdf.

Parameters:
  • gdf (geopandas.GeoDataframe) – vector data can only contain a single geometry type.

  • gwf (flopy groundwater flow model) – model with a structured grid.

  • field (str, optional) – column name in the geodataframe. The default is ‘values’.

  • method (str or None, optional) – method to obtain values in cells without geometry by interpolating between cells with values. Options are ‘nearest’ and ‘linear’.

Returns:

arr – numpy array with interpolated data.

Return type:

np.array

nlmod.dims.grid.lcid_to_reclist(layers, cellids, ds, col1=None, col2=None, col3=None, aux=None)[source]

Create a reclist for stress period data from a set of cellids.

Used for vertex grids.

Parameters:
  • layers (list or numpy.ndarray) – list with the layer for each cell in the reclist.

  • cellids (tuple of numpy arrays) –

    tuple with indices of the cells that will be used to create the list with values for a column. There are 2 options:

    1. cellids contains (layers, cids)

    2. cellids contains (cids)

  • ds (xr.Dataset) – dataset with model data. Should have dimensions (layer, icell2d).

  • col1 (str, int or float, optional) – 1st column of the reclist, if None the reclist will be a list with ((layer,icell2d)) for each row. col1 should be the following value for each package (can also be the name of a timeseries): - rch: recharge [L/T] - ghb: head [L] - drn: drain level [L] - chd: head [L] - riv: stage [L]

  • col2 (str, int or float, optional) – 2nd column of the reclist, if None the reclist will be a list with ((layer,icell2d), col1) for each row. col2 should be the following value for each package (can also be the name of a timeseries): - ghb: conductance [L^2/T] - drn: conductance [L^2/T] - riv: conductacnt [L^2/T]

  • col3 (str, int or float, optional) – 3th column of the reclist, if None the reclist will be a list with ((layer,icell2d), col1, col2) for each row. col3 should be the following value for each package (can also be the name of a timeseries): - riv: bottom [L]

  • aux (str or list of str) – list of auxiliary variables to include in reclist

Raises:

ValueError – Question: will this error ever occur?.

Returns:

reclist – every row consist of ((layer, icell2d), col1, col2, col3) grids.

Return type:

list of tuples

nlmod.dims.grid.lrc_to_reclist(layers, rows, columns, cellids, ds, col1=None, col2=None, col3=None, aux=None)[source]

Create a reclist for stress period data from a set of cellids.

Used for structured grids.

Parameters:
  • layers (list or numpy.ndarray) – list with the layer for each cell in the reclist.

  • rows (list or numpy.ndarray) – list with the rows for each cell in the reclist.

  • columns (list or numpy.ndarray) – list with the columns for each cell in the reclist.

  • cellids (tuple of numpy arrays) – tuple with indices of the cells that will be used to create the list with values.

  • ds (xr.Dataset) – dataset with model data. Can have dimension (layer, y, x) or (layer, icell2d).

  • col1 (str, int or float, optional) –

    1st column of the reclist, if None the reclist will be a list with ((layer,row,column)) for each row.

    col1 should be the following value for each package (can also be the

    name of a timeseries): rch: recharge [L/T] ghb: head [L] drn: drain level [L] chd: head [L]

  • col2 (str, int or float, optional) –

    2nd column of the reclist, if None the reclist will be a list with ((layer,row,column), col1) for each row.

    col2 should be the following value for each package (can also be the

    name of a timeseries): ghb: conductance [L^2/T] drn: conductance [L^2/T]

  • col3 (str, int or float, optional) –

    3th column of the reclist, if None the reclist will be a list with ((layer,row,column), col1, col2) for each row.

    col3 should be the following value for each package (can also be the

    name of a timeseries):

  • aux (str or list of str) – list of auxiliary variables to include in reclist

Raises:

ValueError – Question: will this error ever occur?.

Returns:

reclist – every row consist of ((layer,row,column), col1, col2, col3).

Return type:

list of tuples

nlmod.dims.grid.mask_model_edge(ds, idomain=None, cachedir=None, cachename=None)[source]

Get data array which is 1 for every active cell (defined by idomain) at the boundaries of the model (xmin, xmax, ymin, ymax). Other cells are 0.

Parameters:
  • ds (xr.Dataset) – dataset with model data.

  • idomain (xr.DataArray, optional) – idomain used to get active cells and shape of DataArray. Calculate from ds when None. The default is None.

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

ds_out – dataset with edge mask array

Return type:

xr.Dataset

nlmod.dims.grid.modelgrid_from_ds(ds, rotated=True, nlay=None, top=None, botm=None, **kwargs)[source]

Get flopy modelgrid from ds.

Parameters:

ds (xr.Dataset) – model dataset.

Returns:

modelgrid – grid information.

Return type:

StructuredGrid, VertexGrid

nlmod.dims.grid.modelgrid_to_ds(mg=None, grbfile=None)[source]

Create Dataset from flopy modelgrid object.

Parameters:
  • mg (flopy.discretization.Grid) – flopy modelgrid object

  • grbfile (str) – path to a binary grid file

Returns:

ds – Dataset containing grid information

Return type:

xr.Dataset

nlmod.dims.grid.modelgrid_to_vertex_ds(mg, ds, nodata=-1)[source]

Add information about the calculation-grid to a model dataset.

nlmod.dims.grid.node_to_lrc(node, shape)[source]

Convert a node number to (layer,) row and column.

Layer is only returned if shape is 3D.

Parameters:
  • node (int) – node number.

  • shape (tuple) – shape of the model grid.

Returns:

lrc – layer, row and column.

Return type:

tuple

nlmod.dims.grid.polygon_to_area(modelgrid, polygon, da, gridtype='structured')[source]

Create a grid with the surface area in each cell based on a polygon value.

Parameters:
  • modelgrid (flopy.discretization.structuredgrid.StructuredGrid) – grid.

  • polygon (shapely.geometry.polygon.Polygon) – polygon feature.

  • da (xr.DataArray) – data array that is filled with polygon data

Returns:

area_array – area of polygon within each modelgrid cell

Return type:

xr.DataArray

nlmod.dims.grid.polygons_from_ds(ds)[source]

Create polygons of each cell in a model dataset.

Parameters:

ds (xr.Dataset) – Dataset with model data.

Raises:

ValueError – for wrong gridtype or inconsistent grid definition.

Returns:

polygons – list with polygon of each raster cell.

Return type:

list of shapely Polygons

nlmod.dims.grid.refine(ds, model_ws=None, refinement_features=None, exe_name=None, remove_nan_layers=True, model_coordinates=False, version_tag=None)[source]

Refine the grid (discretization by vertices, disv), using Gridgen.

Parameters:
  • ds (xr.Dataset) – A structured model Dataset.

  • model_ws (str, optional) – The working directory for GridGen. Get from ds when model_ws is None. The default is None.

  • refinement_features (list of tuples of length 2 or 3, optional) – List of tuples containing refinement features. Each tuple must be of the form (GeoDataFrame, level) or (geometry, shape_type, level). When refinement_features is None, no refinement is added, but the structured model Dataset is transformed to a Vertex Dataset. The default is None.

  • exe_name (str, optional) – Filepath to the gridgen executable. The file path within nlmod is chose if exe_name is None. The default is None.

  • remove_nan_layers (bool, optional) – if True layers that are inactive everywhere are removed from the model. If False nan layers are kept which might be usefull if you want to keep some layers that exist in other models. The default is True.

  • model_coordinates (bool, optional) – When model_coordinates is True, the features supplied in refinement_features are already in model-coordinates. Only used when a grid is rotated. The default is False.

  • version_tag (str, default None) – GitHub release ID: for example “18.0” or “latest”. If version_tag is provided, the most recent installation location of MODFLOW is found in flopy metadata that respects version_tag. If not found, the executables are downloaded. Not compatible with exe_name.

Returns:

A Vertex model Dataset

Return type:

xr.Dataset

nlmod.dims.grid.snap_extent(extent, delr, delc)[source]

Snap the extent in such a way that an integer number of columns and rows fit in the extent. The new extent is always equal to, or bigger than the original extent.

Parameters:
  • extent (list, tuple or np.array) – original extent (xmin, xmax, ymin, ymax)

  • delr (int or float,) – cell size along rows, equal to dx

  • delc (int or float,) – cell size along columns, equal to dy

Returns:

extent – adjusted extent

Return type:

list, tuple or np.array

nlmod.dims.grid.update_ds_from_layer_ds(ds, layer_ds, method='nearest', **kwargs)[source]

Add variables from a layer Dataset to a model Dataset. Keep de grid- information from the model Dataset (x and y or icell2d), but update the layer dimension when neccesary.

Parameters:
  • ds (xr.Dataset) – dataset with model data. Can have dimension (layer, y, x) or (layer, icell2d).

  • layer_ds (xr.Dataset) – dataset with layer data.

  • method (str) – The method used for resampling layer_ds to the grid of ds.

  • **kwargs (keyword arguments) – keyword arguments are passed to the fill_nan_top_botm_kh_kv-method.

Returns:

ds – Dataset with variables from layer_ds.

Return type:

xr.Dataset

nlmod.dims.grid.xy_to_icell2d(xy, ds)[source]

Get the icell2d value of a point defined by its x and y coordinates.

Parameters:
  • xy (list, tuple) – coordinates of ta point.

  • ds (xr.Dataset) – model dataset.

Returns:

icell2d – number of the icell2d value of a cell containing the xy point.

Return type:

int

nlmod.dims.grid.xy_to_row_col(xy, ds)[source]

Get the row and column values of a point defined by its x and y coordinates.

Parameters:
  • xy (list, tuple) – coordinates of ta point.

  • ds (xr.Dataset) – model dataset.

Returns:

  • row (int) – number of the row value of a cell containing the xy point.

  • col (int) – number of the column value of a cell containing the xy point.

nlmod.dims.grid.xyz_to_cid(xyz, ds=None, modelgrid=None)[source]

Get the icell2d value of a point defined by its x and y coordinates.

Parameters:
  • xyz (list, tuple) – coordinates of a point.

  • ds (xr.Dataset) – model dataset.

  • modelgrid (StructuredGrid, VertexGrid, optional) – A flopy grid-object

Returns:

cid – (layer, cid) for vertex grid, (layer, row, column) for structured grid.

Return type:

tuple

nlmod.dims.layers._add_bathymetry_to_top_bot_kh_kv(ds, bathymetry, fill_mask, kh_sea, kv_sea)[source]

Add bathymetry to the top and bot of each layer for all cells with fill_mask.

This method sets the top of the model at fill_mask to 0 m, and changes the first layer to sea, by setting the botm of this layer to bathymetry, kh to kh_sea and kv to kv_sea. If deeper layers are above bathymetry. the layer depth is set to bathymetry.

Parameters:
  • ds (xarray.Dataset) – dataset with model data, should

  • bathymetry (xarray DataArray) – bathymetry data

  • fill_mask (xr.DataArray) – cell value is 1 if you want to add bathymetry

  • kh_sea (int or float, optional) – the horizontal conductance of the sea cells

  • kv_sea (int or float, optional) – the vertical conductance of the sea cells

Returns:

ds – dataset with model data where the top, bot, kh and kv are changed

Return type:

xarray.Dataset

nlmod.dims.layers._fill_var(var, by, idomain, msg_suffix='')[source]
nlmod.dims.layers._get_empty_layered_da(da, nlay)[source]

Get empty DataArray with number of layer specified by nlay.

Parameters:
  • da (xarray.DataArray) – original data array

  • nlay (int) – number of layers

Returns:

da_new – new empty data array with updated number of layers

Return type:

xarray.DataArray

nlmod.dims.layers._get_isosurface_1d_numba(da: ndarray, z: ndarray, value: float, left: float, right: float) ndarray[source]

Typed wrapper so linters see the correct signature.

nlmod.dims.layers._get_isosurface_1d_numpy(da_arr, z_arr, value, left=nan, right=nan)[source]

Vectorized numpy implementation of get_isosurface_1d.

da_arr, z_arr: ndarrays with layer as the LAST axis. When called via xr.apply_ufunc with input_core_dims=[[“layer”], [“layer”], []], xarray automatically moves the layer dimension to the last axis before calling this function.

Returns array of shape da_arr.shape[:-1].

nlmod.dims.layers._get_isosurface_numba(da, z, value, left=nan, right=nan, **kwargs)[source]

Wrapper for numba implementation of get_isosurface_1d.

This wrapper is needed to move the layer dimension to the last position, as required by the gufunc, and to move the result back to an xarray DataArray with the correct dimensions and coordinates.

Parameters:
  • da (xr.DataArray) – 3D or 4D DataArray with values, e.g. concentration, pressure

  • z (xr.DataArray) – 3D DataArray with elevations

  • value (float) – value at which to compute the elevations of the isosurface

  • left (float, optional) – value to return when value is above the maximum of da. The default is np.nan.

  • right (float, optional) – value to return when value is below the minimum of da. The default is np.nan.

  • kwargs (dict) – additional arguments passed to xarray.apply_ufunc, not used in this function but included for consistency with get_isosurface.

Returns:

2D/3D DataArray with elevations of the isosurface

Return type:

xr.DataArray

nlmod.dims.layers._get_modellayers_dsobs(ds_obs, dimname='n_obs')[source]

Get modellayers from a dataset of observation point data.

Parameters:
  • ds_obs (xr.Dataset) – typically a subset of a model dataset with only data for observations.

  • dimname (str, optional) – name of the observation dimension. ‘icell2d’ is used for vertex grids and ‘n_obs’ for structured grid, by default ‘n_obs’.

Returns:

zero-based indices of modellayers nan if screen above or below model boundaries if screen spans multiple layers, choose layer with most screen length.

Return type:

list of floats

Raises:

ValueError – If any screen top is lower or equal to screen bottom.

nlmod.dims.layers._get_new_layer_name(name, isplit)[source]
nlmod.dims.layers._insert_layer_above(ds, above_layer, name, isplit, mask, top, bot, kh, kv, copy)[source]
nlmod.dims.layers._insert_layer_below(ds, below_layer, name, isplit, mask, top, bot, kh, kv, copy)[source]
nlmod.dims.layers._set_new_layer_values(ds, new_layer_name, mask, top, bot, kh, kv)[source]
nlmod.dims.layers._split_var(ds, var, layer, thickness, fctrs, top, bot, start_suffix_at)[source]

Internal method to split a variable of one layer in multiple layers.

nlmod.dims.layers.add_bathymetry_to_layer_model(ds, datavar_sea='northsea', datavar_bathymetry='bathymetry', kh_sea=10, kv_sea=10)[source]

Add bathymetry to a layer model.

Performs the following steps:

  1. fill top, bot, kh and kv add northsea cell by extrapolation

  2. add bathymetry to the layer model.

  3. remove inactive layers

Parameters:
  • ds (xr.Dataset) – model dataset

  • datavar_northsea (str, optional) – checks if this datavar is available in ds, if not it will create the datavar by calling ‘get_northsea’.

nlmod.dims.layers.add_kh_kv_from_ml_layer_to_ds(ml_layer_ds, ds, anisotropy, fill_value_kh, fill_value_kv)[source]

Add kh and kv from a model layer dataset to the model dataset.

Supports structured and vertex grids.

Parameters:
  • ml_layer_ds (xarray.Dataset) – dataset with model layer data with kh and kv

  • ds (xarray.Dataset) – dataset with model data where kh and kv are added to

  • anisotropy (int or float) – factor to calculate kv from kh or the other way around

  • fill_value_kh (int or float, optional) – use this value for kh if there is no data in the layer model. The default is 1.0.

  • fill_value_kv (int or float, optional) – use this value for kv if there is no data in the layer model. The default is 1.0.

Returns:

ds – dataset with model data with new kh and kv

Return type:

xarray.Dataset

Notes

some model dataset, such as regis, also have ‘c’ and ‘kd’ values. These are ignored at the moment

nlmod.dims.layers.add_layer_dim_to_top(ds, set_non_existing_layers_to_nan=True, copy=True)[source]

Change top from 2d to 3d, setting top and botm to NaN for non-existent layers.

Parameters:
  • ds (xr.DataSet) – model DataSet

  • set_non_existing_layers_to_nan (bool, optional) – If True, set the values of top and botm to . The default is True.

  • copy (bool, optional) – If copy=True, data in the return value is always copied, so the original Dataset is not altered. The default is True.

Returns:

ds – dataset with a layer dimension in top.

Return type:

xarray.Dataset

nlmod.dims.layers.aggregate_by_weighted_mean_to_ds(ds, source_ds, var_name)[source]

Aggregate source data to a model dataset using the weighted mean.

The weighted average per model layer is calculated for the variable in the source dataset. The datasets must have the same grid.

Parameters:
  • ds (xr.Dataset) – model dataset containing layer information (x, y, top, botm)

  • source_ds (xr.Dataset) – dataset containing x, y, top, botm and a data variable to aggregate.

  • var_name (str) – name of the data array to aggregate

Returns:

da – data array containing aggregated values from source dataset

Return type:

xarray.DataArray

Raises:

ValueError – if source_ds does not have a layer dimension

nlmod.dims.layers.calculate_resistance(ds, kv='kv', thickness='thickness', top='top', botm='botm', between_layers=None)[source]

Calculate vertical resistance (c) of model layers from the vertical conductivity (kv) and the thickness.

Parameters:
  • ds (xarray.Dataset) – dataset containing information about top and bottom elevations of layers

  • kv (str, optional) – name of data variable containing vertical conductivity, by default ‘kv’

  • thickness (str, optional) – name of data variable containing thickness, if this data variable does not exists thickness is calculated using top and botm. By default ‘thickness’

  • top (str, optional) – name of data variable containing tops, only used to calculate thickness if not available in dataset. By default “top”

  • botm (str, optional) – name of data variable containing bottoms, only used to calculate thickness if not available in dataset. By default “botm”

  • between_layers (bool, optional) – If True, calculate the resistance between the layers, which MODFLOW uses to calculate the flow. The resistance between two layers is then assigned to the top layer, and the bottom model layer gets a resistance of infinity. If False, calculate the resistance of the layers themselves. The default is True. However, in a future version of nlmod the default will be changed to False.

Returns:

c – DataArray containing vertical resistance (c). NaN where layer thickness is zero

Return type:

xarray.DataArray

nlmod.dims.layers.calculate_thickness(ds, top='top', bot='botm')[source]

Calculate thickness from dataset.

Parameters:
  • ds (xarray.Dataset) – dataset containing information about top and bottom elevations of layers

  • top (str, optional) – name of data variable containing tops, by default “top”

  • bot (str, optional) – name of data variable containing bottoms, by default “botm”

Returns:

thickness – DataArray containing thickness information

Return type:

xarray.DataArray

nlmod.dims.layers.calculate_transmissivity(ds, kh='kh', thickness='thickness', top='top', botm='botm')[source]

Calculate the transmissivity (T) as the product of the horizontal conductance (kh) and the thickness (D).

Parameters:
  • ds (xarray.Dataset) – dataset containing information about top and bottom elevations of layers

  • kh (str, optional) – name of data variable containing horizontal conductivity, by default ‘kh’

  • thickness (str, optional) – name of data variable containing thickness, if this data variable does not exists thickness is calculated using top and botm. By default ‘thickness’

  • top (str, optional) – name of data variable containing tops, only used to calculate thickness if not available in dataset. By default “top”

  • botm (str, optional) – name of data variable containing bottoms, only used to calculate thickness if not available in dataset. By default “botm”

Returns:

T – DataArray containing transmissivity (T). NaN where layer thickness is zero

Return type:

xarray.DataArray

nlmod.dims.layers.check_elevations_consistency(ds)[source]
nlmod.dims.layers.combine_layers_ds(ds, combine_layers, layer='layer', top='top', bot='botm', kh='kh', kv='kv', kD='kD', c='c', return_reindexer=False)[source]

Combine layers in Dataset.

Parameters:
  • ds (xarray.Dataset) – xarray Dataset containing information about layers (layers, top and bot)

  • combine_layers (dict or list of iterables) – dictionary with new layer names as keys and a collection of layer names or indices to merge as values. Alternatively a list of iterables, with each iterable containing strings or integers indicating layers to merge. The new layer will be named after the first of the layers to be merged.

  • layer (str, optional) – name of layer dimension, by default ‘layer’

  • top (str, optional) – name of data variable containing top of layers, by default ‘top’

  • bot (str, optional) – name of data variable containing bottom of layers, by default ‘botm’

  • kh (str, optional) – name of data variable containg horizontal hydraulic conductivity, by default ‘kh’. Not parsed if set to None.

  • kv (str, optional) – name of data variable containg vertical hydraulic conductivity, by default ‘kv’. Not parsed if set to None.

  • kD (str, optional) – name of data variable containg transmissivity or kD, by default ‘kD’. Not parsed if set to None.

  • c (str, optional) – name of data variable containg resistance or c, by default ‘c’. Not parsed if set to None.

  • return_reindexer (bool, optional) – Return a dictionary that can be used to reindex variables from the original layer-dimension to the new layer-dimension when True. The default is False.

Returns:

ds_combine – Dataset with new tops and bottoms taking into account combined layers, and recalculated values for parameters (kh, kv, kD, c).

Return type:

xarray.Dataset

Examples

Given some layer model Dataset with named layers. Specifying which layers to merge can be done the following ways.

As a dictionary:

>>> combine_layers = {
        "new_layer_name": [0, 1]  # as layer indices
        "PZWAz": ["PZWAz2", "PZWAz3", "PZWAz4"],  # as strings
    }

As a list of iterables:

>>> combine_layers = [
        (0, 1),  # as layer indices
        ("PZWAz2", "PZWAz3", "PZWAz4"),  # as strings
    ]

Note

When passing integers to combine_layers, these are always intepreted as the layer index (i.e. starting at 0 and numbered consecutively), and not the layer “name”. If the dataset layer names are integers, only the layer index can be used to specify which layers to merge.

nlmod.dims.layers.convert_to_modflow_top_bot(ds, **kwargs)[source]

Removes the layer dimension from top and fills nans in top and botm.

Alias to remove_layer_dim_from_top

nlmod.dims.layers.convert_to_regis_top_bot(ds, **kwargs)[source]

Adds a layer dimension to top and sets non-existing cells to nan in top and botm.

Alias to add_layer_dim_to_top

nlmod.dims.layers.fill_nan_top_botm(ds)[source]

Remove Nans in non-existent layers in botm and top variables.

The NaNs are removed by setting the value to the top and botm of higher/lower layers that do exist.

Parameters:

ds (xr.DataSet) – model DataSet

Returns:

ds – dataset with filled top and botm data

Return type:

xarray.Dataset

nlmod.dims.layers.fill_nan_top_botm_kh_kv(ds, anisotropy=10.0, fill_value_kh=1.0, fill_value_kv=0.1, remove_nan_layers=True)[source]

Update a model dataset, by removing nans and adding necessary info.

Steps:

  1. Compute top and botm values, by filling nans by data from other layers

  2. Remove inactive layers, with no positive thickness anywhere

  3. Compute kh and kv, filling nans with anisotropy or fill_values

nlmod.dims.layers.fill_top_bot_kh_kv_at_mask(ds, fill_mask)[source]

Fill values in top, bot, kh and kv.

Fill where: 1. the cell is True in fill_mask 2. the cell thickness is greater than 0

Fill values: - top: 0 - bot: minimum of bottom_filled or top - kh: kh_filled if thickness is greater than 0 - kv: kv_filled if thickness is greater than 0

Parameters:
  • ds (xr.DataSet) – model dataset

  • fill_mask (xr.DataArray) – 1 where a cell should be replaced by masked value.

Returns:

ds – model dataset with adjusted data variables: ‘top’, ‘botm’, ‘kh’, ‘kv’

Return type:

xr.DataSet

nlmod.dims.layers.get_first_active_layer(ds, **kwargs)[source]

Get the first active layer in each cell from a model ds.

Parameters:
  • ds (xr.DataSet) – Model Dataset

  • **kwargs (dict) – Kwargs are passed on to get_first_active_layer_from_idomain.

Returns:

first_active_layer – raster in which each cell has the zero based number of the first active layer. Shape can be (y, x) or (icell2d)

Return type:

xr.DataArray

nlmod.dims.layers.get_first_active_layer_from_idomain(idomain, nodata=-999)[source]

Get the first (top) active layer in each cell from the idomain.

Parameters:
  • idomain (xr.DataArray) – idomain. Shape can be (layer, y, x) or (layer, icell2d)

  • nodata (int, optional) – nodata value. used for cells that are inactive in all layers. The default is -999.

Returns:

first_active_layer – raster in which each cell has the zero based number of the first active layer. Shape can be (y, x) or (icell2d)

Return type:

xr.DataArray

nlmod.dims.layers.get_idomain(ds, active_domain='active_domain')[source]

Get idomain from a model Dataset.

Idomain is calculated from the thickness of the layers, and will be 1 for all layers with a positive thickness, and -1 (pass-through) otherwise. Additionally, an “active_domain” DataArray can be applied which can be 2d or 3d. Idomain is set to 0 where “active_domain” is False or 0.

Parameters:
  • ds (xr.Dataset) – The model Dataset.

  • active_domain (str or xr.DataArray, optional) – boolean array indicating the active model domain (True = active). Idomain is set to 0 everywhere else. If passed as str, this variable is taken from ds, if passed as xr.DataArray, this DataArray is used. The default is “active_domain”.

Returns:

ds – DataArray of idomain-variable.

Return type:

xr.DataArray

nlmod.dims.layers.get_isosurface(da, z, value, left=nan, right=nan, method='numba', input_core_dims=None, exclude_dims=None, **kwargs)[source]

Linear interpolation to compute the elevation of an isosurface.

Currently only supports linear interpolation.

Parameters:
  • da (xr.DataArray) – 3D or 4D DataArray with values, e.g. concentration, pressure, etc.

  • z (xr.DataArray) – 3D DataArray with elevations

  • value (float) – value at which to compute the elevations of the isosurface

  • left (float, optional) – value to return when value is above the maximum of da. The default is np.nan.

  • right (float, optional) – value to return when value is below the minimum of da. The default is np.nan.

  • method (str, optional) – method to compute the isosurface. The default is “numba”. Other option is “numpy”. The numba method is usually faster than the numpy method, but the numpy method can be faster for small datasets, and does not require numba to be installed.

  • input_core_dims (list of lists, optional) – list of core dimensions for each input, if not provided assumes core dimensions are any dimensions that are not x, y or icell2d. Example input_core_dims for structured datasets da (“time”, “layer”, “y”, “x”) and z (“layer”, “y”, “x”), and value (float) would be: input_core_dims=[[“layer”], [“layer”], []].

  • exclude_dims (set, optional) – set of dimensions that can change shape in computation. The default is None, which assumes the layer dimension is allowed to change. In the example above this would mean exclude_dims={“layer”}. This will result in the an isosurface for each time step.

  • kwargs (dict) – additional arguments passed to xarray.apply_ufunc

Returns:

2D/3D DataArray with elevations of the isosurface

Return type:

xr.DataArray

nlmod.dims.layers.get_isosurface_1d(da, z, value, left=nan, right=nan)[source]

Linear interpolation to get the first elevation corresponding to value.

This function interpolates linearly along z, if da crosses the given interpolation value at multiple depths, the first elevation is returned.

Note

This function is no longer used in nlmod, but is kept for backward compatibility, and as a reference for the implementations in get_isosurface.

Parameters:
  • da (1d-array) – array of values, e.g. concentration, pressure, etc.

  • z (1d-array) – array of elevations

  • value (float) – value for which to compute the elevations corresponding to value

  • left (float, optional) – value to return when value is below the minimum of da. The default is np.nan.

  • right (float, optional) – value to return when value is above the maximum of da. The default is np.nan.

Returns:

first elevation at which data crosses value

Return type:

float

See also

get_isosurface

generalization of this function to 3D and 4D DataArrays, with

support

_get_isosurface_1d_numpy

vectorized numpy implementation of this function

_get_isosurface_1d_numba

numba implementation of this function

nlmod.dims.layers.get_kh_kv(kh, kv, anisotropy, fill_value_kh=1.0, fill_value_kv=0.1, idomain=None)[source]

Create kh and kv grid data for flopy from existing kh, kv and anistropy grids with nan values (typically from REGIS).

fill nans in kh grid in these steps: 1. take kv and multiply by anisotropy, if this is nan: 2. take fill_value_kh

fill nans in kv grid in these steps: 1. take kh and divide by anisotropy, if this is nan: 2. take fill_value_kv

Supports structured and vertex grids.

Parameters:
  • kh (xarray.DataArray) – kh from regis with nan values shape(nlay, nrow, ncol) or shape(nlay, len(icell2d))

  • kv (xarray.DataArray) – kv from regis with nan values shape(nlay, nrow, ncol) or shape(nlay, len(icell2d))

  • anisotropy (int or float) – factor to calculate kv from kh or the other way around

  • fill_value_kh (int or float, optional) – use this value for kh if there is no data in kh, kv and anisotropy. The default is 1.0.

  • fill_value_kv (int or float, optional) – use this value for kv if there is no data in kv, kh and anisotropy. The default is 1.0.

  • idomain (xarray.DataArray, optional) – The idomain DataArray, used in log-messages, to report the number of active cells that are filled. When idomain is None, the total number of cells that are filled is reported, and not just the active cells. The default is None.

Returns:

  • kh (np.ndarray) – kh without nan values (nlay, nrow, ncol) or shape(nlay, len(icell2d))

  • kv (np.ndarray) – kv without nan values (nlay, nrow, ncol) or shape(nlay, len(icell2d))

nlmod.dims.layers.get_last_active_layer(ds, **kwargs)[source]

Get the last active layer in each cell from a model ds.

Parameters:
  • ds (xr.DataSet) – Model Dataset

  • **kwargs (dict) – Kwargs are passed on to get_last_active_layer_from_idomain.

Returns:

last_active_layer – raster in which each cell has the zero based number of the last active layer. Shape can be (y, x) or (icell2d)

Return type:

xr.DataArray

nlmod.dims.layers.get_last_active_layer_from_idomain(idomain, nodata=-999)[source]

Get the last (bottom) active layer in each cell from the idomain.

Parameters:
  • idomain (xr.DataArray) – idomain. Shape can be (layer, y, x) or (layer, icell2d)

  • nodata (int, optional) – nodata value. used for cells that are inactive in all layers. The default is -999.

Returns:

last_active_layer – raster in which each cell has the zero based number of the last active layer. Shape can be (y, x) or (icell2d)

Return type:

xr.DataArray

nlmod.dims.layers.get_layer_of_z(ds, z, above_model=-999, below_model=-999)[source]

Get the layer of a certain z-value in all cells from a model ds.

Parameters:
  • ds (xr.DataSet) – Model Dataset

  • z (float or xr.DataArray) – The z-value for which the layer is determined

  • above_model (int, optional) – value used for cells where z is above the top of the model. The default is -999.

  • below_model (int, optional) – value used for cells where z is below the top of the model. The default is -999.

Returns:

layer – DataArray with values representing the integer layer index. Shape can be (y, x) or (icell2d)

Return type:

xr.DataArray

nlmod.dims.layers.get_modellayers_indexer(ds, df, x='x', y='y', screen_top='screen_top', screen_bottom='screen_bottom', full_output=False, drop_nan_layers=True)[source]

Get a model layer indexer (dataset) for a dataframe with observation wells.

The dataframe must contain the spatial data of the observation wells, i.e. the x, y coordinates and the screen top and bottom values. The column names corresponding to these values can be specified with keyword arguments.

Note that the returned x and y data (if grid is structured) are the cell coordinates and not the original coordinates of the observation points!

Parameters:
  • ds (xr.Dataset) – Model Dataset containing layer elevations (top, botm).

  • df (pd.DataFrame) – DataFrame with x, y coordinates and screen_top and screen_bottom values.

  • x (str, optional) – name of the x-coordinate column in df. The default is “x”.

  • y (str, optional) – name of the y-coordinate column in df. The default is “y”.

  • screen_top (str, optional) – name of the screen top column in df. The default is “screen_top”.

  • screen_bottom (str, optional) – name of the screen bottom column in df. The default is “screen_bottom”.

  • full_output (bool, optional) – If True, return all variables needed to construct the model layer indexer. If False, return only the variables needed to index a data array.

  • drop_nan_layers (bool, optional) – If True, drop the observation wells for which the model layer cannot be determined. These probably lie above or below the model. The default is True. If False, the model layer indexer will contain NaN values for these wells.

Returns:

obs_ds – Dataset with the following variables:

  • x: x-coordinate of the cell in which the observation point is located

  • y: y-coordinate of the cell in which the observation point is located

  • layer : model layer of each observation point

Optionally returns the following variables if full_output is True:

  • icell2d: cell id of the observation point

  • screen_top: top of the screen

  • screen_bottom: bottom of the screen

  • top: top elevation of the model at the observation point

  • botm: bottom elevations of the model layers at the observation point

  • x_obs_local: local x-coordinate of the observation point (if grid is rotated)

  • y_obs_local: local y-coordinate of the observation point (if grid is rotated)

Return type:

xr.Dataset

Examples

Given some model dataset ds and a dataframe df containing the metadata for multiple observation wells:

>>> idx = nlmod.layers.get_modellayer_indexer(ds, df)

This indexer dataset can be used to obtain the simulated heads for each observation well:

>>> heads.sel(**idx)
nlmod.dims.layers.get_modellayers_screens(ds, screen_top, screen_bottom, xy=None, icell2d=None)[source]

Get the model layer of a well.

Parameters:
  • ds (xarray.Dataset) – xarray Dataset with with top and bottoms, can be structured or vertex.

  • screen_top (np.ndarray of shape(nobs)) – collection of screen top values).

  • screen_bottom (np.ndarray of shape(nobs)) – collection of screen bottom values. Should have the same length as screen_top and xy.

  • xy (np.ndarray of shape(nobs, 2), optional) – list of x,y coordinates

  • icell2d (np.ndarray of shape(nobs), optional) – To speed up the process for vertex grids a list of icell2d indices can be given instead of a list of xy coordinates.

Returns:

zero-based indices of modellayers nan if screen above or below model boundaries if screen spans multiple layers, choose layer with most screen length.

Return type:

list of floats

nlmod.dims.layers.get_zcellcenters(ds)[source]

Calculate the z-coordinates of cell centers. Equivalent of modelgrid.zcellcenters in flopy.

Parameters:

ds (xarray.Dataset) – dataset containing information about layers, including top and botm

Returns:

zcellcenters – data array containing the z-coordinates of cell centers

Return type:

xarray.DataArray

nlmod.dims.layers.insert_layer(ds, name, top, bot, kh=None, kv=None, copy=True)[source]

Inserts a layer in a model Dataset, burning it in an existing layer model.

This method loops over the existing layers, and checks if (part of) the new layer needs to be inserted above the existing layer, and if the top or bottom of the existing layer needs to be altered.

When comparing the height of the new layer with an existing layer, there are 7 options:

1 The new layer is entirely above the existing layer: layer is added completely above existing layer. When the bottom of the new layer is above the top of the existing layer (which can happen for the first layer), this creates a gap in the layer model. In the returned Dataset, this gap is added to the existing layer.

2 part of the new layer lies within an existing layer, but the bottom of the new layer is never below the bottom of the existing layer: layer is added above the existing layer, and the top of existing layer is lowered.

3 there are locations where the new layer is above the bottom of the existing layer, but below the top of the existing layer. The new layer splits the existing layer into two sub-layers. This is not supported (yet) and raises an Exception.

4 part of the new layer lies above the bottom of the existing layer, while at other locations the new layer is below the existing layer. The new layer is split, part of the layer is added above the existing layer, and part of the new layer is added to the layer model in the next iteration(s) (above the next layer).

5 Only the upper part of the new layer overlaps with the existing layer: the layer is not added above the extsing layer, but the bottom of the existing layer is raised because of the overlap.

6 The new layer is below the existing layer everywhere. Nothing happens, move on to the next existing layer.

7 When (part of) the new layer is not added to the layer model after comparison with the last existing layer, the (remaining part of) the new layer is added below the existing layers, at the bottom of the model.

Parameters:
  • ds (xarray.Dataset) – xarray Dataset containing information about layers

  • name (string) – The name of the new layer.

  • top (xr.DataArray) – The top of the new layer.

  • bot (xr.DataArray) – The bottom of the new layer..

  • kh (xr.DataArray, optional) – The horizontal conductivity of the new layer. The default is None.

  • kv (xr.DataArray, optional) – The vertical conductivity of the new layer. The default is None.

  • copy (bool, optional) – If copy=True, data in the return value is always copied, so the original Dataset is not altered. The default is True.

Returns:

ds – xarray Dataset containing the new layer(s)

Return type:

xarray.Dataset

nlmod.dims.layers.kheq_combined_layers(kh, thickness, reindexer)[source]

Calculate equivalent horizontal hydraulic conductivity.

Parameters:
  • kh (xarray.DataArray) – data array containing horizontal hydraulic conductivity

  • thickness (xarray.DataArray) – data array containing thickness per layer

  • reindexer (OrdererDict) – dictionary mapping new layer indices to old layer indices

Returns:

da_kh – data array containing equivalent horizontal hydraulic conductivity for combined layers and original hydraulic conductivity in unmodified layers

Return type:

xarray.DataArray

nlmod.dims.layers.kveq_combined_layers(kv, thickness, reindexer)[source]

Calculate equivalent vertical hydraulic conductivity.

Parameters:
  • kv (xarray.DataArray) – data array containing vertical hydraulic conductivity

  • thickness (xarray.DataArray) – data array containing thickness per layer

  • reindexer (OrdererDict) – dictionary mapping new layer indices to old layer indices

Returns:

da_kv – data array containing equivalent vertical hydraulic conductivity for combined layers and original hydraulic conductivity in unmodified layers

Return type:

xarray.DataArray

nlmod.dims.layers.layer_combine_top_bot(ds, combine_layers, layer='layer', top='top', bot='botm')[source]

Calculate new tops and bottoms for combined layers.

Parameters:
  • ds (xarray.Dataset) – xarray Dataset containing information about layers (layers, top and bot)

  • combine_layers (list of tuple of ints) – list of tuples, with each tuple containing integers indicating layer indices to combine into one layer. E.g. [(0, 1), (2, 3)] will combine layers 0 and 1 into a single layer (with index 0) and layers 2 and 3 into a second layer (with index 1).

  • layer (str, optional) – name of layer dimension, by default ‘layer’

  • top (str, optional) – name of data variable containing top of layers, by default ‘top’

  • bot (str, optional) – name of data variable containing bottom of layers, by default ‘botm’

Returns:

  • new_top, new_bot (xarray.DataArrays) – DataArrays containing new tops and bottoms after splitting layers.

  • reindexer (dict) – dictionary mapping new to old layer indices.

nlmod.dims.layers.remove_inactive_layers(ds)[source]

Remove layers which only contain inactive cells.

Parameters:

ds (xr.Dataset) – The model Dataset.

Returns:

ds – The model Dataset without inactive layers.

Return type:

xr.Dataset

nlmod.dims.layers.remove_layer(ds, layer)[source]

Removes a layer from a Dataset, without changing elevations of other layers.

This will create gaps in the layer model.

nlmod.dims.layers.remove_layer_dim_from_top(ds, check=True, set_non_existing_layers_to_nan=False, inconsistency_threshold=0.0, return_inconsistencies=False, copy=True)[source]

Change top from 3d to 2d, removing NaNs in top and botm in the process.

This method sets variable top to the top of the upper layer (like the definition in MODFLOW). This removes redundant data, as the top of all layers (except the first layer) is equal to the botm of the layer above.

Parameters:
  • ds (xr.DataSet) – model DataSet

  • check (bool, optional) – If True, checks for inconsistencies in the layer model and report to logger as warning. The default is True.

  • bool (set_non_existing_layers_to_nan) – If True, sets the value of the botm-variable to NaN for non-existent layers. This is not recommended, as this might break some procedures in nlmod. The default is False.

  • optional – If True, sets the value of the botm-variable to NaN for non-existent layers. This is not recommended, as this might break some procedures in nlmod. The default is False.

  • inconsistency_threshold (float, optional) – The threshold, above which to report inconsistencies in the layer model (where the botm of a layer is not equal to top of a deeper layer). The default is 0.0.

  • return_inconsistencies (bool, optional) – if True, return the difference between the top of layers and the botm of a deeper layer, at the location where inconsitencies have been reported.

  • copy (bool, optional) – If copy=True, data in the return value is always copied, so the original Dataset is not altered. The default is True.

Returns:

ds – dataset without a layer dimension in top, if return_inconsistencies is False. If return_inconsistencies is True, a numpy arrays with non-equal top and bottoms is returned.

Return type:

xarray.Dataset or numpy.array

nlmod.dims.layers.remove_thin_layers(ds, min_thickness=0.1, update_thickness_every_layer=False, copy=True)[source]

Remove cells with a thickness less than min_thickness (setting the thickness to 0)

The thickness of the removed cells is added to the first active layer below

Parameters:
  • ds (xr,Dataset) – Dataset containing information about layers.

  • min_thickness (float, optional) – THe minimum thickness of a layer. The default is 0.1.

  • update_thickness_every_layer (bool, optional) – If True, loop over the layers, from the top down, and remove thin layers, adding the thickness to the first active layer below. If the thickness of this layer is increased more than min_thickness (even if it originally was thinner than min_thickness), the layer is kept. If update_thickness_every_layer is False, all cells that originally were thinner than min_thickness are removed. The default is False.

  • copy (bool, optional) – If copy=True, data in the return value is always copied, so the original Dataset is not altered. The default is True.

Returns:

ds – Dataset containing information about layers.

Return type:

xr.Dataset

nlmod.dims.layers.set_layer_botm(ds, layer, botm, copy=True)[source]

Set the bottom of a layer.

nlmod.dims.layers.set_layer_thickness(ds, layer, thickness, change='botm', copy=True)[source]

Set the layer thickness by changing the bottom of the layer.

nlmod.dims.layers.set_layer_top(ds, layer, top, copy=True)[source]

Set the top of a layer.

nlmod.dims.layers.set_minimum_layer_thickness(ds, layer, min_thickness, change='botm', copy=True)[source]

Make sure layer has a minimum thickness by lowering the botm of the layer where neccesary.

nlmod.dims.layers.set_model_top(ds, top, min_thickness=0.0, copy=True)[source]

Set the model top, changing layer bottoms when necessary as well.

If the new top is higher than the previous top, the extra thickness is added to the highest layer with a thickness larger than 0.

Parameters:
  • ds (xarray.Dataset) – The model dataset, containing the current top.

  • top (xarray.DataArray) – The new model top, with the same shape as the current top.

  • min_thickness (float, optional) – The minimum thickness of top layers. When the thickness is less than min_thickness, the thicjness is added to a deeper layer. The default is 0.0.

  • copy (bool, optional) – If copy=True, data in the return value is always copied, so the original Dataset is not altered. The default is True.

Returns:

ds – The model dataset, containing the new top.

Return type:

xarray.Dataset

nlmod.dims.layers.set_nan_top_and_botm(ds, copy=True)[source]

Sets Nans for non-existent layers in botm and top variables.

Nans are only added to top when it contains a layer dimension.

Parameters:
  • ds (xr.DataSet) – model DataSet

  • copy (bool, optional) – If copy=True, data in the return value is always copied, so the original Dataset is not altered. The default is True.

Returns:

ds – dataset with nans in top and botm at non-existent layers.

Return type:

xarray.Dataset

nlmod.dims.layers.split_layers_ds(ds, split_dict, layer='layer', top='top', bot='botm', return_reindexer=False, start_suffix_at=1)[source]

Split layers based in Dataset.

Parameters:
  • ds (xarray.Dataset) – xarray Dataset containing information about layers (layers, top and bot)

  • split_dict (dict) – dictionary with name (string) or index (integer) of layers to split as keys. There are two options for the values of the dictionary, to indicate how to split up layer: an iterable of factors. E.g. {‘BXk1’: [1, 3]} will split layer ‘BXk1’ into 2 layers, with the first layer equal to 0.25 of the original thickness and the second layer equal to 0.75 of the original thickness. The second option would be to set the value to the number of layers to split the layer into, e.g. {‘BXk1’: 2}, which is equal to {‘BXk1’: [0.5, 0.5]}.

  • layer (str, optional) – name of layer dimension, by default ‘layer’

  • top (str, optional) – name of data variable containing top of layers, by default ‘top’

  • bot (str, optional) – name of data variable containing bottom of layers, by default ‘botm’

  • return_reindexer (bool, optional) – Return a dictionary that can be used to reindex variables from the original layer-dimension to the new layer-dimension when True. The default is False.

  • start_suffix_at (int, optional) – The suffix that the first splitted layer will receive, for layers that were splitted into multiple sub-layers. The default is 1.

Returns:

ds – Dataset with new tops and bottoms taking into account split layers, and filled data for other variables.

Return type:

xarray.Dataset

nlmod.dims.layers.sum_param_combined_layers(da, reindexer)[source]

Calculate combined layer parameter with sum.

Parameters:
  • da (xarray.DataArray) – data array to calculate combined parameters for

  • reindexer (dict) – dictionary mapping new layer indices to old layer indices

Returns:

da_new – data array containing new parameters for combined layers and old parameters for unmodified layers.

Return type:

xarray.DataArray

nlmod.dims.layers.update_idomain_from_thickness(idomain, thickness, mask)[source]

Get new idomain from thickness in the cells where mask is 1 (or True).

Idomain becomes: - 1: if cell thickness is bigger than 0 - 0: if cell thickness is 0 and it is the top layer - -1: if cell thickness is 0 and the layer is in between active cells

Parameters:
  • idomain (xr.DataArray) – raster with idomain of each cell. dimensions should be (layer, y, x) or (layer, icell2d).

  • thickness (xr.DataArray) – raster with thickness of each cell. dimensions should be (layer, y, x) or (layer, icell2d).

  • mask (xr.DataArray) – raster with ones in cell where the ibound is adjusted. dimensions should be (y, x) or (icell2d).

Returns:

idomain – raster with adjusted idomain of each cell. dimensions should be (layer, y, x) or (layer, icell2d).

Return type:

xr.DataArray

nlmod.dims.time._pd_timestamp_to_cftime(time_pd)[source]

Convert a pandas timestamp into a cftime stamp

Parameters:

time_pd (pd.Timestamp or list of pd.Timestamp) – datetimes

Return type:

cftime.datetime or list of cftime.datetime

nlmod.dims.time.dataframe_to_flopy_timeseries(df, ds=None, package=None, filename=None, time_series_namerecord=None, interpolation_methodrecord='stepwise', append=False, shift_label_from_right_to_left=False)[source]

Convert a pandas DataFrame to a FloPy time series.

This function converts a time-indexed DataFrame into a list of tuples formatted for use with FloPy’s time series input. If a FloPy package is provided, it will either initialize or append a time series package.

Parameters:
  • df (pandas.DataFrame) – DataFrame containing the time series data. The index must be datetime-like if ds is provided.

  • ds (xarray.Dataset, optional) – Model Dataset.

  • package (flopy.mf6.ModflowGwfwel, flopy.mf6.ModflowGwfdrn or similar, optional) – The FloPy package to which the time series should be attached.

  • filename (str, optional) – Filename for the time series package. If not provided, defaults to {package.filename}_ts.

  • time_series_namerecord (list of str, optional) – List of time series names corresponding to DataFrame columns. If not provided, defaults to the DataFrame column names.

  • interpolation_methodrecord (str or list of str, default "stepwise") – Interpolation method(s) to use for the time series. Can be a single string or a list of methods (one for each column).

  • append (bool, default False) – If True, the time series will be appended to an existing time series package. Otherwise, a new package will be initialized.

  • shift_label_from_right_to_left (bool, default False) – If True, sets the index of df from the end of the timestep to the start of the timestep (which is required if interpolation_methodrecord=”stepwise”). The default is False.

Returns:

If package is None, returns a list of (time, value1, value2, …) tuples. If package is provided, returns the created or updated time series package.

Return type:

list of tuple or flopy.mf6.ModflowUtlts

Raises:

AssertionError – If DataFrame contains NaNs or if ds.time is not datetime64.

nlmod.dims.time.ds_time_from_model(gwf)[source]
nlmod.dims.time.ds_time_from_modeltime(modeltime)[source]
nlmod.dims.time.ds_time_idx(t, start_datetime=None, time_units='D', dtype='datetime')[source]

Get time index variable from elapsed time array.

Parameters:
  • t (np.array) – array containing elapsed time, usually in days

  • start_datetime (str, pd.Timestamp, optional) – starting datetime

  • time_units (str, optional) – time units, default is days

  • dtype (str, optional) – dtype of time index. Can be ‘datetime’ or ‘float’. If ‘datetime’ try to create a pandas datetime index if that fails use cftime lib. Default is ‘datetime’.

Returns:

time coordinate for xarray data-array or dataset

Return type:

IndexVariable

nlmod.dims.time.ds_time_idx_from_model(gwf)[source]

Get time index variable from model (gwf or gwt).

Parameters:

gwf (flopy MFModel object) – groundwater flow or groundwater transport model

Returns:

time coordinate for xarray data-array or dataset

Return type:

IndexVariable

nlmod.dims.time.ds_time_idx_from_modeltime(modeltime)[source]

Get time index variable from modeltime object.

Parameters:

modeltime (flopy ModelTime object) – modeltime object (e.g. gwf.modeltime)

Returns:

time coordinate for xarray data-array or dataset

Return type:

IndexVariable

nlmod.dims.time.ds_time_idx_from_tdis_settings(start, perlen, nstp=1, tsmult=1.0, time_units='D')[source]

Get time index from TDIS perioddata: perlen, nstp, tsmult.

Parameters:
  • start (str, pd.Timestamp) – start datetime

  • perlen (array-like) – array of period lengths

  • nstp (int, or array-like optional) – number of steps per period, by default 1

  • tsmult (float or array-like, optional) – timestep multiplier per period, by default 1.0

  • time_units (str, optional) – time units, by default “D”

Returns:

time coordinate for xarray data-array or dataset

Return type:

IndexVariable

nlmod.dims.time.ds_time_to_pandas_index(da_or_ds, include_start=True)[source]

Convert xarray time index to pandas datetime index.

Parameters:
  • da_or_ds (xarray.Dataset or xarray.DataArray) – dataset with time index or time data array

  • include_start (bool, optional) – include the start time in the index, by default True

Returns:

pandas datetime index

Return type:

pd.DatetimeIndex

nlmod.dims.time.estimate_nstp(forcing, perlen=1, tsmult=1.1, nstp_min=1, nstp_max=25, return_dt_arr=False)[source]

Scale the nstp’s linearly between the min and max of the forcing.

Ensures that the first time step of this stress period connects to the last time step of the previous stress period. The ratio between the two time-step durations can be at most tsmult.

Parameters:
  • forcing (array-like) – Array with a forcing value for each stress period. Forcing can be for example a pumping rate or a rainfall intensity.

  • perlen (float or array of floats (nper)) – An array of the stress period lengths.

  • tsmult (float or array of floats (nper)) – Time step multiplier.

  • nstp (int or array of ints (nper)) – Number of time steps in each stress period.

  • nstp_min (int) – nstp value for the stress period with the smallest forcing.

  • nstp_max (int) – nstp value for the stress period with the largest forcing.

Returns:

  • nstp (np.ndarray) – Array with a nstp for each stress period.

  • dt_arr (np.ndarray (optional)) – if return_dt_arr is True returns the durations of the timesteps corresponding with the returned nstp.

nlmod.dims.time.get_perlen(ds)[source]

Get perlen from ds.

Parameters:

ds (xarray.Dataset) – dataset with time variant model data

Returns:

perlen – length of a stress period.

Return type:

list of floats

nlmod.dims.time.get_time_step_length(perlen, nstp, tsmult)[source]

Get the length of the timesteps within a single stress-period.

Parameters:
  • perlen (float) – The length of the stress period, in the time unit of the model (generally days).

  • nstp (int) – The numer of timesteps within the stress period.

  • tsmult (float) – THe time step multiplier, generally equal or lager than 1.

Returns:

t – An array with the length of each of the timesteps within the stress period, in the same unit as perlen.

Return type:

np.ndarray

nlmod.dims.time.set_ds_time(ds, start, time=None, perlen=None, steady=False, steady_start=True, time_units='DAYS', nstp=1, tsmult=1.0)[source]

Set time discretisation for model dataset.

Parameters:
  • ds (xarray.Dataset) – model dataset

  • start (int, float, str, pandas.Timestamp or cftime.datetime) – model start. When start is an integer or float it is interpreted as the number of days of the first stress-period. When start is a string, pandas Timestamp or cftime datetime it is the start datetime of the simulation. Use cftime datetime when you get an OutOfBounds error using pandas.

  • time (float, int or array-like, optional) – float(s) (indicating elapsed time) or timestamp(s) corresponding to the end of each stress period in the model. When time is a single value, the model will have only one stress period. When time is None, the stress period lengths have to be supplied via perlen. The default is None.

  • perlen (float, int or array-like, optional) – length of each stress-period. Only used when time is None. When perlen is a single value, the model will have only one stress period. The default is None.

  • steady (arraylike or bool, optional) – arraylike indicating which stress periods are steady-state, by default False, which sets all stress periods to transient with the first period determined by value of steady_start.

  • steady_start (bool, optional) – whether to set the first period to steady-state, default is True, only used when steady is passed as single boolean.

  • time_units (str, optional) – time units, by default “DAYS”

  • nstp (int or array-like, optional) – number of steps per stress period, stored in ds.attrs, default is 1

  • tsmult (float, optional) – timestep multiplier within stress periods, stored in ds.attrs, default is 1.0

Returns:

ds – model dataset with added time coordinate

Return type:

xarray.Dataset

nlmod.dims.time.set_ds_time_deprecated(ds, time=None, steady_state=False, steady_start=True, steady_start_perlen=3652.0, time_units='DAYS', start_time=None, transient_timesteps=0, perlen=1.0, nstp=1, tsmult=1.0)[source]

Set timing for a model dataset.

Parameters:
  • ds (xarray.Dataset) – Dataset to add time information to

  • time (list, array or DatetimeIndex of pandas.Timestamps, optional) – an array with the start-time of the model and the ending times of the stress-periods. When steady_start is True, the first value of time is the start of the transient model period. The default is None.

  • steady_state (bool, optional) – if True the model is steady state. The default is False.

  • steady_start (bool, optional) – if True the model is transient with a steady state start time step. The default is True.

  • steady_start_perlen (float, optional) – stress-period length of the first steady state stress period. Only used if steady_start is True. The period is used to determine the recharge in this steady state stress period. Default is 3652 days (approximately 10 years).

  • time_units (str, optional) – time unit of the model. The default is ‘DAYS’, which is the only allowed value for now.

  • start_time (str or datetime, optional) – start time of the model. When steady_start is True, this is the start_time of the transient model period. Input is ignored when time is assigned. The default is January 1, 2000.

  • transient_timesteps (int, optional) – number of transient time steps. Only used when steady_state is False. Input is ignored when time is assigned or perlen is an iterable. The default is 0.

  • perlen (float, int, list or np.array, optional) –

    length of each stress-period depending on the type: - float or int: this is the length of all the stress periods. - list or array: the items are the length of the stress-periods in

    days. The length of perlen should match the number of stress-periods (including steady-state stress-periods).

    When steady_start is True, the length of the first steady state stress period is defined by steady_start_perlen. Input for perlen is ignored when time is assigned. The default is 1.0.

  • nstp (int, optional) – number of steps. The default is 1.

  • tsmult (float, optional) – timestep multiplier. The default is 1.0.

Returns:

ds – dataset with time variant model data

Return type:

xarray.Dataset

nlmod.dims.time.set_ds_time_numeric(ds, start, time=None, perlen=None, steady=False, steady_start=True, time_units='DAYS', nstp=1, tsmult=1.0)[source]

Set a numerical time discretisation for a model dataset.

Parameters:
  • ds (xarray.Dataset) – model dataset

  • start (int, float, str, pandas.Timestamp or cftime.datetime) – model start. When start is an integer or float it is interpreted as the number of days of the first stress-period. When start is a string, pandas Timestamp or cftime datetime it is the start datetime of the simulation. Use cftime datetime when you get an OutOfBounds error using pandas.

  • time (float, int or array-like, optional) – float(s) (indicating elapsed time) corresponding to the end of each stress period in the model. When time is a single value, the model will have only one stress period. When time is None, the stress period lengths have to be supplied via perlen. The default is None.

  • perlen (float, int or array-like, optional) – length of each stress-period. Only used when time is None. When perlen is a single value, the model will have only one stress period. The default is None.

  • steady (arraylike or bool, optional) – arraylike indicating which stress periods are steady-state, by default False, which sets all stress periods to transient with the first period determined by value of steady_start.

  • steady_start (bool, optional) – whether to set the first period to steady-state, default is True, only used when steady is passed as single boolean.

  • time_units (str, optional) – time units, by default “DAYS”

  • nstp (int or array-like, optional) – number of steps per stress period, stored in ds.attrs, default is 1

  • tsmult (float, optional) – timestep multiplier within stress periods, stored in ds.attrs, default is 1.0

Returns:

ds – model dataset with added time coordinate

Return type:

xarray.Dataset

nlmod.dims.time.set_time_variables(ds, start, time, steady, steady_start, time_units, nstp, tsmult)[source]

Add data variables: steady, nstp and tsmult, set attributes: start, time_units.

Parameters:
  • ds (xarray.Dataset) – model dataset

  • start (int, float, str, pandas.Timestamp or cftime.datetime) – model start. When start is an integer or float it is interpreted as the number of days of the first stress-period. When start is a string, pandas Timestamp or cftime datetime it is the start datetime of the simulation. Use cftime datetime when you get an OutOfBounds error using pandas.

  • time (array-like, optional) – numerical (indicating elapsed time) or timestamps corresponding to the end of each stress period in the model.

  • steady (arraylike or bool, optional) – arraylike indicating which stress periods are steady-state, by default False, which sets all stress periods to transient with the first period determined by value of steady_start.

  • steady_start (bool, optional) – whether to set the first period to steady-state, default is True, only used when steady is passed as single boolean.

  • time_units (str, optional) – time units, by default “DAYS”

  • nstp (int or array-like, optional) – number of steps per stress period, stored in ds.attrs, default is 1

  • tsmult (float, optional) – timestep multiplier within stress periods, stored in ds.attrs, default is 1.0

Returns:

ds – model dataset with attributes added to the time coordinates

Return type:

xarray.Dataset

nlmod.dims.resample._set_angrot_attributes(extent, xorigin, yorigin, angrot, attrs)[source]

Internal method to set the properties of the grid in an attribute dictionary.

Parameters:
  • extent (list, tuple or np.array of length 4) – extent (xmin, xmax, ymin, ymax) of the desired grid.

  • xorigin (int or float, optional) – lower left x coordinate of the model grid. When angrot == 0, xorigin is added to the first two values of extent. Otherwise it is the x-coordinate of the point the grid is rotated around, and xorigin is added to the Dataset-attributes. The default is 0.0.

  • yorigin (int or float, optional) – lower left y coordinate of the model grid. When angrot == 0, yorigin is added to the last two values of extent. Otherwise it is the y-coordinate of the point the grid is rotated around, and yorigin is added to the Dataset-attributes. The default is 0.0.

  • angrot (int or float, optinal) – the rotation of the grid in counter clockwise degrees. When angrot != 0 the grid is rotated, and all coordinates of the model are in model coordinates. See https://nlmod.readthedocs.io/en/stable/examples/11_grid_rotation.html for more infomation. The default is 0.0.

  • attrs (dict) – Attributes of a model dataset.

Return type:

None.

nlmod.dims.resample.affine_transform_gdf(gdf, affine)[source]

Apply an affine transformation to a geopandas GeoDataFrame.

nlmod.dims.resample.ds_to_structured_grid(ds_in, extent, delr, delc=None, xorigin=0.0, yorigin=0.0, angrot=0.0, method='nearest')[source]

Resample a dataset from a structured grid to a different structured grid.

Parameters:
  • ds_in (xarray.Dataset) – dataset with dimensions (layer, y, x). y and x are from the original grid

  • extent (list, tuple or np.array of length 4) – extent (xmin, xmax, ymin, ymax) of the desired grid.

  • delr (int or float) – cell size along rows of the desired grid (dx).

  • delc (int or float) – cell size along columns of the desired grid (dy).

  • xorigin (int or float, optional) – lower left x coordinate of the model grid. When angrot == 0, xorigin is added to the first two values of extent. Otherwise it is the x-coordinate of the point the grid is rotated around, and xorigin is added to the Dataset-attributes. The default is 0.0.

  • yorigin (int or float, optional) – lower left y coordinate of the model grid. When angrot == 0, yorigin is added to the last two values of extent. Otherwise it is the y-coordinate of the point the grid is rotated around, and yorigin is added to the Dataset-attributes. The default is 0.0.

  • angrot (int or float, optinal) – the rotation of the grid in counter clockwise degrees. When angrot != 0 the grid is rotated, and all coordinates of the model are in model coordinates. See https://nlmod.readthedocs.io/en/stable/examples/11_grid_rotation.html for more infomation. The default is 0.0.

  • method (str, optional) – type of interpolation used to resample. Sea structured_da_to_ds for possible values of method. The default is ‘nearest’.

Returns:

ds_out – dataset with dimensions (layer, y, x). y and x are from the new grid.

Return type:

xarray.Dataset

nlmod.dims.resample.extent_to_polygon(extent)[source]

Convert an extent to a shapely Polygon.

nlmod.dims.resample.fillnan_da(da, ds=None, method='nearest')[source]

Fill not-a-number values in a DataArray.

Note that extrapolation is only applied if method is “nearest”.

Parameters:
  • da (xr.DataArray) – data array with nan values.

  • ds (xr.Dataset, optional) – Dataset containing grid-properties. Needed when a Vertex grid is used.

  • method (str, optional) – method used in scipy.interpolate.griddata to resample. The default is nearest.

Returns:

xar_out – data array without nan values.

Return type:

xr.DataArray

Notes

can be slow if the xar_in is a large raster

nlmod.dims.resample.fillnan_da_structured_grid(xar_in, method='nearest')[source]

Fill not-a-number values in a structured grid, DataArray.

The nans are replaced with the interpolated values using distance_transform_edt when all cells have the same shape and method is nearest, otherwise the values are interpolated using griddata.

Note that extrapolation is only applied if method is “nearest”.

Parameters:
  • xar_in (xarray DataArray) – DataArray with nan values. DataArray should at least have dimensions x and y.

  • method (str, optional) – method used in scipy.interpolate.griddata to resample, default is nearest.

Returns:

xar_out – DataArray without nan values. DataArray has 2 dimensions (y and x)

Return type:

xarray DataArray

nlmod.dims.resample.fillnan_da_vertex_grid(xar_in, ds=None, x=None, y=None, method='nearest')[source]

Fill not-a-number values in a vertex grid, DataArray.

Note that extrapolation is only applied if method is “nearest”.

Parameters:
  • xar_in (xr.DataArray) – data array with nan values. Shape is (icell2d)

  • ds (xr.Dataset) – Dataset containing grid-properties

  • x (np.array, optional) – x coördinates of the cell centers shape(icell2d), if not defined use x from ds.

  • y (np.array, optional) – y coördinates of the cell centers shape(icell2d), if not defined use y from ds.

  • method (str, optional) – method used in scipy.interpolate.griddata to resample, default is nearest.

Returns:

xar_out – data array without nan values. Shape is (icell2d)

Return type:

xr.DataArray

Notes

If x is provided, x will be used over ds and the x coordinate part of xar_in. If x is not provided, ds will be used to get the x coordinates. If x is not provided and “x” is not in xar_in.coords, an error will be raised.

nlmod.dims.resample.get_affine(ds, sx=None, sy=None)[source]

Get the affine-transformation, from pixel to real-world coordinates.

nlmod.dims.resample.get_affine_mod_to_world(ds)[source]

Get the affine-transformation from model to real-world coordinates.

nlmod.dims.resample.get_affine_world_to_mod(ds)[source]

Get the affine-transformation from real-world to model coordinates.

nlmod.dims.resample.get_extent(ds, rotated=True)[source]

Get the model extent, corrected for angrot if necessary.

nlmod.dims.resample.get_extent_polygon(ds, rotated=True)[source]

Get the model extent, as a shapely Polygon.

nlmod.dims.resample.get_xy_mid_structured(extent, delr, delc, descending_y=True, tiny=1e-10)[source]

Calculates the x and y coordinates of the cell centers of a structured grid.

Parameters:
  • extent (list, tuple or np.array) – extent (xmin, xmax, ymin, ymax)

  • delr (int, float, list, tuple or array, optional) – The gridsize along columns (dx). The default is 100. meter.

  • delc (None, int, float, list, tuple or array, optional) – The gridsize along rows (dy). Set to delr when None. If None delc=delr

  • descending_y (bool, optional) – if True the resulting ymid array is in descending order. This is the default for MODFLOW models. default is True.

  • tiny (float, optional) – A small value (for floating point reasons) to check if the extent is valid. The default is 1e-10.

Returns:

  • x (np.array) – x-coordinates of the cell centers shape (ncol)

  • y (np.array) – y-coordinates of the cell centers shape (nrow)

nlmod.dims.resample.structured_da_to_ds(da, ds, method='average', nodata=nan)[source]

Resample a DataArray to the coordinates of a model dataset.

Parameters:
  • da (xarray.DataArray) – The data-array to be resampled, with dimensions x and y.

  • ds (xarray.Dataset) – The model dataset.

  • method (string or rasterio.enums.Resampling, optional) – The method to resample the DataArray. Possible values are “linear”, “nearest” and all the values in rasterio.enums.Resampling. These values can be provided as a string (‘average’) or as an attribute of rasterio.enums.Resampling (rasterio.enums.Resampling.average). When method is ‘linear’ or ‘nearest’ da.interp() is used. Otherwise da.rio.reproject_match() is used. The default is “average”.

  • nodata (float, optional) – The nodata value in input and output. The default is np.NaN.

Returns:

da_out – The resampled DataArray

Return type:

xarray.DataArray

nlmod.dims.resample.vertex_da_to_ds(da, ds, method='nearest')[source]

Resample a vertex DataArray to a model dataset.

Parameters:
  • da (xaray.DataArray) – A vertex DataArray. When the DataArray does not have ‘icell2d’ as a dimension, the original DataArray is retured. The DataArray da can contain other dimensions as well (for example ‘layer’ or time’’ ).

  • ds (xarray.Dataset) – The model dataset to which the DataArray needs to be resampled, with coordinates x and y.

  • method (str, optional) – The interpolation method, see griddata. The default is “nearest”.

Returns:

A DataArray, with the same gridtype as ds.

Return type:

xarray.DataArray

sim

nlmod.sim.sim.ems(sim, pname='ems', model=None, **kwargs)[source]

Create explicit model solution (EMS) package.

Parameters:
  • sim (flopy MFSimulation) – simulation object.

  • pname (str, optional) – package name

nlmod.sim.sim.get_tdis_perioddata(ds, nstp='nstp', tsmult='tsmult')[source]

Get tdis_perioddata from ds.

dsxarray.Dataset

dataset with time variant model data

tdis_perioddata[perlen, nstp, tsmult]
  • perlen (double) is the length of a stress period.

  • nstp (integer) is the number of time steps in a stress period.

  • tsmult (double) is the multiplier for the length of successive time steps. The length of a time step is calculated by multiplying the length of the previous time step by TSMULT. The length of the first time step, \(\Delta t_1\), is related to PERLEN, NSTP, and TSMULT by the relation :math:`Delta t_1= perlen

rac{tsmult -

1}{tsmult^{nstp}-1}`.

nlmod.sim.sim.ims(sim, complexity='MODERATE', pname='ims', model=None, **kwargs)[source]

Create implicit model solution (IMS) package.

Parameters:
  • sim (flopy MFSimulation) – simulation object.

  • complexity (str, optional) – solver complexity for default settings. The default is “MODERATE”.

  • pname (str, optional) – package name

Returns:

ims – ims object.

Return type:

flopy ModflowIms

nlmod.sim.sim.register_ims_package(sim, model, ims)[source]
nlmod.sim.sim.register_solution_package(sim, model, solver)[source]
nlmod.sim.sim.sim(ds, exe_name=None, version_tag=None, **kwargs)[source]

Create sim from the model dataset.

Parameters:
  • ds (xarray.Dataset) – dataset with model data. Should have the dimension ‘time’ and the attributes: model_name, mfversion, model_ws, time_units, start, perlen, nstp, tsmult

  • exe_name (str, optional) – path to modflow executable, default is None. If None, the path is obtained from the flopy metadata that respects version_tag. If not found, the executables are downloaded. Not compatible with version_tag.

  • version_tag (str, default None) – GitHub release ID: for example “18.0” or “latest”. If version_tag is provided, the most recent installation location of MODFLOW is found in flopy metadata that respects version_tag. If not found, the executables are downloaded. Not compatible with exe_name.

Returns:

sim – simulation object.

Return type:

flopy MFSimulation

nlmod.sim.sim.tdis(ds, sim, pname='tdis', nstp='nstp', tsmult='tsmult', **kwargs)[source]

Create tdis package from the model dataset.

Parameters:
  • ds (xarray.Dataset) – dataset with model data. Should have the dimension ‘time’ and the attributes: time_units, start, perlen, nstp, tsmult

  • sim (flopy MFSimulation) – simulation object.

  • pname (str, optional) – package name

  • **kwargs – passed on to flopy.mft.ModflowTdis

Returns:

dis – tdis object.

Return type:

flopy TDis

nlmod.sim.sim.write_and_run(sim, ds, write_ds=True, script_path=None, silent=False)[source]

Write modflow files and run the model. Extra options include writing the model dataset to a netcdf file in the model workspace and copying the modelscript to the model workspace.

Parameters:
  • sim (flopy.mf6.MFSimulation or flopy.mf6.ModflowGwf) – MF6 Simulation or MF6 Groundwater Flow object.

  • ds (xarray.Dataset) – dataset with model data.

  • write_ds (bool, optional) – if True the model dataset is written to a NetCDF-file (.nc) in the model workspace the name of the .nc file is used from the attribute “model_name”. The default is True.

  • script_path (str or None, optional) – full path of the Jupyter Notebook (.ipynb) or the module (.py) with the modelscript. The default is None. Preferably this path does not have to be given manually but there is currently no good option to obtain the filename of a Jupyter Notebook from within the notebook itself.

  • silent (bool, optional) – write and run model silently

gwf

nlmod.gwf.gwf._dis(ds, model, length_units='METERS', pname='dis', **kwargs)[source]

Create discretisation package from the model dataset.

Parameters:
  • ds (xarray.Dataset) – dataset with model data.

  • model (flopy ModflowGwf or flopy ModflowGwt) – groundwaterflow or groundwater transport object

  • length_units (str, optional) – length unit. The default is ‘METERS’.

  • pname (str, optional) – package name

Returns:

dis – discretisation package.

Return type:

flopy ModflowGwfdis or flopy ModflowGwtdis

nlmod.gwf.gwf._disv(ds, model, length_units='METERS', pname='disv', **kwargs)[source]

Create discretisation vertices package from the model dataset.

Parameters:
  • ds (xarray.Dataset) – dataset with model data.

  • model (flopy ModflowGwf or flopy ModflowGwt) – groundwater flow or groundwater transport object.

  • length_units (str, optional) – length unit. The default is ‘METERS’.

  • pname (str, optional) – package name

Returns:

disv – disv package

Return type:

flopy ModflowGwfdisv or flopy ModflowGwtdisv

nlmod.gwf.gwf._set_record(out, budget, output='head')[source]
nlmod.gwf.gwf.buy(ds, gwf, pname='buy', gwt_modelname=None, **kwargs)[source]

Create buoyancy package from model dataset.

Parameters:
  • ds (xarray.Dataset) – dataset with model data.

  • gwf (flopy ModflowGwf) – groundwaterflow object.

  • pname (str, optional) – package name, by default “buy”

  • gwt_modelname (str, optional) – name of the groundwater transport model, by default None, which uses gwf modelname with ‘_gwt’ appended.

Returns:

buy – buy package

Return type:

flopy ModflowGwfbuy

Raises:

ValueError – if transport is not

nlmod.gwf.gwf.chd(ds, gwf, mask='chd_mask', head='chd_head', pname='chd', auxiliary=None, layer=0, **kwargs)[source]

Create constant head package from model dataset.

Parameters:
  • ds (xarray.Dataset) – dataset with model data.

  • gwf (flopy ModflowGwf) – groundwaterflow object.

  • mask (str, optional) – name of data variable in ds that is 1 for cells with a constant head and zero for all other cells. The default is ‘chd_mask’.

  • head (str, optional) – name of data variable in ds that is used as the head in the chd cells. By default, assumes head data is stored as ‘chd_head’.

  • pname (str, optional) – package name

  • auxiliary (str or list of str) – name(s) of data arrays to include as auxiliary data to reclist

  • layer (int or None) – The layer in which the boundary is added. It is added to the first active layer when layer is None. The default is 0.

  • chd (str, optional) – deprecated, the new argument is ‘mask’

Returns:

chd – chd package

Return type:

flopy ModflowGwfchd

nlmod.gwf.gwf.dis(ds, gwf, length_units='METERS', pname='dis', **kwargs)[source]

Create discretisation package from the model dataset.

Parameters:
  • ds (xarray.Dataset) – dataset with model data.

  • gwf (flopy ModflowGwf) – groundwaterflow object

  • length_units (str, optional) – length unit. The default is ‘METERS’.

  • pname (str, optional) – package name, ignored if ds has a vertex grid (disv)

Returns:

dis – discretisation package.

Return type:

flopy ModflowGwfdis

nlmod.gwf.gwf.disv(ds, gwf, length_units='METERS', pname='disv', **kwargs)[source]

Create discretisation vertices package from the model dataset.

Parameters:
  • ds (xarray.Dataset) – dataset with model data.

  • model (flopy ModflowGwf) – groundwater flow object.

  • length_units (str, optional) – length unit. The default is ‘METERS’.

  • pname (str, optional) – package name

Returns:

disv – disv package

Return type:

flopy ModflowGwfdisv

nlmod.gwf.gwf.drn(ds, gwf, elev='drn_elev', cond='drn_cond', pname='drn', layer=None, **kwargs)[source]

Create drain from model dataset.

Parameters:
  • ds (xarray.Dataset) – dataset with model data.

  • gwf (flopy ModflowGwf) – groundwaterflow object.

  • elev (str or xarray.DataArray, optional) – drain elevation, either as string pointing to data array in ds or as data array. By default assumes data array is stored under “drn_elev”.

  • cond (str or xarray.DataArray, optional) – drain conductance, either as string pointing to data array in ds or as data array. By default assumes data array is stored under “drn_cond”.

  • da_name (str, deprecated) – this is deprecated, name of the drn files in the model dataset

  • pname (str, optional) – package name

  • layer (int or None) – The layer in which the boundary is added. It is added to the first active layer when layer is None. The default is None.

Returns:

drn – drn package

Return type:

flopy ModflowGwfdrn

nlmod.gwf.gwf.ds_to_gwf(ds, complexity='SIMPLE', icelltype=0, under_relaxation=False)[source]

Generate Simulation and GWF model from model DataSet.

Builds the following packages: - sim - tdis - ims - gwf

  • dis

  • npf

  • ic

  • oc

  • rch if “recharge” is present in DataSet

  • evt if “evaporation” is present in DataSet

Parameters:

ds (xarray.Dataset) – model dataset

Returns:

MODFLOW6 GroundwaterFlow model object.

Return type:

flopy.mf6.ModflowGwf

nlmod.gwf.gwf.evt(ds, gwf, pname='evt', **kwargs)[source]

Create evapotranspiration package from model dataset.

Parameters:
  • ds (xarray.Dataset) – dataset with model data.

  • gwf (flopy ModflowGwf) – groundwaterflow object.

  • pname (str, optional) – package name

Returns:

evt – rch package

Return type:

flopy ModflowGwfevt

nlmod.gwf.gwf.ghb(ds, gwf, bhead=None, cond=None, pname='ghb', auxiliary=None, layer=None, **kwargs)[source]

Create general head boundary from model dataset.

Parameters:
  • ds (xarray.Dataset) – dataset with model data.

  • gwf (flopy ModflowGwf) – groundwaterflow object.

  • bhead (str or xarray.DataArray, optional) – ghb boundary head, either as string pointing to data array in ds or as data array. By default None, which assumes data array is stored under “ghb_bhead”.

  • cond (str or xarray.DataArray, optional) – ghb conductance, either as string pointing to data array in ds or as data array. By default None, which assumes data array is stored under “ghb_cond”.

  • pname (str, optional) – package name

  • auxiliary (str or list of str) – name(s) of data arrays to include as auxiliary data to reclist

  • layer (int or None) – The layer in which the boundary is added. It is added to the first active layer when layer is None. The default is None.

Raises:

ValueError – raised if gridtype is not structured or vertex.

Returns:

ghb – ghb package

Return type:

flopy ModflowGwfghb

nlmod.gwf.gwf.gwf(ds, sim, under_relaxation=False, **kwargs)[source]

Create groundwater flow model from the model dataset.

Parameters:
  • ds (xarray.Dataset) – dataset with model data. Should have the dimension ‘time’ and the attributes: model_name, mfversion, model_ws, time_units, start, perlen, nstp, tsmult

  • sim (flopy MFSimulation) – simulation object.

Returns:

gwf – groundwaterflow object.

Return type:

flopy ModflowGwf

nlmod.gwf.gwf.ic(ds, gwf, starting_head='starting_head', pname='ic', **kwargs)[source]

Create initial condictions package from model dataset.

Parameters:
  • ds (xarray.Dataset) – dataset with model data.

  • gwf (flopy ModflowGwf) – groundwaterflow object.

  • starting_head (str, float or int, optional) – if type is int or float this is the starting head for all cells If the type is str the data variable from ds is used as starting head. The default is ‘starting_head’.

  • pname (str, optional) – package name

Returns:

ic – ic package

Return type:

flopy ModflowGwfic

nlmod.gwf.gwf.npf(ds, gwf, k='kh', k33='kv', icelltype=0, save_flows=False, pname='npf', **kwargs)[source]

Create node property flow package from model dataset.

Parameters:
  • ds (xarray.Dataset) – dataset with model data.

  • gwf (flopy ModflowGwf) – groundwaterflow object.

  • icelltype (int or str, optional) – celltype, if int the icelltype for all layer, if str the icelltype from the model ds is used. The default is 0.

  • k (str or array-like) – horizontal hydraulic conductivity, when passed as string, the array is obtained from ds. By default assumes data is stored as “kh”.

  • k33 (str or array-like) – vertical hydraulic conductivity, when passed as string, the array is obtained from ds. By default assumes data is stored as “kv”.

  • save_flows (bool, optional) – value is passed to flopy.mf6.ModflowGwfnpf() to determine if cell by cell flows should be saved to the cbb file. Default is False

  • pname (str, optional) – package name

Raises:

NotImplementedError – only icelltype 0 is implemented.

Returns:

npf – npf package.

Return type:

flopy ModflowGwfnpf

nlmod.gwf.gwf.oc(ds, gwf, save_head=True, save_budget=True, print_head=False, print_budget=False, pname='oc', **kwargs)[source]

Create output control package from model dataset.

Parameters:
  • ds (xarray.Dataset) – dataset with model data.

  • gwf (flopy ModflowGwf) – groundwaterflow object.

  • save_head (bool or str) – Saves the head to the output-file. If save_head is a string, it needs to be “all”, “first” or “last”. If save_head is True, it is set to “last”. The default is True.

  • save_budget (bool or str) – Saves the budgets to the output-file. If save_budget is a string, it needs to be “all”, “first” or “last”. If save_budget is True, it is set to “last”. The default is True.

  • print_head (bool or str) – Prints the head to the list-file. If print_head is a string, it needs to be “all”, “first” or “last”. If print_head is True, it is set to “last”. The default is False.

  • print_budget (bool or str) – Prints the budgets to the list-file. If print_budget is a string, it needs to be “all”, “first” or “last”. If print_budget is True, it is set to “last”. The default is False.

  • pname (str, optional) – package name

Returns:

oc – oc package

Return type:

flopy ModflowGwfoc

nlmod.gwf.gwf.rch(ds, gwf, pname='rch', **kwargs)[source]

Create recharge package from model dataset.

Parameters:
  • ds (xarray.Dataset) – dataset with model data.

  • gwf (flopy ModflowGwf) – groundwaterflow object.

  • pname (str, optional) – package name

Returns:

rch – rch package

Return type:

flopy ModflowGwfrch

nlmod.gwf.gwf.riv(ds, gwf, stage='riv_stage', cond='riv_cond', rbot='riv_rbot', pname='riv', auxiliary=None, layer=None, **kwargs)[source]

Create river package from model dataset.

Parameters:
  • ds (xarray.Dataset) – dataset with model data.

  • gwf (flopy ModflowGwf) – groundwaterflow object.

  • stage (str or xarray.DataArray, optional) – river stage, either as string pointing to data array in ds or as data array. By default assumes data array is stored under “riv_stage”.

  • cond (str or xarray.DataArray, optional) – river conductance, either as string pointing to data array in ds or as data array. By default assumes data array is stored under “riv_cond”.

  • rbot (str or xarray.DataArray, optional) – river bottom elevation, either as string pointing to data array in ds or as data array. By default assumes data array is stored under “riv_rbot”.

  • pname (str, optional) – package name

  • auxiliary (str or list of str) – name(s) of data arrays to include as auxiliary data to reclist

  • layer (int or None) – The layer in which the boundary is added. It is added to the first active layer when layer is None. The default is None.

Returns:

riv – riv package

Return type:

flopy ModflowGwfriv

nlmod.gwf.gwf.sto(ds, gwf, sy='sy', ss='ss', iconvert=1, save_flows=False, pname='sto', **kwargs)[source]

Create storage package from model dataset.

Parameters:
  • ds (xarray.Dataset) – dataset with model data.

  • gwf (flopy ModflowGwf) – groundwaterflow object.

  • sy (str or float, optional) – specific yield. The default is “sy”, or 0.2 if “sy” is not in ds.

  • ss (str or float, optional) – specific storage. The default is “ss”, or 0.00001 if “ss” is not in ds.

  • iconvert (int, optional) – See description in ModflowGwfsto. The default is 1 (differs from FloPY).

  • save_flows (bool, optional) – value is passed to flopy.mf6.ModflowGwfsto() to determine if flows should be saved to the cbb file. Default is False

  • pname (str, optional) – package name

Returns:

sto – sto package

Return type:

flopy ModflowGwfsto

nlmod.gwf.gwf.surface_drain_from_ds(ds, gwf, resistance, elev='ahn', pname='drn', **kwargs)[source]

Create surface level drain (maaivelddrainage in Dutch) from the model dataset.

Parameters:
  • ds (xarray.Dataset) – dataset with model data.

  • gwf (flopy ModflowGwf) – groundwaterflow object.

  • resistance (int or float) – resistance of the surface drain. This value is used to calculate drain conductance by scaling with cell area.

  • elev (str or xarray.DataArray) – name pointing to the data array containing surface drain elevation data, or pass the data array directly. By default assumes the elevation data is stored under “ahn”.

  • pname (str, optional) – package name

Returns:

drn – drn package

Return type:

flopy ModflowGwfdrn

nlmod.gwf.gwf.uzf(ds, gwf, pname='uzf', **kwargs)[source]

Create unsaturated zone flow package from model dataset.

Parameters:
  • ds (xarray.Dataset) – dataset with model data.

  • gwf (flopy ModflowGwf) – groundwaterflow object.

  • pname (str, optional) – package name

Returns:

uzf – uzf package

Return type:

flopy ModflowGwfuzf

nlmod.gwf.recharge._add_time_series(package, df, ds)[source]

Add time series to a package.

Parameters:
  • rch (mfpackage.MFPackage) – The Flopy package to which to add the timeseries.

  • df (pd.DataFrame) – A pandas DataFrane that contains the time series values.

  • ds (xr.Dataset) – The model Dataset. It is used to get the time of the time series.

Return type:

None.

nlmod.gwf.recharge._get_meteo_da_from_input(recharge, ds, pname, stn_var)[source]

Normalize meteorological input for use in package stress-period data.

This helper accepts several input forms for a meteorological variable (e.g. recharge, evaporation, finf/pet for UZF) and returns:

  • a spatial DataArray that describes either a per-cell scalar value or a per-cell reference to a timeseries name,

  • a boolean mask DataArray indicating which cells should receive the meteorological input,

  • a pandas DataFrame of unique timeseries values when time series are used (or None otherwise).

Supported input types and behavior - str:

Treated as the name of a variable in ds (i.e. recharge = ds[recharge]). Further processing follows the xr.DataArray rules below.

  • xr.DataArray:
    • 1-D time series (dims == (“time”,)):

      Interpreted as a single time series applied to all active cells. The returned recharge DataArray contains the timeseries name for every cell (string values), mask_recharge marks active cells, and rch_unique_df is a DataFrame with that single series.

    • 2-D time-by-station (dims == (“time”, “stn_*”)):

      Interpreted as a set of timeseries, one per station. ds[stn_var] is expected to map each spatial cell to a station index. The returned recharge contains per-cell timeseries names constructed from the station index, mask_recharge marks cells with a valid station mapping, and rch_unique_df is the DataFrame of station timeseries (column names are prefixed with pname_).

    • Per-cell (spatial) array with or without time dimension:

      If the array contains a time dimension (transient, time in dims and more than one period), unique time series are identified using _get_unique_series. The returned recharge will then be a per-cell DataArray of timeseries names and rch_unique_df contains the corresponding series. If no time dimension (or single time index), recharge is reduced to a per-cell scalar array and rch_unique_df is None. NaN values in the active domain are not allowed (ValueError).

  • float:

    Treated as a spatially-constant value applied to all active cells; mask_recharge marks active cells and no timeseries DataFrame is returned.

Parameters:
  • recharge (float, str, or xr.DataArray) – The input meteorological data. See “Supported input types and behavior” above for interpretation rules.

  • ds (xr.Dataset) – Model dataset. Used to determine active cells and grid dimensions.

  • pname (str) – Package name used as prefix when constructing timeseries names.

  • stn_var (str) – Name of the DataArray in ds that maps cells to station indices. Used only when recharge is a time-by-station DataArray.

Returns:

  • recharge (xr.DataArray) – If timeseries are used, a spatial DataArray of dtype string containing the timeseries name for each cell (e.g. “rch_0”, “rch_1”, …). Otherwise a per-cell scalar DataArray (no time dimension) with the meteorological value(s).

  • mask_recharge (xr.DataArray (bool)) – Boolean mask indicating which cells should receive the input (True means apply input). Typically this is based on the first active layer.

  • rch_unique_df (pd.DataFrame or None) – When one or more unique timeseries are detected, a DataFrame indexed by model ds.time containing each unique series is returned. If no timeseries are used, returns None.

Raises:
  • ValueError – If per-cell data contains NaN values in the active model domain.

  • NotImplementedError – If recharge is of an unsupported type.

Notes

  • Timeseries names are constructed using pname as a prefix (e.g. “evt_0”, “finf_1”).

nlmod.gwf.recharge._get_unique_series(ds, da, pname)[source]

Get the location and values of unique time series from a variable var in ds.

Parameters:
  • ds (xr.Dataset) – The model Dataset.

  • da (xr.DataArray) – The 3d (structured) or 2d (vertext) DataArray that contains the timeseries.

  • pname (str) – Package name, which is used for the name of the time series.

Raises:

ValueError – DESCRIPTION.

Returns:

  • rch_name_da (xr.DataArray) – The name of the recharge series for each of the cells.

  • rch_unique_df (pd.DataFrame) – The values of each of the time series.

nlmod.gwf.recharge.ds_to_evt(gwf, ds, mask=None, pname='evt', rate='evaporation', nseg=1, surface=None, depth=None, auxiliary=None, **kwargs)[source]

Convert evaporation data in the model dataset to a evt package with time series.

Parameters:
  • gwf (flopy.mf6.modflow.mfgwf.ModflowGwf) – groundwater flow model.

  • ds (xr.DataSet) – dataset containing relevant model grid information

  • mask (xr.DataArray) – data array containing mask, evt is only added where mask is True

  • pname (str, optional) – package name. The default is ‘evt’.

  • rate (str, xr.DataArray or float, optional) – When rate is a string, it is the name of the variable in ds that contains the maximum ET flux rate. When rate is a float, it is interpreted as the constant evaporation rate that is applied in all active cells (within mask if mask is supplied). The default is “evaporation”.

  • nseg (int, optional) – number of ET segments. Only 1 is supported for now. The default is 1.

  • surface (str, float or xr.DataArray, optional) – The elevation of the ET surface. Set to 1 meter below top when None. The default is None.

  • depth (str, float or xr.DataArray, optional) – The ET extinction depth. Set to 1 meter (below surface) when None. The default is None.

  • auxiliary (str or list of str) – name(s) of data arrays to include as auxiliary data to reclist

  • **kwargs (TYPE) – DESCRIPTION.

Raises:

DESCRIPTION.

ValueError

DESCRIPTION.

Returns:

evt – evapotranspiiration package.

Return type:

flopy.mf6.ModflowGwfevt

nlmod.gwf.recharge.ds_to_rch(gwf, ds, mask=None, pname='rch', recharge='recharge', auxiliary=None, **kwargs)[source]

Convert recharge data in the model dataset to a rch package with time series.

Parameters:
  • gwf (flopy.mf6.modflow.mfgwf.ModflowGwf) – groundwater flow model.

  • ds (xr.DataSet) – dataset containing relevant model grid information

  • mask (xr.DataArray) – data array containing mask, recharge is only added where mask is True

  • pname (str, optional) – package name. The default is ‘rch’.

  • recharge (str, xr.DataArray or float, optional) – When recharge is a string, it is the name of the variable in ds that contains the recharge flux rate. When recharge is a float, it is interpreted as the constant recharge rate that is applied in all active cells (within mask if mask is supplied). The default is “recharge”.

  • auxiliary (str or list of str) – name(s) of data arrays to include as auxiliary data to reclist

Returns:

rch – recharge package

Return type:

flopy.mf6.ModflowGwfrch

nlmod.gwf.recharge.ds_to_uzf(gwf, ds, mask=None, pname='uzf', surfdep=0.05, vks='kv', thtr=0.1, thts=0.3, thti=0.1, eps=3.5, landflag=None, finf='recharge', pet='evaporation', extdp=None, extwc=None, ha=None, hroot=None, rootact=None, simulate_et=True, linear_gwet=True, unsat_etwc=False, unsat_etae=False, obs_depth_interval=None, obs_z=None, mask_obs=None, **kwargs)[source]

Create a unsaturated zone flow package for modflow 6.

This method adds uzf-cells to all active Modflow cells (unless mask is specified).

Parameters:
  • gwf (flopy.mf6.modflow.mfgwf.ModflowGwf) – groundwater flow model.

  • ds (xr.DataSet) – dataset containing relevant model grid information

  • mask (xr.DataArray) – data array containing mask, uzf is only added where mask is True

  • pname (str, optional) – package name. The default is ‘uzf’.

  • surfdep (float, str or array-like) – surfdep is the surface depression depth of the UZF cell. When passed as string, the array is obtained from ds. The default is 0.05 m.

  • vks (float, str or array-like) – vks is the saturated vertical hydraulic conductivity of the UZF cell. This value is used with the Brooks-Corey function and the simulated water content to calculate the partially saturated hydraulic conductivity. When passed as string, the array is obtained from ds. The default is ‘kv’.

  • thtr (float, str or array-like) – thtr is the residual (irreducible) water content of the UZF cell. This residual water is not available to plants and will not drain into underlying aquifer cells. When passed as string, the array is obtained from ds. The default is 0.1.

  • thts (float, str or array-like) – thts is the saturated water content of the UZF cell. The values for saturated and residual water content should be set in a manner that is consistent with the specific yield value specified in the Storage Package. The saturated water content must be greater than the residual content. When passed as string, the array is obtained from ds. The default is 0.3.

  • thti (float, str or array-like) – thti is the initial water content of the UZF cell. The value must be greater than or equal to the residual water content and less than or equal to the saturated water content. When passed as string, the array is obtained from ds. The default is 0.1.

  • eps (float, str or array-like) – eps is the exponent used in the Brooks-Corey function. The Brooks-Corey function is used by UZF to calculated hydraulic conductivity under partially saturated conditions as a function of water content and the user-specified saturated hydraulic conductivity. Values must be between 3.5 and 14.0. When passed as string, the array is obtained from ds. The default is 3.5.

  • landflag (xr.DataArray, optional) – A DataArray with integer values, set to one for land surface cells indicating that boundary conditions can be applied and data can be specified in the PERIOD block. A value of 0 specifies a non-land surface cell. Landflag is set to one for the most upper active layer in each vertical column (determined form ds) when it is None. The default is None.

  • finf (float, str or array-like) – The applied infiltration rate of the UZF cell (\(LT^{-1}\)). When passed as string, the array is obtained from ds. The default is “recharge”.

  • pet (float, str or array-like) – The potential evapotranspiration rate of the UZF cell and specified GWF cell. Evapotranspiration is first removed from the unsaturated zone and any remaining potential evapotranspiration is applied to the saturated zone. only used if simulate_et=True. When passed as string, the array is obtained from ds. The default is “evaporation”.

  • extdp (float, optional) – Value that defines the evapotranspiration extinction depth of the UZF cell, in meters below the top of the model. Set to 2.0 meter when None. The default is None.

  • extwc (float, optional) – The evapotranspiration extinction water content of the UZF cell. Only used if simulate_et=True and unsat_etwc=True. Set to thtr when None. The default is None.

  • ha (float, optional) – The air entry potential (head) of the UZF cell. Only used if simulate_et=True and unsat_etae=True. Set to 0.0 when None. The default is None.

  • hroot (float, optional) – The root potential (head) of the UZF cell. Only used if simulate_et=True and unsat_etae=True. Set to 0.0 when None. The default is None.

  • rootact (float, optional) – the root activity function of the UZF cell. ROOTACT is the length of roots in a given volume of soil divided by that volume. Values range from 0 to about 3 \(cm^{-2}\), depending on the plant community and its stage of development. Only used if simulate_et=True and unsat_etae=True. Set to 0.0 when None. The default is None.

  • simulate_et (bool, optional) – If True, ET in the unsaturated (UZF) and saturated zones (GWF) will be simulated. The default is True.

  • linear_gwet (bool, optional) – If True, groundwater ET will be simulated using the original ET formulation of MODFLOW-2005. When False, and square_gwet is True as an extra argument, no evaporation from the saturated zone will be simulated. The default is True.

  • unsat_etwc (bool, optional) – If True, ET in the unsaturated zone will be simulated as a function of the specified PET rate while the water content (THETA) is greater than the ET extinction water content (EXTWC). The default is False.

  • unsat_etae (bool, optional) – If True, ET in the unsaturated zone will be simulated using a capillary pressure based formulation. Capillary pressure is calculated using the Brooks-Corey retention function. The default is False.

  • obs_depth_interval (float, optional) – The depths at which observations of the water content in each cell are added. When not None, this creates a CSV output file with water content at different z-coordinates in each UZF cell. The default is None.

  • obs_z (array-like, optional) – The z-coordinate at which observations of the water content in each cell are added. When not None, this creates a CSV output file with water content at fixes z-coordinates in each UZF cell. The default is None.

  • mask_obs (xr.DataArray, optional) – Mask with the cells where an observations is added. If None all cells will get an observation. The default is None.

  • kwargs (**) – Kwargs are passed onto flopy.mf6.ModflowGwfuzf

Returns:

uzf – Unsaturated zone flow package for Modflow 6.

Return type:

flopy.mf6.ModflowGwfuzf

nlmod.gwf.surface_water._get_start_summer_and_winter(start_summer, start_winter, summer_months)[source]
nlmod.gwf.surface_water._get_waterboard_selection(gdf=None, extent=None, config=None)[source]

Internal method to select waterboards to get data from.

nlmod.gwf.surface_water._is_in_summer(time, start_summer, start_winter)[source]
nlmod.gwf.surface_water.add_bottom_height_from_waterboards(gdf, wc=None, extent=None, columns=None, config=None, min_total_overlap=0.0)[source]

Add information from watercourses to bgt-polygons.

Parameters:
  • gdf (geopandas.GeoDataFrame) – A GeoDataFrame with surface water features, containing the column “bronhouder”.

  • wc (dict, optional) – A dictionary with the name of the waterboards as keys and GeoDataFrames with watercourses as values. It is generated by download_watercourses when None. The default is None.

  • extent (list, tuple or np.array) – Model extent (xmin, xmax, ymin, ymax). When extent is None, all data of the water boards in gdf are downloaded downloaded.

  • columns (TYPE, optional) – The columns that are added to gdf. Columns defaults to ‘bottom_height’ when None. The default is None.

  • config (dict, optional) – A dictionary with information about the webservices of the water boards. When config is None, it is created with nlmod.read.waterboard.get_configuration(). The default is None.

  • min_total_overlap (float, optional) – Only add data from waterboards to gdf when the total overlap between a feature in gdf with all the features from the waterboard is larger than the fraction min_total_overlap. The default is 0.0.

Returns:

gdf – A GeoDataFrame with surface water features, with the added columns

Return type:

geopandas.GeoDataFrame

nlmod.gwf.surface_water.add_info_to_gdf(*args, **kwargs)[source]
nlmod.gwf.surface_water.add_min_ahn_to_gdf(gdf, ahn, buffer=0.0, column='ahn_min', statistic='min', **kwargs)[source]

Add a column with the minimum surface level height near surface water features.

Parameters:
  • gdf (geopandas.GeoDataFrame) – A GeoDataFrame with surface water features

  • ahn (xr.DataArray) – A DataArray containing the height of the surface level.

  • buffer (float, optional) – The buffer that is applied around surface water features to calculate the minimum surface level near these features. The default is 0.0.

  • column (string, optional) – The name of the new column in gdf containing the minimum surface level height. The default is ‘ahn_min’.

  • statistic (string, optional) – The statistic to calculate at each surface water feature. The default is ‘min’.

  • progressbar (bool, optional) – Show a progressbar when True. The default is False.

Returns:

gdf – A GeoDataFrame with surface water features, with an added column containing the minimum surface level height near the features.

Return type:

geopandas.GeoDataFrame

nlmod.gwf.surface_water.add_season_timeseries(ds, package, summer_months=(4, 5, 6, 7, 8, 9), filename='season.ts', summer_name='summer', winter_name='winter', start_summer=None, start_winter=None)[source]

Add time series indicating which season is active (e.g. summer/winter).

Parameters:
  • ds (xr.Dataset) – Xarray dataset used for time discretization.

  • package (flopy.mf6 package) – Modflow 6 package to add time series to.

  • summer_months (tuple, optional) – The summer months (one-based). The parameter summer_months is used to calculate start_summer or start_winter, when they are None. The default is (4, 5, 6, 7, 8, 9), so from april to september.

  • filename (str, optional) – The name of time series file. The default is “season.ts”.

  • winter_name (str, optional) – The name of the time-series with ones in winter. The default is “winter”.

  • summer_name (str, optional) – The name of the time-series with ones in summer. The default is “summer”.

  • start_summer (str, optional) – A string with the month and day of the start of summer (one-based), seperated by ‘-’. For example ‘4-1’ for april 1st. When start_summer is None it is calculated from the parameter ‘summer_months’. The default is None.

  • start_winter (str, optional) – A string with the month and day of the start of winter (one-based), seperated by ‘-’. For example ‘10-1’ for october 1st. When start_winter is None it is calculated from the parameter ‘summer_months’. The default is None.

nlmod.gwf.surface_water.add_stages_from_waterboards(gdf, la=None, extent=None, columns=None, config=None, min_total_overlap=0.0, silent=False)[source]

Add information from level areas (peilgebieden) to bgt-polygons.

Parameters:
  • gdf (geopandas.GeoDataFrame) – A GeoDataFrame with surface water features, containing the column “bronhouder”.

  • la (dict, optional) – A dictionary with the name of the waterboards as keys and GeoDataFrames with level areas as values. It is generated by download_level_areas when None. The default is None.

  • extent (list, tuple or np.array) – Model extent (xmin, xmax, ymin, ymax). When extent is None, all data of the water boards in gdf are downloaded downloaded.

  • columns (TYPE, optional) – The columns that are added to gdf. Columns defaults to ‘summer_stage’ and ‘winter_stage’ when None. The default is None.

  • config (dict, optional) – A dictionary with information about the webservices of the water boards. When config is None, it is created with nlmod.read.waterboard.get_configuration(). The default is None.

  • min_total_overlap (float, optional) – Only add data from waterboards to gdf when the total overlap between a feature in gdf with all the features from the waterboard is larger than the fraction min_total_overlap. The default is 0.0.

  • silent (bool, optional) – If true, do not show prgressbars. The default is False.

Returns:

gdf – A GeoDataFrame with surface water features, with the added columns

Return type:

geopandas.GeoDataFrame

nlmod.gwf.surface_water.agg_area_weighted(gdf, col)[source]
nlmod.gwf.surface_water.agg_de_lange(group, cid, ds, c1=0.0, c0=1.0, N=0.001, crad_positive=True)[source]
nlmod.gwf.surface_water.agg_max_area(gdf, col)[source]
nlmod.gwf.surface_water.aggregate(gdf, method, ds=None, cachedir=None, cachename=None)[source]

Aggregate surface water features.

Parameters:
  • gdf (geopandas.GeoDataFrame) – GeoDataFrame containing surfacewater polygons per grid cell. Must contain columns ‘stage’ (waterlevel), ‘c0’ (bottom resistance), and ‘botm’ (bottom elevation)

  • method (str, optional) – “area_weighted” for area-weighted params, “max_area” for max area params “de_lange” for De Lange formula for conductance

  • ds (xr.Dataset, optional) – DataSet containing model layer information (only required for method=’de_lange’)

  • cachedir (str or None, optional) – directory to save cache. If None no cache is used. Default is None.

  • cachename (str or None, optional) – filename of netcdf cache. If None no cache is used. Default is None.

Returns:

celldata – DataFrame with aggregated surface water parameters per grid cell

Return type:

pd.DataFrame

nlmod.gwf.surface_water.build_spd(celldata, pkg, ds, layer_method='lay_of_rbot', desc=None, silent=False)[source]

Build stress period data for package (RIV, DRN, GHB).

Parameters:
  • celldata (geopandas.GeoDataFrame) – GeoDataFrame containing data. Cellid must be the index, and must have columns “rbot”, “stage” and “cond”. Optional columns are ‘boundname’ and ‘aux’. These columns should have type str.

  • pkg (str) – Modflow package: RIV, DRN or GHB

  • ds (xr.Dataset) – DataSet containing model layer information

  • layer_method (layer_method : str, optional) – The method used to distribute the conductance over the layers. Possible values are ‘lay_of_rbot’ and ‘distribute_cond_over_lays’. The default is “lay_of_rbot”.

  • desc (string, optional) – The description of the progressbar. The default is None, so desc will be “Building stress period data RIV/DRN/GHB”.

  • silent (bool, optional) – Do not show a progressbar when silent is True. The default is False.

Returns:

spd – list containing stress period data: - RIV: [(cellid), stage, cond, rbot] - DRN: [(cellid), elev, cond] - GHB: [(cellid), elev, cond]

Return type:

list

nlmod.gwf.surface_water.coth(x)[source]
nlmod.gwf.surface_water.de_lange_eqns(A, H0, kv, kh, c1, li, Bin, c0, p, N, crad_positive=True)[source]

Calculates the conductance according to De Lange.

Parameters:
  • A (float) – celoppervlak (m2)

  • H0 (float) – doorstroomde dikte (m)

  • kv (float) – verticale doorlotendheid (m/d)

  • kh (float) – horizontale doorlatendheid (m/d)

  • c1 (float) – deklaagweerstand (d)

  • li (float) – lengte van de waterlopen (m)

  • Bin (float) – bodembreedte (m)

  • c0 (float) – slootbodemweerstand (d)

  • p (float) – water peil

  • N (float) – grondwateraanvulling

  • crad_positive (bool, optional) – whether to allow negative crad values. If True, crad will be set to 0 if it is negative.

Returns:

Conductance (m2/d)

Return type:

float

nlmod.gwf.surface_water.distribute_cond_over_lays(cond, cellid, rivbot, laytop, laybot, idomain=None, kh=None, stage=None)[source]

Distribute the conductance in a cell over the layers in that cell, based on the the river-bottom and the layer bottoms, and optionally based on the stage and the hydraulic conductivity.

nlmod.gwf.surface_water.download_level_areas(gdf=None, extent=None, config=None, raise_exceptions=True, drop_duplicates=True, **kwargs)[source]

Download level areas (peilgebieden) of bronhouders.

Parameters:
  • gdf (geopandas.GeoDataFrame, optional) – A GeoDataFrame with surface water features, containing the column “bronhouder”.

  • extent (list, tuple or np.array, optional) – Model extent (xmin, xmax, ymin, ymax). When extent is None, all data of the water boards in gdf are downloaded downloaded.

  • config (dict, optional) – A dictionary with information about the webservices of the water boards. When config is None, it is created with nlmod.read.waterboard.get_configuration(). The default is None.

  • raise_exceptions (bool, optional) – Raises exceptions, mostly caused by a webservice that is offline. When raise_exceptions is False, the error is raised as a warning. The default is True.

  • drop_duplicates (bool, optional) – Drop features with a duplicate index, keeping the first occurence. The default is True.

Returns:

la – A dictionary with the name of the waterboards as keys and GeoDataFrames with level areas as values.

Return type:

dict

nlmod.gwf.surface_water.download_watercourses(gdf=None, extent=None, config=None, raise_exceptions=True, **kwargs)[source]

Download watercourses of bronhouders.

Parameters:
  • gdf (geopandas.GeoDataFrame, optional) – A GeoDataFrame with surface water features, containing the column “bronhouder”. Determine the required waterboards for this gdf, when not None. The default is None.

  • extent (list, tuple or np.array, optional) – Model extent (xmin, xmax, ymin, ymax). When extent is None, all data of the water boards in gdf are downloaded downloaded. The default is None.

  • config (dict, optional) – A dictionary with information about the webservices of the water boards. When config is None, it is created with nlmod.read.waterboard.get_configuration(). The default is None.

  • raise_exceptions (bool, optional) – Raises exceptions, mostly caused by a webservice that is offline. When raise_exceptions is False, the error is raised as a warning. The default is True.

Returns:

wc – A dictionary with the name of the waterboards as keys and GeoDataFrames with watercourses as values.

Return type:

dict

nlmod.gwf.surface_water.estimate_polygon_length(gdf)[source]
nlmod.gwf.surface_water.gdf_to_seasonal_pkg(gdf, gwf, ds, pkg='DRN', default_water_depth=0.5, boundname_column='identificatie', c0=1.0, summer_months=(4, 5, 6, 7, 8, 9), layer_method='lay_of_rbot', silent=False, start_summer=None, start_winter=None, season_filename='season.ts', **kwargs)[source]

Add a surface water package to a groundwater-model, based on input from a GeoDataFrame. This method adds two boundary conditions for each record in the GeoDataFrame: one for the winter_stage and one for the summer_stage. The conductance of each record is a time-series called ‘winter’ or ‘summer’ with values of either 0 or 1. These conductance values are multiplied by an auxiliary variable that contains the actual conductance.

Parameters:
  • gdf (GeoDataFrame) – A GeoDataFrame with Polygon-data. Cellid must be the index and must have columns ‘winter_stage’ and ‘summer_stage’.

  • gwf (flopy ModflowGwf) – groundwaterflow object.

  • ds (xr.Dataset) – Dataset with model data.

  • pkg (str, optional) – The package to generate. Possible options are ‘DRN’, ‘RIV’ and ‘GHB’. The default is ‘DRN’.

  • default_water_depth (float, optional) – The default water depth, only used when there is no ‘rbot’ column in gdf or when this column contains nans. The default is 0.5.

  • boundname_column (str, optional) – The name of the column in gdf to use for the boundnames. The default is “identificatie”, which is a unique identifier in the BGT.

  • c0 (float, optional) – The resistance of the surface water, in days. Only used when there is no ‘cond’ column in gdf. The default is 1.0.

  • summer_months (list or tuple, optional) – The months in which ‘summer_stage’ is active (one-based). The parameter summer_months is used to calculate start_summer or start_winter, when they are None. The default is (4, 5, 6, 7, 8, 9), which means summer is from april through september.

  • layer_method (str, optional) – The method used to distribute the conductance over the layers. Possible values are ‘lay_of_rbot’ and ‘distribute_cond_over_lays’. The default is “lay_of_rbot”.

  • silent (bool, optional) – Do not show a progressbar when silent is True. The default is False.

  • start_summer (str, optional) – A string with the month and day of the start of summer (one-based), seperated by ‘-’. For example ‘4-1’ for april 1st. When start_summer is None it is calculated from the parameter ‘summer_months’. The default is None.

  • start_winter (str, optional) – A string with the month and day of the start of winter (one-based), seperated by ‘-’. For example ‘10-1’ for october 1st. When start_winter is None it is calculated from the parameter ‘summer_months’. The default is None.

  • season_filename (str, optional) – A string with the name of the file containing the timeseries of the seasons. The default is ‘season.ts’.

  • **kwargs (dict) – Kwargs are passed onto ModflowGwfdrn, ModflowGwfriv or ModflowGwfghb.

Returns:

package – The generated flopy-package

Return type:

ModflowGwfdrn, ModflowGwfriv or ModflowGwfghb

nlmod.gwf.surface_water.get_gdf(ds=None, extent=None, fname_ahn=None, ahn=None, buffer=0.0)[source]

Generate a GeoDataFrame based on BGT-data and data from waterboards.

Parameters:
  • ds (TYPE, optional) – The Model Dataset, used to determine the extent (when None) and to grid the surface level features. The default is None.

  • extent (list, tuple or np.array) – Model extent (xmin, xmax, ymin, ymax). When extent is None, extent is extracted from ds

  • fname_ahn (str, optional) – When not None, fname_ahn is the path to a tiff-file with ahn-data, to calculate the minimum height of the surface level near the surface water features. The default is None.

  • ahn (xr.DataArray, optional) – When not None, ahn is a DataArray containing the height of the surface level and is used to calculate the minimum height of the surface level near the surface water features. The default is None.

  • buffer (float, optional) – The buffer that is applied around surface water features to calculate the minimum surface level near these features. The default is 0.0.

Returns:

gdf – A GeoDataFrame with surface water features, with added columns from waterboards and gridded to the model grid (when ds is aupplied)

Return type:

geopandas.GeoDataFrame

nlmod.gwf.surface_water.get_gdf_stage(gdf, season='winter')[source]

Get the stage from a GeoDataFrame for a specific season.

Parameters:
  • gdf (GeoDataFrame) – A GeoDataFrame of the polygons of the BGT with added information in the columns ‘summer_stage’, ‘winter_stage’, and ‘ahn_min’.

  • season (str, optional) – The season for which the stage needs to be determined. The default is “winter”.

Returns:

stage – The stage for each of the records in the GeoDataFrame.

Return type:

pandas.Series

nlmod.gwf.surface_water.get_seaonal_timeseries(ds, summer_value, winter_value, summer_months=(4, 5, 6, 7, 8, 9), start_summer=None, start_winter=None)[source]

Get a time-series of winter-values and summer-values.

The index of the returned series is set at the start of the period that the value describes. This is the way that flopy uses this value in a timeseries, when interpolation_methodrecord is set to ‘stepwise’.

Parameters:
  • ds (xr.Dataset) – Xarray dataset used for time discretization.

  • summer_value (float) – The value to be used in summer.

  • winter_value (float) – The value to be used in winter.

  • summer_months (tuple, optional) – summer months (one-based). The parameter summer_months is used to calculate start_summer or start_winter, when they are None. The default is (4, 5, 6, 7, 8, 9), so from april to september.

  • start_summer (str, optional) – A string with the month and day of the start of summer (one-based), seperated by ‘-’. For example ‘4-1’ for april 1st. When start_summer is None it is calculated from the parameter ‘summer_months’. The default is None.

  • start_winter (str, optional) – A string with the month and day of the start of winter (one-based), seperated by ‘-’. For example ‘10-1’ for october 1st. When start_winter is None it is calculated from the parameter ‘summer_months’. The default is None.

Returns:

s – A series with stages as the values and time as the index.

Return type:

pd.Series

nlmod.gwf.surface_water.get_subsurface_params_by_cellid(ds, cid)[source]
nlmod.gwf.surface_water.get_surfacewater_params(group, method, cid=None, ds=None, delange_params=None)[source]
nlmod.gwf.surface_water.radial_resistance(L, B, H, kh, kv)[source]
nlmod.gwf.surface_water.rivdata_from_xylist(gwf, xylist, layer, stage, cond, rbot, aux=None)[source]
nlmod.gwf.wells._add_cellid(df, ds=None, gwf=None, x='x', y='y', silent=False)[source]

Intersect a DataFrame of point Data with the model grid, and add cellid-column.

Parameters:
  • df (pd.DataFrame or gpd.GeoDataFrame) – A (Geo)DataFrame containing the properties of the wells.

  • ds (xarray.Dataset) – Dataset with model data. Either supply ds or gwf. The default is None.

  • gwf (flopy ModflowGwf) – Groundwaterflow object. Only used when ds is None. The default is None.

  • x (str, optional) – The column in df that contains the x-coordinate of the well. Only used when df is a DataFrame. The default is ‘x’.

  • y (str, optional) – The column in df that contains the y-coordinate of the well. Only used when df is a DataFrame. The default is ‘y’.

  • silent (bool, optional) – Hide progressbar when silent is True. Default is False.

Returns:

df – A GeoDataFrame with a column named cellid that contains the icell2d-number (vertex-grid) or (row, column) (structured grid).

Return type:

gpd.GeoDataFrame

nlmod.gwf.wells._get_layer_multiplier_for_well(cid, well_top, well_bot, ml_top, ml_bot, ml_kh)[source]

Get a factor (numpy array) for each layer that a well screen intersects with.

Parameters:
  • cid (int or tuple of 2 ints) – THe cellid of the well (either icell2d or (row, column).

  • well_top (float) – The z-coordinate of the top of the well screen.

  • well_bot (float) – The z-coordinate of the top of the well screen.

  • ml_top (numpy array) – The top of the upper layer of the model (1d or 2d)

  • ml_bot (numpy array) – The bottom of all cells of the model (2d or 3d)

  • ml_kh (numpy array) – The horizontal conductivity of all cells of the model (2d or 3d).

Returns:

multiplier – An array with a factor (between 0 and 1) for each of the model layers.

Return type:

numpy array

nlmod.gwf.wells._get_layer_multiplier_for_wells(df, top, botm, ds=None, gwf=None)[source]

Get factors (pandas.DataFrame) for each layer that well screens intersects with.

Parameters:
  • df (pd.DataFrame) – A DataFrame containing the properties of the wells.

  • top (str) – The column in df that contains the z-coordinate of the top of the well screen.

  • botm (str) – The column in df that contains the z-coordinate of the bottom of the well screen.

  • ds (xarray.Dataset) – Dataset with model data. Either supply ds or gwf. The default is None.

  • gwf (flopy ModflowGwf) – Groundwaterflow object. Only used when ds is None. The default is None.

Returns:

multipliers – A DataFrame containg the multiplication factors, with the layers as the index and the name of the well screens (the index of df) as columns.

Return type:

pd.DataFrame

nlmod.gwf.wells.maw_from_df(df, gwf, x='x', y='y', top='top', botm='botm', Q='Q', rw='rw', radius_skin=None, hk_skin=None, condeqn='THIEM', strt=None, group=None, aux=None, boundnames=None, ds=None, silent=False, **kwargs)[source]

Add a Multi-Aquifer Well (MAW) package based on input from a (Geo)DataFrame.

Parameters:
  • df (pd.DataFrame or gpd.GeoDataFrame) – A (Geo)DataFrame containing the properties of the wells.

  • gwf (flopy ModflowGwf) – Groundwaterflow object to add the wel-package to.

  • x (str, optional) – The column in df that contains the x-coordinate of the well. Only used when df is a DataFrame. The default is ‘x’.

  • y (str, optional) – The column in df that contains the y-coordinate of the well. Only used when df is a DataFrame. The default is ‘y’.

  • top (str) – The column in df that contains the z-coordinate of the top of the well screen. The defaults is ‘top’.

  • botm (str) – The column in df that contains the z-coordinate of the bottom of the well screen. The defaults is ‘botm’.

  • Q (str, optional) – The column in df that contains the volumetric well rate. This column can contain floats, or strings belonging to timeseries added later. A positive value indicates recharge (injection) and a negative value indicates discharge (extraction). If wells are grouped, the values refer to the rates of all the wells in the group combined and thus must be the same for all wells in the group. The default is “Q”.

  • rw (str, optional) – The column in df that contains the radius for the multi-aquifer well. The default is “rw”.

  • radius_skin (str, optional) – The column in df that contains the radius of the skin around the well; the distance between the center of the well and the outside of the filter pack. Is larger than rw. Only used if condeqn is SKIN, CUMULATIVE, or MEAN. The default is None, which means that the skin is not used.

  • hk_skin (str, optional) – The column in df that contains the horizontal hydraulic conductivity of the skin around the well. Only used if condeqn is SKIN, CUMULATIVE, or MEAN. The defaultis None, which means that the skin is not used.

  • condeqn (str, optional) – String that defines the conductance equation that is used to calculate the saturated conductance for the multi-aquifer well. The default is “THIEM”.

  • strt (float, optional) – The starting head for the multi-aquifer well. The default is None, which uses model surface level as the strt value.

  • group (str, optional) – The column in df that contains the group name for the wells. If this is not None, wells with the same group name are grouped together such that the rate is divided over the wells in the group. Note that empty strings are treated as unique group names, so wells with an empty string in the group column are treated as a separate group. The default is None, which means that each well is treated as a separate well.

  • aux (str of list of str, optional) – The column(s) in df that contain auxiliary variables. The default is None.

  • boundnames (str, optional) – The column in df that contains the boundary names. The default is None.

  • ds (xarray.Dataset) – Dataset with model data. Needed to determine cellid when grid-rotation is used. The default is None.

  • silent (bool, optional) – Hide progressbar when silent is True. Default is False.

  • **kwargs (TYPE) – Kwargs are passed onto ModflowGwfmaw.

Returns:

wel – maw package.

Return type:

flopy.mf6.ModflowGwfmaw

nlmod.gwf.wells.wel_from_df(df, gwf, x='x', y='y', top='top', botm='botm', Q='Q', aux=None, boundnames=None, ds=None, auxmultname='multiplier', **kwargs)[source]

Add a Well (WEL) package based on input from a (Geo)DataFrame.

Parameters:
  • df (pd.DataFrame or gpd.GeoDataFrame) – A (Geo)DataFrame containing the properties of the wells.

  • gwf (flopy ModflowGwf) – Groundwaterflow object to add the wel-package to.

  • x (str, optional) – The column in df that contains the x-coordinate of the well. Only used when df is a DataFrame. The default is ‘x’.

  • y (str, optional) – The column in df that contains the y-coordinate of the well. Only used when df is a DataFrame. The default is ‘y’.

  • top (str) – The column in df that contains the z-coordinate of the top of the well screen. The defaults is ‘top’.

  • botm (str) – The column in df that contains the z-coordinate of the bottom of the well screen. The defaults is ‘botm’.

  • Q (str, optional) – The column in df that contains the volumetric well rate. This column can contain floats, or strings belonging to timeseries added later. A positive value indicates recharge (injection) and a negative value indicates discharge (extraction) The default is “Q”.

  • aux (str of list of str, optional) – The column(s) in df that contain auxiliary variables. The default is None.

  • boundnames (str, optional) – THe column in df thet . The default is None.

  • ds (xarray.Dataset) – Dataset with model data. Needed to determine cellid when grid-rotation is used. The default is None.

  • auxmultname (str, optional) – The name of the auxiliary varibale that contains the multiplication factors to distribute the well discharge over different layers. When auxmultname is None, this auxiliary variable will not be added, and Q is multiplied by these factors directly. auxmultname cannot be None when df[Q] contains strings (the names of timeseries). The default is “multiplier”.

  • **kwargs (dict) – kwargs are passed to flopy.mf6.ModflowGwfwel.

Returns:

wel – wel package.

Return type:

flopy.mf6.ModflowGwfwel

nlmod.gwf.lake._copy_da_from_ds(gdf, ds, variable, boundname_column=None, set_to_0_in_ds=False)[source]
nlmod.gwf.lake._get_and_check_single_value(lake_gdf, column)[source]
nlmod.gwf.lake._parse_laksetting_value(value, ds, key, iper)[source]
nlmod.gwf.lake.add_lakeno_to_gdf(gdf, boundname_column)[source]
nlmod.gwf.lake.clip_meteorological_data_from_ds(gdf, ds, boundname_column=None)[source]

Clip meteorlogical data from the model dataset, and return rainfall and evaporation. This method retrieves the values of rainfall and evaporation from a model Dataset. It uses the ‘recharge’variable, and optionally the ‘evaporation’-variable, and returns a rainfall- and evaporation-DataFrame. These dataframes contain input for each of the lakes. The columns of this DataFrame are either the boundnames (when boundname_column is specified) or the lake-number (lakeno).

clip_meteorological_data_from_ds sets the meteorological data in the model detaset to 0. It turns out MODFLOW 6 already does this, and so this can create problemsm, see https://github.com/gwmod/nlmod/issues/497. Therefore, this method is deprected, and replaced by copy_meteorological_data_from_ds.

Parameters:
  • gdf (gpd.GeoDataframe) – geodataframe with the cellids as the index

  • ds (xr.Dataset) – dataset containing relevant model grid and time information

  • boundname_column (str, optional) – The name of the column in gdf to use for the boundnames. When boundname_column is None, the lake-number (lakeno) is used to determine which rows in gdf belong to each lake, and the columns of rainfall and evaporation are set by the lake-number. If boundname_column is not None, the boundnames are used instead of the lake-number. The default is None.

Returns:

  • rainfall (pd.DataFrame) – The rainfall of each lake (columns) in time (index).

  • evaporation (pd.DataFrame) – The evaporation of each lake (columns) in time (index).

nlmod.gwf.lake.copy_meteorological_data_from_ds(gdf, ds, boundname_column=None, set_to_0_in_ds=False)[source]

Copy meteorlogical data from the model dataset, and return rainfall and evaporation. This method retrieves the values of rainfall and evaporation from a model Dataset. It uses the ‘recharge’variable, and optionally the ‘evaporation’-variable, and returns a rainfall- and evaporation-DataFrame. These dataframes contain input for each of the lakes. The columns of this DataFrame are either the boundnames (when boundname_column is specified) or the lake-number (lakeno).

Parameters:
  • gdf (gpd.GeoDataframe) – geodataframe with the cellids as the index

  • ds (xr.Dataset) – dataset containing relevant model grid and time information

  • boundname_column (str, optional) – The name of the column in gdf to use for the boundnames. When boundname_column is None, the lake-number (lakeno) is used to determine which rows in gdf belong to each lake, and the columns of rainfall and evaporation are set by the lake-number. If boundname_column is not None, the boundnames are used instead of the lake-number. The default is None.

  • set_to_0_in_ds (bool, optional) – If True, sets the meteorological data to 0 in ds, which is not recommended (see issue https://github.com/gwmod/nlmod/issues/497). The default is False.

Returns:

  • rainfall (pd.DataFrame) – The rainfall of each lake (columns) in time (index).

  • evaporation (pd.DataFrame) – The evaporation of each lake (columns) in time (index).

nlmod.gwf.lake.lake_from_gdf(gwf, gdf, ds, rainfall=None, evaporation=None, claktype='VERTICAL', boundname_column='identificatie', obs_type='STAGE', surfdep=0.05, pname='lak', gwt=None, obs_type_gwt='CONCENTRATION', rainfall_concentration=0.0, evaporation_concentration=0.0, **kwargs)[source]

Add a lake from a geodataframe.

Parameters:
  • gwf (flopy.mf6.ModflowGwf) – groundwater flow model.

  • gdf (gpd.GeoDataframe) –

    geodataframe with the cellids as the index and the columns:

    lakeno : with the number of the lake strt : with the starting head of the lake clake : with the bed resistance of the lake optional columns are ‘STATUS’, ‘STAGE’, ‘RAINFALL’, ‘EVAPORATION’, ‘RUNOFF’, ‘INFLOW’, ‘WITHDRAWAL’, ‘AUXILIARY’, ‘RATE’, ‘INVERT’, ‘WIDTH’, ‘SLOPE’, ‘ROUGH’. These columns should contain the name of a dataarray in ds with the dimension time.

    if the lake has any outlets they should be specified in the column

    lakeout : the lake number of the outlet, if this is -1 the water is removed from the model. optinal columns are ‘couttype’, ‘outlet_invert’, ‘outlet_width’, ‘outlet_rough’ and ‘outlet_slope’. These columns should contain a unique value for each outlet.

  • ds (xr.Dataset) – dataset containing relevant model grid and time information

  • rainfall (int, float, str, np.array or pd.DataFrame, optional) – The rainfall to be applied on the lakes. If rainfall is a DataFrame, there should be one column for each lake, with lakeno as the column index, or the boundnames if boundname_column is specified. There should be one row for each stress period. To generate rainfall and evaporation from a model dataset, the method clip_meteorological_data_from_ds can be used. If rainfall is not a DataFrame, all lakes have the same rainfall-values. When rainfall is a pandas Series or a numpy array, a value for rainfall for each stress period is taken from this series/array. When rainfall is a float, it is directly used as the rainfall-value. The default is None.

  • evaporation (int, float, str, np.array or pd.DataFrame, optional) – The evaporation to be applied on the lakes. If evaporation is a DataFrame, there should be one column for each lake, with lakeno as the column index, or the boundnames if boundname_column is specified. There should be one row for each stress period. To generate rainfall and evaporation from a model dataset, the method clip_meteorological_data_from_ds can be used. If evaporation is not a DataFrame, all lakes have the same evaporation-values. When evaporation is a pandas Series or a numpy array, a value for evaporation for each stress period is taken from this series/array. When evaporation is a float, it is directly used as the evaporation-value. The default is None.

  • claktype (str, optional) – defines the lake-GWF connection type. For now only VERTICAL is supported. The default is ‘VERTICAL’.

  • boundname_column (str, optional) – The name of the column in gdf to use for the boundnames. The default is “identificatie”, which is a unique identifier in the BGT.

  • surfdep (float, optional) – Defines the surface depression depth for VERTICAL lake-GWF connections. The default is 0.05.

  • pname (str, optional) – name of the lake package. The default is ‘lak’.

  • obs_type (str or list of str, optional) – observation or list of observations to add to the lake package. The default is ‘STAGE’.

  • obs_type_gwt (str or list of str, optional) – observation or list of observations to add to the lake transport package. The default is ‘CONCENTRATION’.

  • **kwargs – passed to flopy.mf6.ModflowGwflak.

Raises:

NotImplementedError

Returns:

lak

Return type:

flopy lake package

nlmod.gwf.output._calculate_gxg(head_bimonthly: DataArray, below_surfacelevel: bool = False) DataArray[source]
nlmod.gwf.output.calculate_gxg(head: DataArray, below_surfacelevel: bool = False, tolerance: Timedelta = Timedelta('7 days 00:00:00')) DataArray[source]

Calculate GxG groundwater characteristics from head time series.

GLG and GHG (average lowest and average highest groundwater level respectively) are calculated as the average of the three lowest (GLG) or highest (GHG) head values per Dutch hydrological year (april - april), for head values measured at a semi-monthly frequency (14th and 28th of every month). GVG (average spring groundwater level) is calculated as the average of groundwater level on 14th and 28th of March, and 14th of April. Supplied head values are resampled (nearest) to the 14/28 frequency.

Hydrological years without all 24 14/28 dates present are discarded for glg and ghg. Years without the 3 dates for gvg are discarded.

This method is copied from imod-python, and edited so that head-DataArray does not need to contain dimensions ‘x’ and ‘y’, so this method also works for refined grids. The original method can be found in: https://gitlab.com/deltares/imod/imod-python/-/blob/master/imod/evaluate/head.py

Parameters:
  • head (xr.DataArray of floats) – Head relative to sea level, in m, or m below surface level if below_surfacelevel is set to True. Must have dimenstion ‘time’.

  • below_surfacelevel (boolean, optional, default: False.) – False (default) if heads are relative to a datum (e.g. sea level). If True, heads are taken as m below surface level.

  • tolerance (pd.Timedelta, default: 7 days.) – Maximum time window allowed when searching for dates around the 14th and 28th of every month.

Returns:

gxg – Dataset containing glg: average lowest head, ghg: average highest head, gvg: average spring head, n_years_gvg: numbers of years used for gvg, n_years_gxg: numbers of years used for glg and ghg.

Return type:

xr.Dataset

Examples

Load the heads, and calculate groundwater characteristics for the simulation period:

>>> import nlmod
>>> head = nlmod.gwf.get_heads_da(ds)
>>> gxg = nlmod.gwf.output.calculate_gxg(head)
nlmod.gwf.output.get_budget_da(text, ds=None, gwf=None, fname=None, grb_file=None, column='q', delayed=False, chunked=False, precision='auto', **kwargs)[source]

Read binary budget file.

Parameters:
  • text (str) – record to get from budget file

  • ds (xarray.Dataset, optional) – xarray dataset with model data. One of ds or gwf must be provided.

  • gwf (flopy ModflowGwf, optional) – Flopy groundwaterflow object. One of ds or gwf must be provided.

  • fname (path, optional) – specify the budget file to load, if not provided budget file will be obtained from ds or gwf.

  • grb_file (str) – path to file containing binary grid information, only needed if reading output from file using fname

  • column (str) – name of column in rec-array to read, default is ‘q’ which contains the fluxes for most budget datasets.

  • delayed (bool, optional) – if delayed is True, do not load output data into memory, default is False.

  • chunked (bool, optional) – chunk data array containing output, default is False.

  • precision (str, optional) – precision of floating point data in the budget-file. Accepted values are ‘auto’, ‘single’ or ‘double’. When precision is ‘auto’, it is determined from the budget-file. Default is ‘auto’.

Returns:

da – budget data array.

Return type:

xarray.DataArray

nlmod.gwf.output.get_cellbudgetfile(ds=None, gwf=None, fname=None, grb_file=None, **kwargs)[source]

Get flopy CellBudgetFile object.

Provide one of ds, gwf or fname.

Parameters:
  • ds (xarray.Dataset, optional) – model dataset, by default None

  • gwf (flopy.mf6.ModflowGwf, optional) – groundwater flow model, by default None

  • fname_cbc (str, optional) – path to cell budget file, by default None

  • grb_file (str, optional) – path to file containing binary grid information, only needed if fname_cbc is passed as only argument.

Returns:

CellBudgetFile object handle

Return type:

flopy.utils.CellBudgetFile

nlmod.gwf.output.get_flow_lower_face(ds, gwf=None, fname=None, grb_file=None, kstpkper=None, lays=None)[source]

Get the flow over the lower face of all model cells.

The flow Lower Face (flf) used to be written to the budget file in previous versions of MODFLOW. In MODFLOW 6 we determine these flows from the flow-ja-face-records.

Parameters:
  • ds (xarray.Dataset) – Xarray dataset with model data.

  • gwf (flopy ModflowGwf, optional) – Flopy groundwaterflow object. One of ds or gwf must be provided.

  • fname (path, optional) – specify the budget file to load, if not provided budget file will be obtained from ds or gwf.

  • grb_file (str, optional) – The location of the grb-file. grb_file is determied from ds when None. The default is None.

  • kstpkper (tuple of 2 ints, optional) – The index of the timestep and the stress period to include in the result. Include all data in the budget-file when None. The default is None.

  • lays (int or list of ints, optional) – The layers to include in the result. When lays is None, all layers are included. The default is None.

Returns:

da – The flow over the lower face of each cell, in m3/d.

Return type:

xr.DataArray

nlmod.gwf.output.get_flow_residuals(ds, gwf=None, fname=None, grb_file=None, kstpkper=None)[source]

Get the flow residuals of a MODFLOW 6 simulation.

Parameters:
  • ds (xarray.Dataset) – Xarray dataset with model data.

  • gwf (flopy ModflowGwf, optional) – Flopy groundwaterflow object. One of ds or gwf must be provided.

  • fname (path, optional) – specify the budget file to load, if not provided budget file will be obtained from ds or gwf.

  • grb_file (str) – The location of the grb-file. grb_file is determied from ds when None. The default is None.

  • kstpkper (tuple of 2 ints, optional) – The index of the timestep and the stress period to include in the result. Include all data in the budget-file when None. The default is None.

Returns:

da – The flow residual in each cell, in m3/d.

Return type:

xr.DataArray

nlmod.gwf.output.get_gwl_from_wet_cells(head, layer='layer', botm=None)[source]

Get the groundwater level from a multi-dimensional head array where dry cells are NaN. This methods finds the most upper non-nan-value of each cell or timestep.

Parameters:
  • head (xarray.DataArray or numpy array) – A multi-dimensional array of head values. NaN-values represent inactive or dry cells.

  • layer (string or int, optional) – The name of the layer dimension of head (if head is a DataArray) or the integer of the layer dimension of head (if head is a numpy array). The default is ‘layer’.

  • botm (xarray.DataArray, optional) – A DataArray with the botm of each model-cell. It can be used to set heads below the botm of the cells to NaN. botm is only used when head is a DataArray.

Returns:

gwl – An array of the groundwater-level, without the layer-dimension.

Return type:

numpy-array

nlmod.gwf.output.get_head_at_point(head, x, y, ds=None, gi=None, drop_nan_layers=True, rotated=True)[source]

Get the head at a certain point from a head DataArray for all cells.

Parameters:
  • head (xarray.DataArray) – A DataArray of heads, with dimensions (time, layer, y, x) or (time, layer, icell2d).

  • x (float) – The x-coordinate of the requested head.

  • y (float) – The y-coordinate of the requested head.

  • ds (xarray.Dataset, optional) – Xarray dataset with model data. Only used when a Vertex grid is used. The default is None.

  • gi (flopy.utils.GridIntersect, optional) – A GridIntersect class, to determine the cell at point x,y. Only used when a Vertex grid is used, and it is determined from ds when None. The default is None.

  • drop_nan_layers (bool, optional) – Drop layers that are NaN at all timesteps. The default is True.

  • rotated (bool, optional) – If the model grid has a rotation, and rotated is False, x and y are in model coordinates. Otherwise x and y are in real world coordinates. The default is True.

Returns:

head_point – A DataArray with dimensions (time, layer).

Return type:

xarray.DataArray

nlmod.gwf.output.get_headfile(ds=None, gwf=None, fname=None, grb_file=None, **kwargs)[source]

Get flopy HeadFile object.

Provide one of ds, gwf or fname.

Parameters:
  • ds (xarray.Dataset, optional) – model dataset, by default None

  • gwf (flopy.mf6.ModflowGwf, optional) – groundwater flow model, by default None

  • fname (str, optional) – path to heads file, by default None

  • grb_file (str) – path to file containing binary grid information

Returns:

HeadFile object handle

Return type:

flopy.utils.HeadFile

nlmod.gwf.output.get_heads_da(ds=None, gwf=None, fname=None, grb_file=None, delayed=False, chunked=False, precision='auto', **kwargs)[source]

Read binary heads file.

Parameters:
  • ds (xarray.Dataset) – Xarray dataset with model data.

  • gwf (flopy ModflowGwf) – Flopy groundwaterflow object.

  • fname (path, optional) – path to a binary heads file

  • grb_file (str, optional) – path to file containing binary grid information, only needed if reading output from file using fname

  • delayed (bool, optional) – if delayed is True, do not load output data into memory, default is False.

  • chunked (bool, optional) – chunk data array containing output, default is False.

  • precision (str, optional) – precision of floating point data in the head-file. Accepted values are ‘auto’, ‘single’ or ‘double’. When precision is ‘auto’, it is determined from the head-file. Default is ‘auto’.

Returns:

da – heads data array.

Return type:

xarray.DataArray

gwt

nlmod.gwt.gwt.adv(ds, gwt, scheme=None, **kwargs)[source]

Create advection package for groundwater transport model.

Parameters:
  • ds (xarray.Dataset) – dataset with model data.

  • gwt (flopy ModflowGwt) – groundwater transport object.

  • scheme (str, optional) – advection scheme to use, by default “UPSTREAM”, options are (“UPSTREAM”, “CENTRAL”, “TVD”).

Returns:

adv – adv package

Return type:

flopy ModflowGwtadv

nlmod.gwt.gwt.cnc(ds, gwt, da_mask, da_conc, pname='cnc', **kwargs)[source]

Create constant concentration package for groundwater transport model.

Parameters:
  • ds (xarray.Dataset) – dataset with model data

  • gwt (flopy ModflowGwt) – groundwater transport object

  • da_mask (str) – data array containing mask where to create constant concentration cells

  • da_conc (str) – data array containing concentration data

Returns:

cnc – cnc package

Return type:

flopy ModflowGwtcnc

nlmod.gwt.gwt.dis(ds, gwt, length_units='METERS', pname='dis', **kwargs)[source]

Create discretisation package from the model dataset.

Parameters:
  • ds (xarray.Dataset) – dataset with model data.

  • gwt (flopy ModflowGwf) – groundwater transport object

  • length_units (str, optional) – length unit. The default is ‘METERS’.

  • pname (str, optional) – package name, ignored if ds has a vertex grid (disv)

Returns:

dis – discretisation package.

Return type:

flopy ModflowGwtdis

nlmod.gwt.gwt.disv(ds, gwt, length_units='METERS', pname='disv', **kwargs)[source]

Create discretisation vertices package from the model dataset.

Parameters:
  • ds (xarray.Dataset) – dataset with model data.

  • gwt (flopy ModflowGwt) – groundwater transport object.

  • length_units (str, optional) – length unit. The default is ‘METERS’.

  • pname (str, optional) – package name

Returns:

disv – disv package

Return type:

flopy ModflowGwtdisv

nlmod.gwt.gwt.dsp(ds, gwt, **kwargs)[source]

Create dispersion package for groundwater transport model.

Parameters:
  • ds (xarray.Dataset) – dataset with model data

  • gwt (flopy ModflowGwt) – groundwater transport object

Returns:

dsp – dsp package

Return type:

flopy ModflowGwtdsp

nlmod.gwt.gwt.gwfgwt(ds, sim, exgtype='GWF6-GWT6', **kwargs)[source]

Create GWF-GWT exchange package for mf6 simulation.

Parameters:
  • ds (xarray.Dataset) – dataset with model data.

  • sim (flopy MFSimulation) – simulation object

  • exgtype (str, optional) – exchange type, by default “GWF6-GWT6”

Returns:

_description_

Return type:

gwfgwt

nlmod.gwt.gwt.gwt(ds, sim, modelname=None, **kwargs)[source]

Create groundwater transport model from the model dataset.

Parameters:
  • ds (xarray.Dataset) – dataset with model data. Should have the dimension ‘time’ and the attributes: model_name, mfversion, model_ws, time_units, start, perlen, nstp, tsmult

  • sim (flopy MFSimulation) – simulation object.

  • modelname (str) – name of the transport model

Returns:

gwt – groundwater transport object.

Return type:

flopy ModflowGwt

nlmod.gwt.gwt.ic(ds, gwt, strt, pname='ic', **kwargs)[source]

Create initial condictions package for groundwater transport model.

Parameters:
  • ds (xarray.Dataset) – dataset with model data.

  • gwt (flopy ModflowGwf) – groundwater transport object.

  • strt (str, float or int) – if type is int or float this is the starting concentration for all cells If the type is str the data variable from ds is used as starting concentration.

  • pname (str, optional) – package name

Returns:

ic – ic package

Return type:

flopy ModflowGwtic

nlmod.gwt.gwt.mst(ds, gwt, porosity=None, **kwargs)[source]

Create mass storage transfer package for groundwater transport model.

Parameters:
  • ds (xarray.Dataset) – dataset with model data

  • gwt (flopy ModflowGwt) – groundwater transport object

  • porosity (Any, optional) – porosity, can be passed as float, array-like or string. If passed as string data is taken from dataset.

Returns:

mst – mst package

Return type:

flopy ModflowGwtmst

nlmod.gwt.gwt.oc(ds, gwt, save_concentration=True, save_budget=True, print_concentration=False, print_budget=False, pname='oc', **kwargs)[source]

Create output control package for groundwater transport model.

Parameters:
  • ds (xarray.Dataset) – dataset with model data.

  • gwt (flopy ModflowGwt) – groundwater transport object.

  • pname (str, optional) – package name

Returns:

oc – oc package

Return type:

flopy ModflowGwtoc

nlmod.gwt.gwt.ssm(ds, gwt, sources=None, **kwargs)[source]

Create source-sink mixing package for groundwater transport model.

Parameters:
  • ds (xarray.Dataset) – dataset with model data

  • gwt (flopy ModflowGwt) – groundwater transport object

  • sources (list of tuple, None) – list of tuple(s) with packages that function as source in model, e.g. [(“GHB”, “AUX”, “CONCENTRATION”)]. If None, sources is derived from model dataset attribute ds.ssm_sources.

Returns:

ssm – ssm package

Return type:

flopy ModflowGwtssm

nlmod.gwt.prepare.set_default_transport_parameters(ds, transport_type)[source]

Set default transport parameters based on type of transport model.

Convenience function for setting several variables at once for which default values are often used.

Parameters:
  • ds (xarray.Dataset) – dataset

  • transport_type (str) – type of transport model, currently supports “chloride” or “tracer”.

Returns:

ds – dataset with transport parameters added to attributes.

Return type:

xarray.Dataset

nlmod.gwt.output.freshwater_head(ds, hp, conc, denseref=None, drhodc=None)[source]

Calculate equivalent freshwater head from point water heads. Heads file produced by mf6 contains point water heads.

Parameters:
  • ds (xarray.Dataset) – model dataset containing layer elevation/thickness data, and reference density (denseref) relationship between concentration and density (drhodc) if not provided separately

  • hp (xarray.DataArray) – data array containing point water heads

  • conc (xarray.DataArray) – data array containing concentration

  • denseref (float, optional) – reference density, by default None, which will use denseref attribute in model dataset.

  • drhodc (float, optional) – density-concentration gradient, by default None, which will use drhodc attribute in model dataset.

Returns:

hf – data array containing equivalent freshwater heads.

Return type:

xarray.DataArray

nlmod.gwt.output.get_concentration_at_gw_surface(conc, layer='layer')[source]

Get the concentration level from a multi-dimensional concentration array where dry or inactive cells are NaN. This methods finds the most upper non- nan-value of each cell or timestep.

Parameters:
  • conc (xarray.DataArray or numpy array) – A multi-dimensional array of conc values. NaN-values represent inactive or dry cells.

  • layer (string or int, optional) – The name of the layer dimension of conc (if conc is a DataArray) or the integer of the layer dimension of conc (if conc is a numpy array). The default is ‘layer’.

Returns:

ctop – an array of the top level concentration, without the layer-dimension.

Return type:

numpy-array or xr.DataArray

nlmod.gwt.output.get_concentration_da(ds=None, gwt=None, fname=None, grbfile=None, delayed=False, chunked=False, **kwargs)[source]

Reads binary concentration file.

Parameters:
  • ds (xarray.Dataset) – Xarray dataset with model data.

  • gwt (flopy ModflowGwt) – Flopy groundwater transport object.

  • fname (path, optional) – Instead of loading the binary concentration file corresponding to ds or gwf load the concentration from this file.

  • grbfile (str) – path to file containing binary grid information, only needed if reading output from file using fname

  • delayed (bool, optional) – if delayed is True, do not load output data into memory, default is False.

  • chunked (bool, optional) – chunk data array containing output, default is False.

Returns:

da – concentration data array.

Return type:

xarray.DataArray

nlmod.gwt.output.get_concentration_obj(ds=None, gwt=None, fname=None, grbfile=None)[source]

Get flopy HeadFile object connected to the file with the concetration of cells.

Provide one of ds, gwf or fname.

Parameters:
  • ds (xarray.Dataset, optional) – model dataset, by default None

  • gwt (flopy.mf6.ModflowGwt, optional) – groundwater transport model, by default None

  • fname (str, optional) – path to heads file, by default None

  • grbfile (str) – path to file containing binary grid information

Returns:

HeadFile object handle

Return type:

flopy.utils.HeadFile

nlmod.gwt.output.pointwater_head(ds, hf, conc, denseref=None, drhodc=None)[source]

Calculate point water head from freshwater heads. Heads file produced by mf6 contains point water heads.

Parameters:
  • ds (xarray.Dataset) – model dataset containing layer elevation/thickness data, and reference density (denseref) relationship between concentration and density (drhodc) if not provided separately

  • hf (xarray.DataArray) – data array containing freshwater heads

  • conc (xarray.DataArray) – data array containing concentration

  • denseref (float, optional) – reference density, by default None, which will use denseref attribute in model dataset.

  • drhodc (float, optional) – density-concentration gradient, by default None, which will use drhodc attribute in model dataset.

Returns:

hf – data array containing point water heads.

Return type:

xarray.DataArray

plot

gis

nlmod.gis._break_down_dimension(ds, variables, dim, add_dim_name=False, add_one_based_index=False)[source]

Internal method to split a dimension of a variable into multiple variables.

Copied and altered from imod-python.

nlmod.gis.dataarray_to_shapefile(model_ds, data_variables, fname, polygons=None)[source]

Save one or more DataArrays from a model dataset as a shapefile.

Parameters:
  • model_ds (xr.DataSet) – xarray with model data

  • data_variables (list, tuple or set) – data_variables in model_ds that will be stored as attributes in the shapefile.

  • fname (str) – filename of the shapefile.

  • polygons (list of shapely Polygons, optional) – geometries used for the GeoDataframe, if None the polygons are created from the data variable ‘vertices’ in model_ds. The default is None.

Return type:

None.

nlmod.gis.ds_to_ugrid_nc_file(model_ds, fname=None, variables=None, dummy_var='mesh_topology', xv='xv', yv='yv', face_node_connectivity='icvert', split_layer_dimension=True, split_time_dimension=False, for_imod_qgis_plugin=False)[source]

Save a model dataset to a UGRID NetCDF file, so it can be opened as a Mesh Layer in qgis.

Parameters:
  • model_ds (xr.DataSet) – xarray Dataset with model data

  • fname (str, optional) – filename of the UGRID NetCDF-file, preferably with the extension .nc. When fname is None,only a ugird-ready Dataset is created, without saving this Dataset to file. The defaults is None.

  • variables (str or list of str, optional) – THe variables to be saved in the NetCDF file. The default is None, which means all variables will be saved in the file.

  • dummy_var (str, optional) – The name of the new dummy-variable that contains the gridinformation in its attributes. The default is ‘mesh_topology’.

  • xv (str, optional) – The name of the variable that contains the x-coordinate of the vertices that together form the edges of the faces. The default is ‘xv’.

  • yv (str, optional) – The name of the variable that contains the y-coordinate of the vertices that together form the edges of the faces. The default is ‘yv’.

  • face_node_connectivity (str, optional) – The name of the variable that contains the indexes of the vertices for each face. The default is ‘icvert’.

  • split_layer_dimension (bool, optional) – Splits the layer dimension into seperate variables when True. The defaults is True.

  • split_time_dimension (bool, optional) – Splits the time dimension into seperate variables when True. The defaults is False.

  • for_imod_qgis_plugin (bool, optional) – When True, set some properties of the netcdf file to improve compatibility with the iMOD-QGIS plugin. Layers are renamed to ‘layer_i’ until ‘layer_n’, a variable ‘top’ is added for each layer, and the variable ‘botm’ is renamed to ‘bottom’. The default is False.

Returns:

ds – The dataset that was saved to a NetCDF-file. Can be used for debugging.

Return type:

xr.DataSet

nlmod.gis.ds_to_vector_file(model_ds, gisdir=None, driver='GPKG', combine_dic=None, exclude=('x', 'y', 'time_steps', 'area', 'vertices', 'rch_name', 'icvert'), crs='+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725 +units=m +no_defs')[source]

Save all data variables in a model dataset to multiple shapefiles.

Parameters:
  • model_ds (xr.DataSet) – xarray with model data

  • gisdir (str, optional) – gis directory to save shapefiles, if None a subdirectory ‘gis’ of the model_ws (which is an attribute of model_ds) is used. The default is None.

  • driver (str, optional) – determines if the data variables are saved as seperate “ESRI Shapefile” or a single “GPKG” (geopackage). The default is geopackage.

  • combine_dic (dictionary of str:set(), optional) – The items in this dictionary are data variables in model_ds that should be combined in one shapefile. The key defines the name of the shapefile. If None the default combine_dic is used. The default is None.

  • exclude (tuple of str, optional) – data variables that are not exported to shapefiles. The default is (‘x’, ‘y’, ‘time_steps’, ‘area’, ‘vertices’).

  • crs (str, optional) – coordinate reference system for the vector file. The default is EPSG:28992 (RD).

Returns:

fnames – filename(s) of exported geopackage or shapefiles.

Return type:

str or list of str

nlmod.gis.struc_da_to_gdf(model_ds, data_variables, polygons=None, dealing_with_time='mean', crs='+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725 +units=m +no_defs')[source]

Convert one or more DataArrays from a structured model dataset to a Geodataframe.

Parameters:
  • model_ds (xr.DataSet) – xarray with model data

  • data_variables (list, tuple or set) – data_variables in model_ds that will be stored as attributes in the shapefile.

  • polygons (list of shapely Polygons, optional) – geometries used for the GeoDataframe, if None the polygons are created from the data variable ‘vertices’ in model_ds. The default is None.

  • crs (str, optional) – coordinate reference system for the geodataframe. The default is EPSG:28992 (RD).

Raises:
  • ValueError – for DataArrays with unexpected dimension names.

  • NotImplementedError – for DataArrays with more than 2 dimensions or dealing_with_time!=’mean’

Returns:

gdf – geodataframe of one or more DataArrays.

Return type:

geopandas.GeoDataframe

nlmod.gis.vertex_da_to_gdf(model_ds, data_variables, polygons=None, dealing_with_time='mean', crs='+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725 +units=m +no_defs')[source]

Convert one or more DataArrays from a vertex model dataset to a Geodataframe.

Parameters:
  • model_ds (xr.DataSet) – xarray with model data

  • data_variables (list, tuple or set) – data_variables in model_ds that will be stored as attributes in the shapefile.

  • polygons (list of shapely Polygons, optional) – geometries used for the GeoDataframe, if None the polygons are created from the data variable ‘vertices’ in model_ds. The default is None.

  • dealing_with_time (str, optional) – when there is time variant data in the model dataset this function becomes very slow. For now only the time averaged data will be saved in the geodataframe. Later this can be extended with multiple possibilities. The default is ‘mean’.

  • crs (str, optional) – coordinate reference system for the geodataframe. The default is EPSG:28992 (RD).

Raises:
  • ValueError – for DataArrays with unexpected dimension names.

  • NotImplementedError – for DataArrays with more than 2 dimensions or dealing_with_time!=’mean’

Returns:

gdf – geodataframe of one or more DataArrays.

Return type:

geopandas.GeoDataframe

util

class nlmod.util.ColoredFormatter(*args, colors: Dict[str, str] | None = None, **kwargs)[source]

Colored log formatter.

Taken from https://gist.github.com/joshbode/58fac7ababc700f51e2a9ecdebe563ad

format(record) str[source]

Format the specified record as text.

exception nlmod.util.LayerError[source]

Generic error when modifying layers.

exception nlmod.util.MissingValueError[source]

Generic error when an expected value is not defined.

nlmod.util._get_value_from_ds_attr(ds, varname, attr=None, value=None, default=None, warn=True)[source]

Internal function to get value from dataset attributes.

Parameters:
  • ds (xarray.Dataset) – dataset containing model data

  • varname (str) – name of the variable in flopy package

  • attr (str, optional) – name of the attribute in dataset (is sometimes different to varname)

  • value (Any, optional) – variable value, by default None

  • default (Any, optional) – When default is not None, value is None, and attr is not present in ds, this default is returned. The default is None.

  • warn (bool, optional) – log warning if value not found

Returns:

value – returns variable value, if value was None, attempts to obtain variable from dataset attributes.

Return type:

Any

nlmod.util._get_value_from_ds_datavar(ds, varname, datavar=None, default=None, warn=True, return_da=False, checkcoords=True)[source]

Internal function to get value from dataset data variables.

Parameters:
  • ds (xarray.Dataset) – dataset containing model data

  • varname (str) – name of the variable in flopy package (used for logging)

  • datavar (Any, optional) – if str, treated as the name of the data variable (which can be different to varname) in dataset, if not provided is assumed to be the same as varname. If not passed as string, it is treated as data

  • default (Any, optional) – When default is not None, datavar is a string, and not present in ds, this default is returned. The default is None.

  • warn (bool, optional) – log warning if value not found

  • return_da (bool, optional) – if True, a DataArray can be returned. If False, a DataArray is always converted to a numpy array before being returned. The default is False.

  • checkcoords (bool, optional) – if True and datavar is a DataArray, the DataArray coords are checked against the Dataset coordinates. Raises an AssertionError if they do not match.

Returns:

value – returns variable value, if value is None or str, attempts to obtain variable from dataset data variables.

Return type:

Any

Note

For optional data, use warn=False, e.g.:

_get_value_from_ds_datavar(ds, "ss", datavar=None, warn=False)
nlmod.util.check_da_dims_coords(da, ds)[source]

Check if DataArray dimensions and coordinates match those in Dataset.

Only checks overlapping dimensions.

Parameters:
  • da (xarray.DataArray) – dataarray to check

  • ds (xarray.Dataset or xarray.DataArray) – dataset or dataarray to compare coords with

Returns:

if dimensions and coordinates match

Return type:

True

Raises:

AssertionError – error that coordinates do not match

nlmod.util.check_presence_mfbinaries(exe_name='mf6', binpath=None)[source]

Check if exe_name is present in the binpath folder.

Parameters:
  • exe_name (str, optional) – the name of the file that is checked to be present, by default ‘mf6’

  • binpath (str, optional) – path to directory to download binaries to, if it doesnt exist it is created. Default is None which sets dir to nlmod/bin.

nlmod.util.compare_model_extents(extent1, extent2)[source]

Check overlap between two model extents.

Parameters:
  • extent1 (list, tuple or array) – first extent [xmin, xmax, ymin, ymax]

  • extent2 (xr.DataSet) – second extent

Returns:

several outcomes:

1: extent1 is completely within extent2 2: extent2 is completely within extent1

Return type:

int

nlmod.util.download_file_from_google_drive(fid, destination=None)[source]

Download a file from google drive using it’s id.

Parameters:
  • fid (str) – google drive id name of a file.

  • destination (str, optional) – location to save the file to. If destination is None the file is written to the current working directory. The default is None.

nlmod.util.download_mfbinaries(bindir=None, version_tag='latest', repo='executables')[source]

Download and unpack platform-specific modflow binaries.

Source: USGS

Parameters:
  • binpath (str, optional) – path to directory to download binaries to, if it doesnt exist it is created. Default is None which sets dir to nlmod/bin.

  • repo (str, default "executables") – Name of GitHub repository. Choose one of “executables” (default), “modflow6”, or “modflow6-nightly-build”.

  • version_tag (str, default "latest") – GitHub release ID.

nlmod.util.download_modpath_provisional_exe(bindir=None, timeout=120)[source]

Download the provisional version of modpath to the folder with binaries.

nlmod.util.extent_to_gdf(extent, crs='EPSG:28992')[source]

Create a geodataframe with a single polygon with the extent given.

Parameters:
  • extent (tuple, list or array) – extent.

  • crs (str, optional) – coördinate reference system of the extent, default is EPSG:28992 (RD new)

Returns:

gdf_extent – geodataframe with extent.

Return type:

geopandas.GeoDataFrame

nlmod.util.extent_to_polygon(extent)[source]

Generate a shapely Polygon from an extent ([xmin, xmax, ymin, ymax])

Parameters:

extent (tuple, list or array) – extent (xmin, xmax, ymin, ymax).

Returns:

polygon of the extent.

Return type:

shapely.geometry.Polygon

nlmod.util.find_most_recent_file(folder, name, extension='.pklz')[source]

Find the most recent file in a folder.

File must startwith name and end width extension. If you want to look for the most recent folder use extension = ‘’.

Parameters:
  • folder (str) – path of folder to look for files

  • name (str) – find only files that start with this name

  • extension (str) – find only files with this extension

Returns:

newest_file – name of the most recent file

Return type:

str

nlmod.util.gdf_from_extent(extent, crs='EPSG:28992')[source]

Create a geodataframe with a single polygon with the extent given.

Parameters:
  • extent (tuple, list or array) – extent.

  • crs (str, optional) – coördinate reference system of the extent, default is EPSG:28992 (RD new)

Returns:

gdf_extent – geodataframe with extent.

Return type:

GeoDataFrame

nlmod.util.gdf_intersection_join(gdf_from, gdf_to, columns=None, desc='', silent=False, min_total_overlap=0.5, geom_type='Polygon', add_index_from_column=None)[source]

Add information from ‘gdf_from’ to ‘gdf_to’, based on the spatial intersection.

Parameters:
  • gdf_from (gpd.GeoDataFrame) – The GeoDataFrame containing the information to be added.

  • gdf_to (gpd.GeoDataFrame) – The GeoDataFrame to add information to.

  • columns (list, optional) – A list of the columns to add from gdf_from to gdf_to. When columns is None, columns is set to all columns that are present in gdf_from but not in gdf_to. The default is None.

  • desc (string, optional) – The description of the progressbar. The default is “”.

  • silent (bool, optional) – If true, do not show the prgressbar. The default is False.

  • min_total_overlap (float, optional) – The minimum required total overlap between geometries. If the total overlap of a feature in gdf_to with the features of gdf_from, no information is added to the feature of gdf_to. The default is 0.5.

  • geom_type (string, optional) – The type of Geometries to evaluate. Can be “Polygon” or “LineString”. When geom_type is “Polygon” the overlap of features are ranked by the area of features. When geom_type is “LineString”, the overlap is ranked by the length of features. The default is “Polygon”.

  • add_index_from_column (string, optional) – The name of the column to add the index of gdf_from to gdf_to. The default is None.

Raises:

TypeError – DESCRIPTION.

Returns:

gdf_to – gdf_to with extra columns from gdf_from.

Return type:

gpd.GeoDataFrame

nlmod.util.gdf_within_extent(gdf, extent)[source]

Select only parts of the geodataframe within the extent.

Only accepts Polygon and Linestring geometry types.

Parameters:
  • gdf (geopandas GeoDataFrame) – dataframe with polygon features.

  • extent (list or tuple) – extent to slice gdf, (xmin, xmax, ymin, ymax).

Returns:

gdf – dataframe with only polygon features within the extent.

Return type:

geopandas GeoDataFrame

nlmod.util.get_bin_directory(exe_name='mf6', bindir=None, download_if_not_found=True, version_tag=None, repo='executables') Path[source]

Get the directory where the executables are stored.

Searching for the executables is done in the following order: 0. If exe_name is a full path, return the full path of the executable. 1. The directory specified with bindir. Raises error if exe_name is provided

and not found. Requires enable_version_check to be False.

  1. The directory used by nlmod installed in this environment.

  2. If the executables were downloaded with flopy/nlmod from an other env,

    most recent installation location of MODFLOW is found in flopy metadata

Else: 4. Download the executables using version_tag and repo.

The returned directory is checked to contain exe_name if exe_name is provided. If exe_name is set to None only the existence of the directory is checked.

Parameters:
  • exe_name (str, optional) – The name of the executable, by default mf6.

  • bindir (Path, optional) – The directory where the executables are stored, by default “None”.

  • download_if_not_found (bool, optional) – Download the executables if they are not found, by default True.

  • repo (str, default "executables") – Name of GitHub repository. Choose one of “executables” (default), “modflow6”, or “modflow6-nightly-build”. If repo and version_tag are provided the most recent installation location of MODFLOW is found in flopy metadata that respects version_tag and repo. If not found, the executables are downloaded using repo and version_tag.

  • version_tag (str, default None) – GitHub release ID: for example “18.0” or “latest”. If repo and version_tag are provided the most recent installation location of MODFLOW is found in flopy metadata that respects version_tag and repo. If not found, the executables are downloaded using repo and version_tag.

Returns:

The directory where the executables are stored.

Return type:

Path

Raises:

FileNotFoundError – If the executables are not found in the specified directories.

nlmod.util.get_color_logger(level='INFO', logger_name=None)[source]

Get a logger with colored output.

Parameters:

level (str, optional) – The logging level to set for the logger. Default is “INFO”.

Returns:

logger – The configured logger object.

Return type:

logging.Logger

nlmod.util.get_da_from_da_ds(da_ds, dims=('y', 'x'), data=None)[source]

Get a dataarray from ds with certain dimensions.

Parameters:
  • da_ds (xr.Dataset or xr.DataArray) – Dataset or DataArray with at least the dimensions defined in dims

  • dims (tuple of string, optional) – dimensions of the desired data array. Default is ‘y’, ‘x’

  • data (None, int, float, array_like, optional) – data to fill data array with. if None DataArray is filled with nans. Default is None.

Returns:

da – DataArray with coordinates from ds

Return type:

xr.DataArray

nlmod.util.get_ds_empty(ds, keep_coords=None)[source]

Get a copy of a dataset with only coordinate information.

Parameters:
  • ds (xr.Dataset) – dataset with coordinates

  • keep_coords (tuple or None, optional) – the coordinates in ds the you want keep in your empty ds. If None all coordinates are kept from original ds. The default is None.

Returns:

empty_ds – dataset with only coordinate information

Return type:

xr.Dataset

nlmod.util.get_exe_path(exe_name='mf6', bindir=None, download_if_not_found=True, version_tag=None, repo='executables')[source]

Get the full path of the executable.

Searching for the executables is done in the following order: 0. If exe_name is a full path, return the full path of the executable. 1. The directory specified with bindir. Raises error if exe_name is provided

and not found.

  1. The directory used by nlmod installed in this environment.

  2. If the executables were downloaded with flopy/nlmod from an other env, most recent installation location of MODFLOW is found in flopy metadata

Else: 4. Download the executables using version_tag and repo.

The returned directory is checked to contain exe_name if it is provided.

Parameters:
  • exe_name (str, optional) – The name of the executable, by default “mf6”.

  • bindir (Path, optional) – The directory where the executables are stored, by default None

  • download_if_not_found (bool, optional) – Download the executables if they are not found, by default True.

  • repo (str, default "executables") – Name of GitHub repository. Choose one of “executables” (default), “modflow6”, or “modflow6-nightly-build”. If repo and version_tag are provided the most recent installation location of MODFLOW is found in flopy metadata that respects version_tag and repo. If not found, the executables are downloaded using repo and version_tag.

  • version_tag (str, default None) – GitHub release ID: for example “18.0” or “latest”. If repo and version_tag are provided the most recent installation location of MODFLOW is found in flopy metadata that respects version_tag and repo. If not found, the executables are downloaded using repo and version_tag.

Returns:

exe_full_path – full path of the executable.

Return type:

str

nlmod.util.get_flopy_bin_directories(version_tag=None, repo='executables')[source]

Get the directories where the executables are stored.

Obtain the bin directory installed with flopy. If enable_version_check is True, all installation location of MODFLOW are found in flopy metadata that respects version_tag and repo.

Parameters:
  • repo (str, default "executables") – Name of GitHub repository. Choose one of “executables” (default), “modflow6”, or “modflow6-nightly-build”. If repo and version_tag are provided the most recent installation location of MODFLOW is found in flopy metadata that respects version_tag and repo. If not found, the executables are downloaded using repo and version_tag.

  • version_tag (str, default None) – GitHub release ID: for example “18.0” or “latest”. If repo and version_tag are provided the most recent installation location of MODFLOW is found in flopy metadata that respects version_tag and repo. If not found, the executables are downloaded using repo and version_tag.

Returns:

list of directories where the executables are stored.

Return type:

list

nlmod.util.get_google_drive_filename(fid, timeout=120)[source]

Get the filename of a google drive file.

Parameters:

fid (str) – google drive id name of a file.

Returns:

file_name – filename.

Return type:

str

nlmod.util.get_model_dirs(model_ws)[source]

Creates a new model workspace directory, if it does not exists yet. Within the model workspace directory a few subdirectories are created (if they don’t exist yet):

  • figure

  • cache

Parameters:

model_ws (str) – model workspace.

Returns:

  • figdir (str) – figure directory inside model workspace.

  • cachedir (str) – cache directory inside model workspace.

nlmod.util.import_hydropandas(method='nlmod.read.knmi.download_knmi()')[source]

Import hydropandas, raise an error if hydopandas is not installed.

Parameters:

method (str, optional) – The name of the method that is displayed in the error message. The default is “nlmod.read.knmi.download_knmi()”.

Raises:

ImportError – When hydropandas is not installed, an ImportError is raised.

Returns:

hpd – The hydropandas module.

Return type:

module

nlmod.util.polygon_from_extent(extent)[source]

Create a shapely polygon from a given extent.

Parameters:

extent (tuple, list or array) – extent (xmin, xmax, ymin, ymax).

Returns:

polygon_ext – polygon of the extent.

Return type:

shapely.geometry.polygon.Polygon

nlmod.util.zonal_statistics(gdf, da, columns=None, buffer=0.0, engine='geocube', all_touched=True, statistics='mean', add_to_gdf=True, progressbar=False)[source]

Calculate raster statistics in the features of a GeoDataFrame.

Parameters:
  • gdf (gpd.GeoDataFrame) – A GeoDataFrame with features for which to calculate the statistics.

  • da (xr.DataArray or str) – A DataArray of the raster. When engine==’rasterstats’, da can also be the filename of the raster.

  • columns (str or list of str, optional) – The columns in gdf to add the statistics to. When columns is None, use the name of the statistics. When columns is a string, and more than one statistic is calculated, the value of columns is used as the prefix before each statistic, to form the column-names. The default is None.

  • buffer (float, optional) – The buffer, in m, that is added to each of the features of gdf, before calculating the statistics. The default is 0.0.

  • engine (str, optional) – The engine to use for the calculation of the statistics. The possible values are ‘geocube’ and ‘rasterstats’. The two engines should approximately result in the same values. If features overlap (which will happen more frequently when buffer > 0), the result will differ. If engine==’geocube’, each raster-cell is only designated to one feature. If engine=’rasterstats’ a raster-cell can be designated to multiple features, if features overlap. The default is “geocube”.

  • all_touched (bool, optional) – If True, include every raster cell touched by a geometry. Otherwise only include those having a center point within the polygon. The default is True.

  • statistics (str or list of str, optional) – The name or a list of names of the statistics to be calculated. The default is “mean”.

  • add_to_gdf (bool, optional) – Add the result to the orignal GeoDataFrame if True. Otherwise return a GeoDataFrame with only the statistics. The default is True.

  • progressbar (bool, optional) – show progressbar when using rasterstats. The default is False.

Returns:

A GeoDataFrame containing the the statistics in some of its columns.

Return type:

gpd.GeoDataFrame

cache

class nlmod.cache.NumpyEncoder(*, skipkeys=False, ensure_ascii=True, check_circular=True, allow_nan=True, sort_keys=False, indent=None, separators=None, default=None)[source]

Special json encoder for numpy types.

default(o)[source]

Implement this method in a subclass such that it returns a serializable object for o, or calls the base implementation (to raise a TypeError).

For example, to support arbitrary iterators, you could implement default like this:

def default(self, o):
    try:
        iterable = iter(o)
    except TypeError:
        pass
    else:
        return list(iterable)
    # Let the base class default method raise the TypeError
    return super().default(o)
nlmod.cache._check_for_data_array(ds)[source]

Check if the saved NetCDF-file represents a DataArray or a Dataset, and return this data-variable.

The file contains a DataArray when a variable called “__xarray_dataarray_variable__” is present in the Dataset. If so, return a DataArray, otherwise return the Dataset.

By saving the DataArray, the coordinate “spatial_ref” was saved as a separate variable. Therefore, add this variable as a coordinate to the DataArray again.

Parameters:

ds (xr.Dataset) – Dataset with dimensions and coordinates.

Returns:

ds – A Dataset or DataArray containing the cached data.

Return type:

xr.Dataset or xr.DataArray

nlmod.cache._explicit_dataset_coordinate_comparison(ds_in, ds_cache)[source]

Perform explicit dataset coordinate comparison.

Uses xarray.testing.assert_identical().

Parameters:
  • ds_in (xr.Dataset) – Input dataset.

  • ds_cache (xr.Dataset) – Cached dataset.

Returns:

True if coordinates are identical, else False.

Return type:

bool

Raises:

AssertionError – If the coordinates are not equal.

nlmod.cache._get_modification_time(func)[source]

Return the modification time of the module where func is defined.

Parameters:

func (function) – function.

Returns:

modification time of module.

Return type:

float

nlmod.cache._same_function_arguments(func_args_dic, func_args_dic_cache)[source]

Checks if two dictionaries with function arguments are identical.

The following items are checked:
  1. if they have the same keys

  2. if the items have the same type

  3. if the items have the same values (only implemented for the types: int, float, bool, str, bytes, list, tuple, dict, np.ndarray, xr.DataArray, flopy.mf6.ModflowGwf).

Parameters:
  • func_args_dic (dictionary) – dictionary with all the args and kwargs of a function call.

  • func_args_dic_cache (dictionary) – dictionary with all the args and kwargs of a previous function call of which the results are cached.

Returns:

if True the dictionaries are identical which means that the cached data was created using the same function arguments as the requested data.

Return type:

bool

Notes

Keys that end with ‘_hash’ are assumed to be hashes and not function arguments. They are checked equally.

nlmod.cache._update_docstring_and_signature(func)[source]

Add function arguments ‘cachedir’ and ‘cachename’ to the docstring and signature of a function.

The function arguments are added before the “Returns” header in the docstring. If the function has no Returns header in the docstring, the function arguments are not added to the docstring.

Parameters:

func (function) – function that is decorated.

Return type:

None

nlmod.cache.cache_netcdf(coords_2d=False, coords_3d=False, coords_time=False, attrs_ds=False, datavars=None, coords=None, attrs=None, nc_hash=True)[source]

Decorator to read/write the result of a function from/to a file to speed up function calls with the same arguments. Should only be applied to functions that:

  • return an Xarray Dataset

  • have no more than one xarray dataset as function argument

  • have functions arguments of types that can be checked using the

_is_valid_cache functions

1. The directory and filename of the cache should be defined by the person calling a function with this decorator. If not defined no cache is created nor used. 2. Create a new cached file if it is impossible to check if the function arguments used to create the cached file are the same as the current function arguments. This can happen if one of the function arguments has a type that cannot be checked using the _is_valid_cache function. 3. Function arguments are pickled together with the cache to check later if the cache is valid. 4. If one of the function arguments is an xarray Dataset it is not pickled. Therefore we cannot check if this function argument is identical for the cached data and the new function call. We do check if the xarray Dataset coördinates correspond to the coördinates of the cached netcdf file. 5. This function uses functools.wraps and some home made magic in _update_docstring_and_signature to add arguments of the decorator to the decorated function. This assumes that the decorated function has a docstring with a “Returns” heading. If this is not the case an error is raised when trying to decorate the function.

If all kwargs are left to their defaults, the function caches the full dataset.

Parameters:
  • ds (xr.Dataset) – Dataset with dimensions and coordinates.

  • coords_2d (bool, optional) – Shorthand for adding 2D coordinates. The default is False.

  • coords_3d (bool, optional) – Shorthand for adding 3D coordinates. The default is False.

  • coords_time (bool, optional) – Shorthand for adding time coordinates. The default is False.

  • attrs_ds (bool, optional) – Shorthand for adding model dataset attributes. The default is False.

  • datavars (list, optional) – List of data variables to check for. The default is an empty list.

  • coords (list, optional) – List of coordinates to check for. The default is an empty list.

  • attrs (list, optional) – List of attributes to check for. The default is an empty list.

  • nc_hash (bool, optional) – check if the pickled function arguments belong to the cached netcdf file. Default is True.

nlmod.cache.cache_pickle(func)[source]

Decorator to read/write the result of a function from/to a file to speed up function calls with the same arguments. Should only be applied to functions that:

  • return a picklable object

  • have functions arguments of types that can be checked using the

_is_valid_cache functions

1. The directory and filename of the cache should be defined by the person calling a function with this decorator. If not defined no cache is created nor used. 2. Create a new cached file if it is impossible to check if the function arguments used to create the cached file are the same as the current function arguments. This can happen if one of the function arguments has a type that cannot be checked using the _is_valid_cache function. 3. Function arguments are pickled together with the cache to check later if the cache is valid. 4. This function uses functools.wraps and some home made magic in _update_docstring_and_signature to add arguments of the decorator to the decorated function. This assumes that the decorated function has a docstring with a “Returns” heading. If this is not the case an error is raised when trying to decorate the function.

nlmod.cache.clear_cache(cachedir, prompt=True)[source]

Clears the cache in a given cache directory by removing all .pklz and corresponding .nc files.

Parameters:
  • cachedir (str) – path to cache directory.

  • prompt (bool, optional) – Ask for confirmation before removing the cache. The default is True.

Return type:

None.

nlmod.cache.ds_contains(ds, coords_2d=False, coords_3d=False, coords_time=False, attrs_ds=False, datavars=None, coords=None, attrs=None)[source]

Returns a Dataset containing only the required data.

If all kwargs are left to their defaults, the function returns the full dataset.

Parameters:
  • ds (xr.Dataset) – Dataset with dimensions and coordinates.

  • coords_2d (bool, optional) – Shorthand for adding 2D coordinates. The default is False.

  • coords_3d (bool, optional) – Shorthand for adding 3D coordinates. The default is False.

  • coords_time (bool, optional) – Shorthand for adding time coordinates. The default is False.

  • attrs_ds (bool, optional) – Shorthand for adding model dataset attributes. The default is False.

  • datavars (list, optional) – List of data variables to check for. The default is an empty list.

  • coords (list, optional) – List of coordinates to check for. The default is an empty list.

  • attrs (list, optional) – List of attributes to check for. The default is an empty list.

Returns:

ds – A Dataset containing only the required data.

Return type:

xr.Dataset

nlmod.cache.hash_xarray_coords(ds, include_metadata: bool = False)[source]

Create a hash of xarray coordinate(s) using array bytes and optionally metadata.

Parameters:
  • coord (xarray.core.coordinates.Coordinate) – The xarray coordinate object to hash

  • include_metadata (bool, optional) – Whether to include metadata in the hash. Default is False.

Returns:

The hexadecimal hash string

Return type:

str

nlmod.cache.hash_xarray_data_vars(ds, include_metadata: bool = False)[source]

Create a hash of xarray data variables using array bytes and optionally metadata.

Parameters:
  • ds (xarray.Dataset or xarray.DataArray) – The xarray coordinate object to hash

  • include_metadata (bool, optional) – Whether to include metadata in the hash. Default is False.

Returns:

The hexadecimal hash string

Return type:

str