gridmarthe.grid.processing.geometry

Module to manage geometry attributes of Marthe grids/domain

gridmarthe.grid.processing.geometry.get_active_mask(ds, varname='permeab', nanval=[-9999.0, 0.0], as_array=False, only_mask=False, shp_file=None, epsg=27572)[source]

Get the mask of active domain from hydraulic conductivity variable

This function (i) get the active mask as a 0/1 array and optionally (ii) filter the dataset on valid values and dissolve results to get a mask shape

Input ds should be the permh dataset (read from permh file, ie Horizontal hydraulic conductivity, without the dropnan option).

Parameters:

dsxarray.Dataset
varnamestr, optional: default is ‘permeab’
nanvalfloat or list, optional: default are ‘permeab’ nan values : 0, -9999.
as_arraybool, optional.: Option to get result as a xr.Dataset and not geodataframe. Default is False.
only_maskbool, optional: filter ds on active mask. Default is False, returns a dataset with ‘ibound’ variable set to 1 (active domain) or 0.
shp_filestr, optional.: if set (and not as_array), used to stored result in a file.
epsgint, optional: if shp_file, use epsg to set projection.

Returns:

xr.Dataset: Dataset with ibound field, or gpd.GeoDataFrame of active domain if as_array is set to False.

gridmarthe.grid.processing.geometry.compute_geometry(topo, hsubs, mask=None, topo_varname='h_topogr', subs_varname='h_substrat')[source]

Compute geometry attributes of Marthe domain

Parameters:

topoxarray.Dataset: Topgraphy of the domain (stored in the first layer, in Marthe Conventions).
hsubsxarray.Dataset: altitude of all the lower boundary in the domain
masknumpy.array, optional: list of indices (zone) to keep, if None (default) not used. It is recommended to use this mask to avoid computing on invalid cells. For example, values may be defined in masked cells of the model domain, which will lead to incorrect results. Using the active domain as mask is a good practice (See example).
topo_varnamestr, optional: name of the variable containing the topography in the corresponding dataset, allow custom name for marthe backward compatibility
subs_varnamestr, optional: name of the variable containing the substratum in the corresponding dataset, allow custom name for marthe backward compatibility

Returns:

xarray.Dataset: A new dataset with layer, depth, thickness, upper/lower altitude.

Notes

If mask is not provided, the input datasets should be sliced on valid cells before calling this function for accurate results;
The time dimension is dropped during process

Examples

>>> import gridmarthe as gm
>>> permh = gm.load_marthe_grid('data/craie_npc.permh', drop_nan=True)
>>> topo  = gm.load_marthe_grid('data/craie_npc.topog', varname='H_TOPOGR')
>>> hsub  = gm.load_marthe_grid('data/craie_npc.hsubs', varname='H_SUBSTRAT')
>>> geom = gm.compute_geometry(topo, hsub, mask=permh.zone.data)

gridmarthe.grid.processing.geometry.get_surface_layer(ds, aquif_layers=None)[source]

Compute surface mask of marthe domain

This function return min layer for every zone of a grimarthe dataset with z coords A subset on specific (aquifers) layers can be performed with aquif_layers. if set, aquif_layers must be a sequence (list, tuple, array) of layer (list of int).

This should be used to get a surface mask, ie get zone to filter a dataset.

Parameters:

dsxarray.Dataset
aquif_layerssequence (list, tuple, array) of int: representing layers to subset ds. Only active domain must be passed to function (ie drop nan first)

Returns:

surface_mask: xarray.Dataset

Examples

>>>    mask = get_surface_layer(ds, [6,8,9])
>>>    ds_surf = ds.sel(zone=mask.zone.data)

gridmarthe.grid.processing.geometry.search_zone(ds, i=None, j=None, x=None, y=None, z=None)[source]

Search zone number in marthe grid, based on xy or ij (col, lig)

This function can be used to search zone number from coordinates or indices. You must provide either (i,j) or (x,y).

Parameters:

dsxarray.Dataset: dataset with zone, x, y, dx, dy variables.
iint, optional: column index to search zone.
jint, optional: row index to search zone.
xfloat, optional: x coordinate to search zone.
yfloat, optional: y coordinate to search zone.
zint, optional: layer index to search zone. If not provided, all layers are considered.

Returns:

zonexarray.Dataset: dataset with zone variable, containing the zone number(s) corresponding to the provided coordinates. If no zone is found, an empty dataset is returned. If multiple zones are found, all of them are returned.

Notes

if ds is multilayered, you need to provide the layer you want (z arg., int type)
ds should contains dx and dy
ds should not have assigned coords (x and y are variables, zone is the

dimension coordinates (with time))