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_gig.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) int64 8B 0
  * zone     (zone) int32 405kB 47402 47403 47763 ... 2381150 2381151 2381457
Data variables:
    permeab  (time, zone) float64 810kB 1e-05 1e-05 1e-05 ... 0.0001 0.0001
    z        (zone) int64 810kB 1 1 1 1 1 1 1 1 1 1 1 ... 9 9 9 9 9 9 9 9 9 9 9
    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
    izone    (zone) int32 405kB 1 2 3 4 5 ... 101205 101206 101207 101208 101209
Attributes: (12/15)
    conventions:          CF-1.10
    title:                
    marthe_grid_version:  9.0
    original_dimensions:  x,y,z [grids]: 387 304 10; 307 442 10
    crs:                  {'crs_wkt': 'PROJCRS["NTF (Paris) / Lambert zone II...
    lon_resolution:       100.0, 500.0
    ...                   ...
    nested_grid:          True
    extend:               xymin : 543750.0 2526350.0; xymax: 736750.0 2677850.0
    creation_date:        Created on 2026-05-29T12:15:40Z UTC
    comment:              Hydrogeological model created with MARTHE code (Thi...
    domain:               FR-France
    institution:          BRGM, French Geological Survey, Orléans, France

xarray.Dataset

• Dimensions:
  • time: 1
  • zone: 101209
• Coordinates: (2)
  • time
    (time)
    int64
    0

    array([0])

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

    compress :

    z y x

    array([  47402,   47403,   47763, ..., 2381150, 2381151, 2381457],
          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

    mart_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))

  • z
    (zone)
    int64
    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,))

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

    axis :

    X

    long_name :

    Easting

    standard_name :

    projection_x_coordinate

    units :

    metre

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

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

    axis :

    Y

    long_name :

    Northing

    standard_name :

    projection_y_coordinate

    units :

    metre

    array([2.61685e+06, 2.61685e+06, 2.61635e+06, ..., 2.59235e+06,
           2.59235e+06, 2.59225e+06], 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)

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

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

• Attributes: (15)

  conventions :

  CF-1.10

  title :

  marthe_grid_version :

  9.0

  original_dimensions :

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

  crs :

  {'crs_wkt': 'PROJCRS["NTF (Paris) / Lambert zone II",BASEGEOGCRS["NTF (Paris)",DATUM["Nouvelle Triangulation Francaise (Paris)",ELLIPSOID["Clarke 1880 (IGN)",6378249.2,293.466021293627,LENGTHUNIT["metre",1]]],PRIMEM["Paris",2.5969213,ANGLEUNIT["grad",0.0157079632679489]],ID["EPSG",4807]],CONVERSION["Lambert zone II",METHOD["Lambert Conic Conformal (1SP)",ID["EPSG",9801]],PARAMETER["Latitude of natural origin",52,ANGLEUNIT["grad",0.0157079632679489],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["grad",0.0157079632679489],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.99987742,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",600000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",2200000,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting (X)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["northing (Y)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Engineering survey, topographic mapping."],AREA["France mainland onshore between 50.5 grads and 53.5 grads North (45°27\'N to 48°09\'N). Also used over all onshore mainland France."],BBOX[42.33,-4.87,51.14,8.23]],ID["EPSG",27572]]', 'semi_major_axis': 6378249.2, 'semi_minor_axis': 6356515.0, 'inverse_flattening': 293.4660212936269, 'reference_ellipsoid_name': 'Clarke 1880 (IGN)', 'longitude_of_prime_meridian': 2.5969213, 'prime_meridian_name': 'Paris', 'geographic_crs_name': 'NTF (Paris)', 'horizontal_datum_name': 'Nouvelle Triangulation Francaise (Paris)', 'projected_crs_name': 'NTF (Paris) / Lambert zone II', 'grid_mapping_name': 'lambert_conformal_conic', 'standard_parallel': 52.0, 'longitude_of_central_meridian': 0.0, 'false_easting': 600000.0, 'false_northing': 2200000.0}

  lon_resolution :

  100.0, 500.0

  lat_resolution :

  100.0, 500.0

  resolution_units :

  metre

  scale_factor :

  1.0

  nested_grid :

  True

  extend :

  xymin : 543750.0 2526350.0; xymax: 736750.0 2677850.0

  creation_date :

  Created on 2026-05-29T12:15:40Z UTC

  comment :

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

  domain :

  FR-France

  institution :

  BRGM, French Geological Survey, Orléans, France

Get surface mask

surf = gm.get_surface_layer(grid)
surf

<xarray.Dataset> Size: 1MB
Dimensions:  (zone: 21025)
Coordinates:
  * zone     (zone) int32 84kB 178171 177784 647989 ... 1389346 1660427 1660120
Data variables:
    time     (zone) int64 168kB 0 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0 0
    permeab  (zone) float64 168kB 1e-06 1e-06 0.003068 ... 1e-06 2e-05 2e-05
    z        (zone) int64 168kB 2 2 6 6 8 8 8 8 5 2 6 ... 2 2 2 2 2 2 2 2 2 4 4
    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
    izone    (zone) int32 84kB 345 331 2139 2111 ... 13132 13062 19985 19965

xarray.Dataset

• Dimensions:
  • zone: 21025
• Coordinates: (1)
  • zone
    (zone)
    int32
    178171 177784 ... 1660427 1660120

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

• Data variables: (8)
  • time
    (zone)
    int64
    0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0

    array([0, 0, 0, ..., 0, 0, 0], shape=(21025,))

  • 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,))

  • z
    (zone)
    int64
    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,))

  • 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([2.59985e+06, 2.60035e+06, 2.60085e+06, ..., 2.60605e+06,
           2.60615e+06, 2.60625e+06], 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)

  • izone
    (zone)
    int32
    345 331 2139 ... 13062 19985 19965

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

# 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

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: 80MB
Dimensions:    (time: 1, zone: 1176480)
Coordinates:
  * time       (time) int64 8B 0
  * 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
    z          (zone) int64 9MB 1 1 1 1 1 1 1 1 1 ... 10 10 10 10 10 10 10 10 10
    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_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/15)
    conventions:          CF-1.10
    title:                
    marthe_grid_version:  9.0
    original_dimensions:  x,y,z [grids]: 387 304 10
    crs:                  {'crs_wkt': 'PROJCRS["NTF (Paris) / Lambert zone II...
    lon_resolution:       500.0
    ...                   ...
    nested_grid:          False
    extend:               xymin : 543750.0 2526350.0; xymax: 736750.0 2677850.0
    creation_date:        Created on 2026-05-29T12:15:44Z UTC
    comment:              Hydrogeological model created with MARTHE code (Thi...
    domain:               FR-France
    institution:          BRGM, French Geological Survey, Orléans, France

xarray.Dataset

• Dimensions:
  • time: 1
  • zone: 1176480
• Coordinates: (2)
  • time
    (time)
    int64
    0

    array([0])

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

    compress :

    z y x

    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))

  • z
    (zone)
    int64
    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,))

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

    axis :

    X

    long_name :

    Easting

    standard_name :

    projection_x_coordinate

    units :

    metre

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

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

    axis :

    Y

    long_name :

    Northing

    standard_name :

    projection_y_coordinate

    units :

    metre

    array([2.67785e+06, 2.67785e+06, 2.67785e+06, ..., 2.52635e+06,
           2.52635e+06, 2.52635e+06], 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_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))

• Attributes: (15)

  conventions :

  CF-1.10

  title :

  marthe_grid_version :

  9.0

  original_dimensions :

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

  crs :

  {'crs_wkt': 'PROJCRS["NTF (Paris) / Lambert zone II",BASEGEOGCRS["NTF (Paris)",DATUM["Nouvelle Triangulation Francaise (Paris)",ELLIPSOID["Clarke 1880 (IGN)",6378249.2,293.466021293627,LENGTHUNIT["metre",1]]],PRIMEM["Paris",2.5969213,ANGLEUNIT["grad",0.0157079632679489]],ID["EPSG",4807]],CONVERSION["Lambert zone II",METHOD["Lambert Conic Conformal (1SP)",ID["EPSG",9801]],PARAMETER["Latitude of natural origin",52,ANGLEUNIT["grad",0.0157079632679489],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["grad",0.0157079632679489],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.99987742,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",600000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",2200000,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting (X)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["northing (Y)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Engineering survey, topographic mapping."],AREA["France mainland onshore between 50.5 grads and 53.5 grads North (45°27\'N to 48°09\'N). Also used over all onshore mainland France."],BBOX[42.33,-4.87,51.14,8.23]],ID["EPSG",27572]]', 'semi_major_axis': 6378249.2, 'semi_minor_axis': 6356515.0, 'inverse_flattening': 293.4660212936269, 'reference_ellipsoid_name': 'Clarke 1880 (IGN)', 'longitude_of_prime_meridian': 2.5969213, 'prime_meridian_name': 'Paris', 'geographic_crs_name': 'NTF (Paris)', 'horizontal_datum_name': 'Nouvelle Triangulation Francaise (Paris)', 'projected_crs_name': 'NTF (Paris) / Lambert zone II', 'grid_mapping_name': 'lambert_conformal_conic', 'standard_parallel': 52.0, 'longitude_of_central_meridian': 0.0, 'false_easting': 600000.0, 'false_northing': 2200000.0}

  lon_resolution :

  500.0

  lat_resolution :

  500.0

  resolution_units :

  metre

  scale_factor :

  1.0

  nested_grid :

  False

  extend :

  xymin : 543750.0 2526350.0; xymax: 736750.0 2677850.0

  creation_date :

  Created on 2026-05-29T12:15:44Z UTC

  comment :

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

  domain :

  FR-France

  institution :

  BRGM, French Geological Survey, Orléans, France

Interpolation, regridding, transformations

Regrid (coarse/refine resolution)

# Tutorial in construction

Projection transformation