imod.select - Get points and cross sections

Get points, cross sections and layers.

imod.select.cross_section_line(data, start, end)

Obtain an interpolated cross-sectional slice through gridded data. Utilizing the interpolation functionality in xarray, this function takes a vertical cross-sectional slice along a line through the given data on a regular (possibly non-equidistant) grid, which is given as an xarray.DataArray so that we can utilize its coordinate data.

Adapted from Metpy: https://github.com/Unidata/MetPy/blob/master/metpy/interpolate/slices.py

Parameters

• data (xarray.DataArray or xarray.Dataset) –

  Three- (or higher) dimensional field(s) to interpolate. The DataArray (or each DataArray in the Dataset) must have been parsed by MetPy and include both an x and y coordinate dimension and the added crs

  coordinate.

• start ((2, ) array_like) – A latitude-longitude pair designating the start point of the cross section.

• end ((2, ) array_like) – A latitude-longitude pair designating the end point of the cross section.

Returns

The interpolated cross section, with new dimension “s” along the cross-section. The cellsizes along “s” are given in the “ds” coordinate.

Return type

xarray.DataArray or xarray.Dataset

imod.select.cross_section_linestring(data, linestring)

Obtain an interpolated cross-sectional slice through gridded data. Utilizing the interpolation functionality in xarray, this function takes a vertical cross-sectional slice along a linestring through the given data on a regular grid, which is given as an xarray.DataArray so that we can utilize its coordinate data.

Adapted from Metpy: https://github.com/Unidata/MetPy/blob/master/metpy/interpolate/slices.py

Parameters

• data (xarray.DataArray or xarray.Dataset) –

  Three- (or higher) dimensional field(s) to interpolate. The DataArray (or each DataArray in the Dataset) must have been parsed by MetPy and include both an x and y coordinate dimension and the added crs

  coordinate.

• linestring (shapely.geometry.LineString) –

  Shapely geometry designating the linestring along which to sample the cross section.

  Note that a LineString can easily be taken from a geopandas.GeoDataFrame using the .geometry attribute. Please refer to the examples.

Returns

The interpolated cross section, with new index dimension along the cross-section.

Return type

xarray.DataArray or xarray.Dataset

Examples

Load a shapefile (that you might have drawn before using a GIS program), take a linestring from it, and use it to extract the data for a cross section.

>>> geodataframe = gpd.read_file("cross_section.shp")
>>> linestring = geodataframe.geometry[0]
>>> section = cross_section_linestring(data, linestring)

Or, construct the linestring directly in Python:

>>> import shapely.geometry as sg
>>> linestring = sg.LineString([(0.0, 1.0), (5.0, 5.0), (7.5, 5.0)])
>>> section = cross_section_linestring(data, linestring)

If you have drawn multiple cross sections within a shapefile, simply loop over the linestrings:

>>> sections = [cross_section_linestring(data, ls) for ls in geodataframe.geometry]

imod.select.points_in_bounds(da, **points)

Returns whether points specified by keyword arguments fall within the bounds of da.

Parameters

• da (xr.DataArray) –

• points (keyword arguments of coordinate=values) – keyword arguments specifying coordinate and values. Please refer to the examples.

Returns

in_bounds

Return type

np.array of bools

Examples

Create the DataArray, then use the keyword arguments to define along which coordinate to check whether the points are within bounds.

>>> nrow, ncol = 3, 4
>>> data = np.arange(12.0).reshape(nrow, ncol)
>>> coords = {"x": [0.5, 1.5, 2.5, 3.5], "y": [2.5, 1.5, 0.5]}
>>> dims = ("y", "x")
>>> da = xr.DataArray(data, coords, dims)
>>> x = [0.4, 2.6]
>>> points_in_bounds(da, x=x)

This works for an arbitrary number of coordinates:

>>> y = [1.3, 2.7]
>>> points_in_bounds(da, x=x, y=y)

imod.select.points_indices(da, **points)

Get the indices for points as defined by the arrays x and y.

This function will raise a ValueError if the points fall outside of the bounds of the DataArray to avoid undefined behavior. Use the imod.select.points_in_bounds function to detect these points.

Parameters

• da (xr.DataArray) –

• points (keyword arguments of coordinates and values) –

Returns

indices

Return type

dict of {coordinate: xr.DataArray with indices}

Examples

To extract values:

>>> x = [1.0, 2.2, 3.0]
>>> y = [4.0, 5.6, 7.0]
>>> indices = imod.select.points_indices(da, x=x, y=y)
>>> ind_y = indices["y"]
>>> ind_x = indices["x"]
>>> selection = da.isel(x=ind_x, y=ind_y)

Or shorter, using dictionary unpacking:

>>> indices = imod.select.points_indices(da, x=x, y=y)
>>> selection = da.isel(**indices)

To set values (in a new array), the following will do the trick:

>>> empty = xr.full_like(da, np.nan)
>>> empty.data[indices["y"].values, indices["x"].values] = values_to_set

Unfortunately, at the time of writing, xarray’s .sel method does not support setting values yet. The method here works for both numpy and dask arrays, but you’ll have to manage dimensions yourself!

The imod.select.points_set_values() function will take care of the dimensions.

imod.select.points_set_values(da, values, **points)

Set values at specified points.

This function will raise a ValueError if the points fall outside of the bounds of the DataArray to avoid undefined behavior. Use the imod.select.points_in_bounds function to detect these points.

Parameters

• da (xr.DataArray) –

• values ((int, float) or array of (int, float)) –

• points (keyword arguments of coordinate=values) – keyword arguments specifying coordinate and values.

Returns

None

Return type

mutates da.

Examples

>>> x = [1.0, 2.2, 3.0]
>>> y = [4.0, 5.6, 7.0]
>>> values = [10.0, 11.0, 12.0]
>>> selection = imod.select.points_values(da, values, x=x, y=y)

imod.select.points_values(da, **points)

Get values from specified points.

This function will raise a ValueError if the points fall outside of the bounds of the DataArray to avoid undefined behavior. Use the imod.select.points_in_bounds function to detect these points.

Parameters

• da (xr.DataArray) –

• points (keyword arguments of coordinate=values) – keyword arguments specifying coordinate and values.

Returns

selection

Return type

xr.DataArray

Examples

>>> x = [1.0, 2.2, 3.0]
>>> y = [4.0, 5.6, 7.0]
>>> selection = imod.select.points_values(da, x=x, y=y)

imod.select.upper_active_layer(da, is_ibound=True, include_constant_head=False)

Function to get the upper active layer from ibound xr.DataArray

Parameters: da : 3d xr.DataArray is_ibound: da is interpreted as ibound, with values 0 inactive, 1, active, -1 chd

if False: upper_active_layer is interpreted as upper layer that has data

include_constant_head : also incluse chd cells? bool, default False

Returns: 2d xr.DataArray of layernumber of upper active model layer