Exemple d’utilisation de gridmarthe [FR]

Le package gridmarthe est conçu pour faciliter la lecture/l’écriture de grille au format Marthe. Ce notebook permet d’explorer quelques fonctionnalités de base du package et de montrer un exemple de traitement des grilles.

# import du package
import gridmarthe as gm

Fonctionnalités de base

Les grilles contiennent certaines métadonnées, dont l’index du pas de temps mais pas les dates en elles-même. Aussi, il est utile de fournir à la fonction de lecture soit une série de date (à construire manuellement) soit un fichier pastp pour que les dates soient lues automatiquement.

En complément, plusieurs arguments peuvent être utiles à la lecture (retrait des valeurs nulles, transformation des xy, ajout de l’indicateur de maillage ou des numéros de colonnes/lignes, etc.).

Pour évaluer les options possibles, on peut afficher l’aide de la fonction.

help(gm.load_marthe_grid)

Help on function load_marthe_grid in module gridmarthe.gridmarthe:

load_marthe_grid(filename: str, varname: Optional[str] = None, dates=None, fpastp: Optional[str] = None, nanval: Union[int, float, NoneType] = None, drop_nan: bool = False, xyfactor: Union[int, float] = 1.0, keepligcol: bool = False, add_id_grid: bool = False, title: Optional[str] = None, var_attrs: dict = {}, model_attrs: dict = {'resolution_units': 'm', 'projection': 'epsg:27572', 'domain': 'FR-France'}, verbose: bool = False)
    Read Marthe Grid File as xarray.Dataset 
    
    The gridfile is read as a sequence: the variable for all layer
    for main grid, then all layer for nested grids, is stored in
    a 1D vector for every timestep. A single spatial identifier
    `zone` is used to map spatial coordinates.
    
    Before plot operations, user can assign coordinates (set x,y
    as dimension coordinates and drop zone) to get 2-D arrays (or
    3D arrays if multilayer) for every timesteps.
    
    
    Parameters
    ----------
    
        filename: str
            A path to marthegrid file (*.permh, *.out, etc.)
        
        varname : str, Optionnal
            variable to access in martgrid file, e.g `CHARGE` for groundwater head. See marthegrid file content.
            if None  is passed (default), function will scan all varnames in filename and keep first only
            if 'all' is passed,  function will scan all varnames in filename and keep all. 
            All datavars are added to dataset, using recursive call to func
            if wrong variable name is passed, empty data will be returned.
        
        dates: sequence, Optionnal
            Can be a pd.date_range, pd.Series, pd.DatetimeIndex, np.array or list of datetime/np.datetime objects.
            If no dates (or no fpastp) is provided, a fake sequence of dates from 1850 to 1900 will
            be used for xarray object
        
        fpastp: str, Optionnal
            A pastp file to read for dates
        
        nanval: float, Optionnal
            A code value for nan values. Default is 9999.
        
        drop_nan: bool, Optionnal
            Drop nan values (corresponding to nanval) in xarray object to return.
            Default is False (keep nan values).
        
        xyfactor: int or float, Optionnal
            factor to transform X and Y values. e.g.: 1000 to convert km XY to meters.
            Default is 1.
        
        keepligcol: bool, Optionnal
            Add columns (col) and rows (lig) index (from 1 to n), Default is False.
        
        add_id_grid: bool, Optionnal
            Add grid id (from 0 to n), useful for nested grids.
            0 is main grid, Default is False
        
        title: str , Optionnal
            Title for grid attributes. Default is None (not used)
        
        var_attrs: dict, Optionnal
            Dictionnary of attributes to add to variable DataArray.
        
        model_attrs: dict, Optionnal
            Dictionnary of attributes to add to Dataset.
            by default, gis attrs are added and can be modified
            >>>    'resolution_units': 'm',
            >>>    'projection'      : 'epsg:27572',
            >>>    'domain'          : 'FR-France',
        
        verbose: bool, Optionnal
            Print some information about execution in stdout.
            Default is False.
    
    Returns
    -------
        ds: xr.Dataset
            A xarray.Dataset object containing values and attributes read from Marthe grid file.

# lecture d'une grille marthe et stockage dans un Dataset (librairie xarray)
ds = gm.load_marthe_grid('./chasim_hallue.out', 'CHARGE', fpastp='./hallue.pastp', drop_nan=True) # si la valeur des nans n'est pas 9999 (ex. permh, NaNs = 0), on peut spécifier la valeur

# affichage des valeurs et attributs
ds

<xarray.Dataset> Size: 784kB
Dimensions:  (time: 205, zone: 927)
Coordinates:
  * time     (time) datetime64[ns] 2kB 1995-07-31 1995-08-01 ... 2012-07-01
  * zone     (zone) int64 7kB 255 256 257 258 259 ... 2722 2723 2724 2725 2726
Data variables:
    charge   (time, zone) float32 760kB 100.0 100.6 101.1 ... 27.0 26.0 26.35
    x        (zone) float32 4kB 617.8 618.2 618.8 619.2 ... 606.8 607.2 607.8
    y        (zone) float32 4kB 2.567e+03 2.567e+03 ... 2.543e+03 2.543e+03
    dx       (zone) float32 4kB 0.5 0.5 0.5 0.5 0.5 0.5 ... 0.5 0.5 0.5 0.5 0.5
    dy       (zone) float32 4kB 0.5 0.5 0.5 0.5 0.5 0.5 ... 0.5 0.5 0.5 0.5 0.5
Attributes: (12/17)
    conventions:          CF-1.10
    title:                Modélisation du bassin de la SOMME Nappe_Libre
    marthe_grid_version:  9.0
    original_dimensions:  x,y,z [grids]: 53 54 1
    lon_resolution:       0.5
    lat_resolution:       0.5
    ...                   ...
    creation_date:        Created on 2025-03-26T12:54:45Z UTC
    institution:          BRGM, French Geological Survey, Orléans, France
    comment:              Hydrogeological model created with MARTHE code (Thi...
    resolution_units:     m
    projection:           epsg:27572
    domain:               FR-France

La grille est chargée dans un objet xarray 2D, de dimension (time, zone). La dimension spatiale est donc 1D, contenue dans zone. Chaque zone possède un couple de coordonnées xy, stocké en variable (et non en dimension). Ce format permet d’effecture de nombreuses opérations de manière plus efficaces qu’en 2 ou 3D + time. Il suffit d’assigner les coordonnées (transformation en 2 ou 3D + temps) avant les fonctionnalités graphiques.

# on ajoute les coordonées
ds_3d = gm.assign_coords(ds)
# a noter que les fonctions sont aussi accessibles par méthodes/attributs:
# ds_3d = ds.mart.assign_coords() # equivalent
ds_3d

<xarray.Dataset> Size: 2MB
Dimensions:  (y: 48, x: 44, time: 205)
Coordinates:
  * y        (y) float32 192B 2.543e+03 2.544e+03 ... 2.566e+03 2.567e+03
  * x        (x) float32 176B 599.2 599.8 600.2 600.8 ... 619.8 620.2 620.8
  * time     (time) datetime64[ns] 2kB 1995-07-31 1995-08-01 ... 2012-07-01
Data variables:
    charge   (time, y, x) float32 2MB nan nan nan nan ... 102.7 103.4 103.5
    dx       (y, x) float32 8kB nan nan nan nan nan nan ... 0.5 0.5 0.5 0.5 0.5
    dy       (y, x) float32 8kB nan nan nan nan nan nan ... 0.5 0.5 0.5 0.5 0.5
Attributes: (12/17)
    conventions:          CF-1.10
    title:                Modélisation du bassin de la SOMME Nappe_Libre
    marthe_grid_version:  9.0
    original_dimensions:  x,y,z [grids]: 53 54 1
    lon_resolution:       0.5
    lat_resolution:       0.5
    ...                   ...
    creation_date:        Created on 2025-03-26T12:54:45Z UTC
    institution:          BRGM, French Geological Survey, Orléans, France
    comment:              Hydrogeological model created with MARTHE code (Thi...
    resolution_units:     m
    projection:           epsg:27572
    domain:               FR-France
# on plot à un instant t
ds_3d.isel(time=-1)['charge'].plot.pcolormesh(x='x',y='y')

<matplotlib.collections.QuadMesh at 0x7d9c66137cd0>
../../_images/66ee1cad3b2912800353be385a23a0ed8a805d8c3e11fd6e9259bc695a3a4979.png 
# ou la chronique dans une maille
ds.isel(zone=255)['charge'].plot()

[<matplotlib.lines.Line2D at 0x7d9c66049f00>]
../../_images/dfd50bba10e1cc8a1dd57648126dbdf96743962a23877c15798a32250b4b0b06.png

Pour les gigognes, une fonction de plot spécifique est proposée (gm.plot_nested_grids()).

Pour les modèles avec plusieurs couches, on peut récupérer les couches de surfaces~: gm.get_min_layer(ds, [subset_layers=[3,5]])

ou même afficher ce résultat~: gm.plot_outcrop()

Les fonctions natives de xarray, numpy permettent d’ores-et-déjà d’accéder à certaines fonctionnalités de Operasem (substitution de valeurs, transformation, opérations entre deux grilles, etc.)

Fonctionnalités avancées

Exemple de calcul du SPLI (package pyspli, https://gitlab.brgm.fr/brgm/si-eau/traitements-et-valorisations/py_spli)

import numpy as np
import spli
from multiprocessing import cpu_count

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[8], line 2
      1 import numpy as np
----> 2 import spli
      3 from multiprocessing import cpu_count

ModuleNotFoundError: No module named 'spli'

ds_spli = spli.compute_spli(ds, n_jobs=cpu_count()) # monthly mean included, dummy ref_period: only for example
#ds_spli = spli.compute_spli(ds, ref_period=[1995, 2010], n_jobs=cpu_count()) # monthly mean included, dummy ref_period: only for example

ds_spli

<xarray.Dataset> Size: 2MB
Dimensions:  (zone: 927, time: 205)
Coordinates:
  * zone     (zone) int64 7kB 255 256 257 258 259 ... 2722 2723 2724 2725 2726
  * time     (time) datetime64[ns] 2kB 1995-07-01 1995-08-01 ... 2012-07-01
Data variables:
    spli     (time, zone) float64 2MB 0.06557 0.06549 0.06554 ... 0.01162 0.6706
Attributes:
    varname:        CHARGE
    units:          m
    missing_value:  9999.0
    standard_name:  water_table_level
    long_name:      groundwater head
# get coords back
ds_spli['x'] = ("zone", ds['x'].data)
ds_spli['y'] = ("zone", ds['y'].data)
if "z" in ds_spli.data_vars.keys():
    ds_spli['z'] = ("zone", ds['z'].data)

# assign coords for plot
ds_xy_spli = ds_spli.assign_coords(
    x=('zone', np.around(ds_spli['x'].data, 1) ),
    y=('zone', np.around(ds_spli['y'].data, 1) ),
)
ds_xy_spli = ds_xy_spli.set_index(zone=['y', 'x']).unstack('zone').sortby('y', 'x')
ds_xy_spli

<xarray.Dataset> Size: 3MB
Dimensions:  (y: 48, x: 44, time: 205)
Coordinates:
  * y        (y) float32 192B 2.543e+03 2.544e+03 ... 2.566e+03 2.567e+03
  * x        (x) float32 176B 599.2 599.8 600.2 600.8 ... 619.8 620.2 620.8
  * time     (time) datetime64[ns] 2kB 1995-07-01 1995-08-01 ... 2012-07-01
Data variables:
    spli     (time, y, x) float64 3MB nan nan nan nan ... 0.1627 0.519 0.6357
Attributes:
    varname:        CHARGE
    units:          m
    missing_value:  9999.0
    standard_name:  water_table_level
    long_name:      groundwater head
# plot res
da=ds_xy_spli['spli'].sel(time=['1997-10-01', '2001-05-01', '2005-02-01','2012-01-01'])
da.plot.pcolormesh(
    x="x", y="y",
    extend='both',
    levels=[-3, -1.2815, -0.8416, -0.253, 0.253, 0.8416, 1.2815, 3],
    colors=['red', 'red', 'orange', 'yellow', 'lime', 'cyan', 'blue', 'navy', 'navy'], # double extrema for extended colors
    col='time',
    cbar_kwargs={
        'label': 'SPLI', 
        # 'ticks_label': ['Très bas', 'Bas', 'Modérément\nbas', 'Autour de la\nmoyenne', 'Modérément\nhaut', 'Haut', 'Très haut']
    }
)

<xarray.plot.facetgrid.FacetGrid at 0x7f7b284d8070>
../../_images/38d04bb6c4e149ce0b6eac19c46499236ecca88e514f8494335b1a7294948154.png