Example of treatments available

This notebook shows how to use some extra features for your gridmarthe treatments

# import modules
import numpy as np
import gridmarthe as gm
import matplotlib.pyplot as plt

grid = gm.load_marthe_grid('./data/craie_npc.permh', drop_nan=True)  # NaN needs to be drop or set (change 0 to nan for permh for example)
grid

<xarray.Dataset> Size: 4MB
Dimensions:   (time: 1, zone: 101209)
Coordinates:
  * time      (time) datetime64[ns] 8B 1850-01-01
  * zone      (zone) int32 405kB 1 2 3 4 5 ... 101206 101207 101208 101209
Data variables:
    permeab   (time, zone) float64 810kB 1e-05 1e-05 1e-05 ... 0.0001 0.0001
    x         (zone) float32 405kB 6.372e+05 6.378e+05 ... 6.382e+05 6.38e+05
    y         (zone) float32 405kB 2.617e+06 2.617e+06 ... 2.592e+06 2.592e+06
    dx        (zone) float32 405kB 500.0 500.0 500.0 500.0 ... 100.0 100.0 100.0
    dy        (zone) float32 405kB 500.0 500.0 500.0 500.0 ... 100.0 100.0 100.0
    z         (zone) int32 405kB 1 1 1 1 1 1 1 1 1 1 1 ... 9 9 9 9 9 9 9 9 9 9 9
    zone_all  (zone) int32 405kB 47402 47403 47763 ... 2381150 2381151 2381457
Attributes: (12/17)
    conventions:          CF-1.10
    title:                
    marthe_grid_version:  9.0
    original_dimensions:  x,y,z [grids]: 387 304 10; 307 442 10
    lon_resolution:       100.0, 500.0
    lat_resolution:       100.0, 500.0
    ...                   ...
    creation_date:        Created on 2025-09-02T09:31:47Z 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

xarray.Dataset

• Dimensions:
  • time: 1
  • zone: 101209
• Coordinates: (2)
  • time
    (time)
    datetime64[ns]
    1850-01-01

    array(['1850-01-01T00:00:00.000000000'], dtype='datetime64[ns]')

  • zone
    (zone)
    int32
    1 2 3 4 ... 101207 101208 101209

    array([     1,      2,      3, ..., 101207, 101208, 101209],
          shape=(101209,), dtype=int32)

• Data variables: (7)
  • permeab
    (time, zone)
    float64
    1e-05 1e-05 1e-05 ... 0.0001 0.0001

    varname :

    PERMEAB

    units :

    m/s

    missing_value :

    0.0

    standard_name :

    long_name :

    aquifer_hydraulic_conductivity

    array([[1.e-05, 1.e-05, 1.e-05, ..., 1.e-04, 1.e-04, 1.e-04]],
          shape=(1, 101209))

  • x
    (zone)
    float32
    6.372e+05 6.378e+05 ... 6.38e+05

    units :

    m

    axis :

    X

    coverage_content_type :

    coordinate

    array([637250., 637750., 624250., ..., 638050., 638150., 638050.],
          shape=(101209,), dtype=float32)

  • y
    (zone)
    float32
    2.617e+06 2.617e+06 ... 2.592e+06

    units :

    m

    axis :

    Y

    coverage_content_type :

    coordinate

    array([2616850., 2616850., 2616350., ..., 2592350., 2592350., 2592250.],
          shape=(101209,), dtype=float32)

  • dx
    (zone)
    float32
    500.0 500.0 500.0 ... 100.0 100.0

    array([500., 500., 500., ..., 100., 100., 100.],
          shape=(101209,), dtype=float32)

  • dy
    (zone)
    float32
    500.0 500.0 500.0 ... 100.0 100.0

    array([500., 500., 500., ..., 100., 100., 100.],
          shape=(101209,), dtype=float32)

  • z
    (zone)
    int32
    1 1 1 1 1 1 1 1 ... 9 9 9 9 9 9 9 9

    units :

    -

    axis :

    Z

    positive :

    down

    standard_name :

    depth

    long_name :

    aquifer_layer

    array([1, 1, 1, ..., 9, 9, 9], shape=(101209,), dtype=int32)

  • zone_all
    (zone)
    int32
    47402 47403 ... 2381151 2381457

    array([  47402,   47403,   47763, ..., 2381150, 2381151, 2381457],
          shape=(101209,), dtype=int32)

• Indexes: (2)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['1850-01-01'], dtype='datetime64[ns]', name='time', freq=None))

  • zone
    PandasIndex

    PandasIndex(Index([     1,      2,      3,      4,      5,      6,      7,      8,      9,
               10,
           ...
           101200, 101201, 101202, 101203, 101204, 101205, 101206, 101207, 101208,
           101209],
          dtype='int32', name='zone', length=101209))

• Attributes: (17)

  conventions :

  CF-1.10

  title :

  marthe_grid_version :

  9.0

  original_dimensions :

  x,y,z [grids]: 387 304 10; 307 442 10

  lon_resolution :

  100.0, 500.0

  lat_resolution :

  100.0, 500.0

  scale_factor :

  1.0

  nested_grid :

  True

  extend :

  xymin : 543750.0 2526350.0; xymax: 736750.0 2677850.0

  period :

  1850-1850

  frequency :

  nan day(s)

  creation_date :

  Created on 2025-09-02T09:31:47Z UTC

  institution :

  BRGM, French Geological Survey, Orléans, France

  comment :

  Hydrogeological model created with MARTHE code (Thiery, D. 2020. Guidelines for MARTHE v7.8 computer codefor hydro-systems modelling. report BRGM/RP-69660-FR).

  resolution_units :

  m

  projection :

  epsg:27572

  domain :

  FR-France

Get surface mask

surf = gm.get_surface_layer(grid)
surf

<xarray.Dataset> Size: 925kB
Dimensions:   (zone: 21025)
Coordinates:
  * zone      (zone) int32 84kB 345 331 2139 2111 ... 13132 13062 19985 19965
Data variables:
    time      (zone) datetime64[ns] 168kB 1850-01-01 1850-01-01 ... 1850-01-01
    permeab   (zone) float64 168kB 1e-06 1e-06 0.003068 ... 1e-06 2e-05 2e-05
    x         (zone) float32 84kB 6.188e+05 6.188e+05 ... 6.494e+05 6.494e+05
    y         (zone) float32 84kB 2.6e+06 2.6e+06 ... 2.606e+06 2.606e+06
    dx        (zone) float32 84kB 500.0 500.0 500.0 500.0 ... 100.0 100.0 100.0
    dy        (zone) float32 84kB 500.0 500.0 500.0 500.0 ... 100.0 100.0 100.0
    z         (zone) int32 84kB 2 2 6 6 8 8 8 8 5 2 6 ... 2 2 2 2 2 2 2 2 2 4 4
    zone_all  (zone) int32 84kB 178171 177784 647989 ... 1389346 1660427 1660120

xarray.Dataset

• Dimensions:
  • zone: 21025
• Coordinates: (1)
  • zone
    (zone)
    int32
    345 331 2139 ... 13062 19985 19965

    array([  345,   331,  2139, ..., 13062, 19985, 19965],
          shape=(21025,), dtype=int32)

• Data variables: (8)
  • time
    (zone)
    datetime64[ns]
    1850-01-01 ... 1850-01-01

    array(['1850-01-01T00:00:00.000000000', '1850-01-01T00:00:00.000000000',
           '1850-01-01T00:00:00.000000000', ...,
           '1850-01-01T00:00:00.000000000', '1850-01-01T00:00:00.000000000',
           '1850-01-01T00:00:00.000000000'],
          shape=(21025,), dtype='datetime64[ns]')

  • permeab
    (zone)
    float64
    1e-06 1e-06 ... 2e-05 2e-05

    array([1.00000e-06, 1.00000e-06, 3.06763e-03, ..., 1.00000e-06,
           2.00000e-05, 2.00000e-05], shape=(21025,))

  • x
    (zone)
    float32
    6.188e+05 6.188e+05 ... 6.494e+05

    array([618750., 618750., 618750., ..., 649350., 649350., 649350.],
          shape=(21025,), dtype=float32)

  • y
    (zone)
    float32
    2.6e+06 2.6e+06 ... 2.606e+06

    array([2599850., 2600350., 2600850., ..., 2606050., 2606150., 2606250.],
          shape=(21025,), dtype=float32)

  • dx
    (zone)
    float32
    500.0 500.0 500.0 ... 100.0 100.0

    array([500., 500., 500., ..., 100., 100., 100.],
          shape=(21025,), dtype=float32)

  • dy
    (zone)
    float32
    500.0 500.0 500.0 ... 100.0 100.0

    array([500., 500., 500., ..., 100., 100., 100.],
          shape=(21025,), dtype=float32)

  • z
    (zone)
    int32
    2 2 6 6 8 8 8 8 ... 2 2 2 2 2 2 4 4

    array([2, 2, 6, ..., 2, 4, 4], shape=(21025,), dtype=int32)

  • zone_all
    (zone)
    int32
    178171 177784 ... 1660427 1660120

    array([ 178171,  177784,  647989, ..., 1389346, 1660427, 1660120],
          shape=(21025,), dtype=int32)

• Indexes: (1)
  • zone
    PandasIndex

    PandasIndex(Index([  345,   331,  2139,  2111,  4730,  4693,  4655,  4616,  1008,   397,
           ...
           13569, 13502, 13430, 13356, 13282, 13207, 13132, 13062, 19985, 19965],
          dtype='int32', name='zone', length=21025))

• Attributes: (0)

# xarray version
toto = gm.assign_coords(surf, add_lay=False) # add lay false because it needs to be a variable to plot, not a dimension
gm.plot_outcrop(toto)
plt.show()

This function also allow user to get the value of a variable in the surface layer, or a subset of aquifer layers.

For example, to get the value of your variable (groundwater head for example) in the first layer encountered between 3 layers:

head_surf = gm.get_surface_layer(grid, aquif_layers=[6,8,9])  # get the values in the first layer between layers 6, 8 and 9

Get mask of active domain

mask = gm.get_active_mask(grid) # return a geopandas geodataframe

/home/docs/checkouts/readthedocs.org/user_builds/gridmarthe/envs/v0.1.3/lib/python3.10/site-packages/gridmarthe/operasem/gis.py:85: FutureWarning: The return type of `Dataset.dims` will be changed to return a set of dimension names in future, in order to be more consistent with `DataArray.dims`. To access a mapping from dimension names to lengths, please use `Dataset.sizes`.
  if 'time' in ds.dims.keys():

mask.boundary.plot(color='k')
plt.show()

Compute depths and thickness

Users can get geometry attributes (depth, thickness, upper/lower altitudes) in a dataset:

topo = gm.load_marthe_grid('./data/example.topog')   # be careful with nan here, topography can be set outside of active area
hsubs = gm.load_marthe_grid('./data/example.hsubs')  # same here, hsubs can be set outside of active area
geom = gm.compute_geometry(topo, hsubs)
# print(geom.where(~geom['depth'].isnull(), drop=True))
geom

<xarray.Dataset> Size: 75MB
Dimensions:    (time: 1, zone: 1176480)
Coordinates:
  * time       (time) datetime64[ns] 8B 1850-01-01
  * zone       (zone) int32 5MB 1 2 3 4 5 ... 1176477 1176478 1176479 1176480
Data variables:
    h_topogr   (time, zone) float64 9MB 9.999e+03 9.999e+03 ... 9.999e+03
    x          (zone) float32 5MB 5.438e+05 5.442e+05 ... 7.362e+05 7.368e+05
    y          (zone) float32 5MB 2.678e+06 2.678e+06 ... 2.526e+06 2.526e+06
    dx         (zone) float32 5MB 500.0 500.0 500.0 500.0 ... 500.0 500.0 500.0
    dy         (zone) float32 5MB 500.0 500.0 500.0 500.0 ... 500.0 500.0 500.0
    z          (zone) int32 5MB 1 1 1 1 1 1 1 1 1 ... 10 10 10 10 10 10 10 10 10
    z_lower    (time, zone) float64 9MB nan nan nan nan nan ... nan nan nan nan
    z_upper    (time, zone) float64 9MB nan nan nan nan nan ... nan nan nan nan
    thickness  (time, zone) float64 9MB nan nan nan nan nan ... nan nan nan nan
    depth      (time, zone) float64 9MB nan nan nan nan nan ... nan nan nan nan
Attributes: (12/17)
    conventions:          CF-1.10
    title:                
    marthe_grid_version:  9.0
    original_dimensions:  x,y,z [grids]: 387 304 10
    lon_resolution:       500.0
    lat_resolution:       500.0
    ...                   ...
    creation_date:        Created on 2025-09-02T09:31:50Z 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

xarray.Dataset

• Dimensions:
  • time: 1
  • zone: 1176480
• Coordinates: (2)
  • time
    (time)
    datetime64[ns]
    1850-01-01

    array(['1850-01-01T00:00:00.000000000'], dtype='datetime64[ns]')

  • zone
    (zone)
    int32
    1 2 3 4 ... 1176478 1176479 1176480

    array([      1,       2,       3, ..., 1176478, 1176479, 1176480],
          shape=(1176480,), dtype=int32)

• Data variables: (10)
  • h_topogr
    (time, zone)
    float64
    9.999e+03 9.999e+03 ... 9.999e+03

    array([[9999., 9999., 9999., ..., 9999., 9999., 9999.]],
          shape=(1, 1176480))

  • x
    (zone)
    float32
    5.438e+05 5.442e+05 ... 7.368e+05

    units :

    m

    axis :

    X

    coverage_content_type :

    coordinate

    array([543750., 544250., 544750., ..., 735750., 736250., 736750.],
          shape=(1176480,), dtype=float32)

  • y
    (zone)
    float32
    2.678e+06 2.678e+06 ... 2.526e+06

    units :

    m

    axis :

    Y

    coverage_content_type :

    coordinate

    array([2677850., 2677850., 2677850., ..., 2526350., 2526350., 2526350.],
          shape=(1176480,), dtype=float32)

  • dx
    (zone)
    float32
    500.0 500.0 500.0 ... 500.0 500.0

    array([500., 500., 500., ..., 500., 500., 500.],
          shape=(1176480,), dtype=float32)

  • dy
    (zone)
    float32
    500.0 500.0 500.0 ... 500.0 500.0

    array([500., 500., 500., ..., 500., 500., 500.],
          shape=(1176480,), dtype=float32)

  • z
    (zone)
    int32
    1 1 1 1 1 1 1 ... 10 10 10 10 10 10

    units :

    -

    axis :

    Z

    positive :

    down

    standard_name :

    depth

    long_name :

    aquifer_layer

    array([ 1,  1,  1, ..., 10, 10, 10], shape=(1176480,), dtype=int32)

  • z_lower
    (time, zone)
    float64
    nan nan nan nan ... nan nan nan nan

    array([[nan, nan, nan, ..., nan, nan, nan]], shape=(1, 1176480))

  • z_upper
    (time, zone)
    float64
    nan nan nan nan ... nan nan nan nan

    array([[nan, nan, nan, ..., nan, nan, nan]], shape=(1, 1176480))

  • thickness
    (time, zone)
    float64
    nan nan nan nan ... nan nan nan nan

    array([[nan, nan, nan, ..., nan, nan, nan]], shape=(1, 1176480))

  • depth
    (time, zone)
    float64
    nan nan nan nan ... nan nan nan nan

    array([[nan, nan, nan, ..., nan, nan, nan]], shape=(1, 1176480))

• Indexes: (2)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['1850-01-01'], dtype='datetime64[ns]', name='time', freq=None))

  • zone
    PandasIndex

    PandasIndex(Index([      1,       2,       3,       4,       5,       6,       7,       8,
                 9,      10,
           ...
           1176471, 1176472, 1176473, 1176474, 1176475, 1176476, 1176477, 1176478,
           1176479, 1176480],
          dtype='int32', name='zone', length=1176480))

• Attributes: (17)

  conventions :

  CF-1.10

  title :

  marthe_grid_version :

  9.0

  original_dimensions :

  x,y,z [grids]: 387 304 10

  lon_resolution :

  500.0

  lat_resolution :

  500.0

  scale_factor :

  1.0

  nested_grid :

  False

  extend :

  xymin : 543750.0 2526350.0; xymax: 736750.0 2677850.0

  period :

  1850-1850

  frequency :

  nan day(s)

  creation_date :

  Created on 2025-09-02T09:31:50Z UTC

  institution :

  BRGM, French Geological Survey, Orléans, France

  comment :

  Hydrogeological model created with MARTHE code (Thiery, D. 2020. Guidelines for MARTHE v7.8 computer codefor hydro-systems modelling. report BRGM/RP-69660-FR).

  resolution_units :

  m

  projection :

  epsg:27572

  domain :

  FR-France

Interpolation, regridding, transformations

Regrid (coarse/refine resolution)

# Tutorial in construction

Projection transformation