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) datetime64[ns] 8B 1850-01-01
  * 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
    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) int64 810kB 1 1 1 1 1 1 1 1 1 1 1 ... 9 9 9 9 9 9 9 9 9 9 9
    izone    (zone) int32 405kB 1 2 3 4 5 ... 101205 101206 101207 101208 101209
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
    crs:                  {'crs_wkt': 'PROJCRS["NTF (Paris) / Lambert zone II...
    lon_resolution:       100.0, 500.0
    ...                   ...
    period:               1850-1850
    frequency:            nan day(s)
    creation_date:        Created on 2025-12-09T16:23:56Z 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)
    datetime64[ns]
    1850-01-01

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

  • 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

    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 :

    meters

    axis :

    X

    standard_name :

    projection_x_coordinate

    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 :

    meters

    axis :

    Y

    standard_name :

    projection_y_coordinate

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

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

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

• Indexes: (2)
  • time
    PandasIndex

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

  • zone
    PandasIndex

    PandasIndex(Index([  47402,   47403,   47763,   47764,   47765,   47766,   47767,   47770,
             47789,   47790,
           ...
           2380539, 2380540, 2380541, 2380843, 2380844, 2380845, 2380846, 2381150,
           2381151, 2381457],
          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

  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

  period :

  1850-1850

  frequency :

  nan day(s)

  creation_date :

  Created on 2025-12-09T16:23:56Z 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) 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) int64 168kB 2 2 6 6 8 8 8 8 5 2 6 ... 2 2 2 2 2 2 2 2 2 4 4
    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)
    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)
    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,))

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

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

• Indexes: (1)
  • zone
    PandasIndex

    PandasIndex(Index([ 178171,  177784,  647989,  647602,  882511,  882124,  881737,  881350,
            528019,  179720,
           ...
           1391495, 1391188, 1390881, 1390574, 1390267, 1389960, 1389653, 1389346,
           1660427, 1660120],
          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.3.0/lib/python3.10/site-packages/gridmarthe/operasem/gis.py:71: 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: 80MB
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) int64 9MB 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
    crs:                  {'crs_wkt': 'PROJCRS["NTF (Paris) / Lambert zone II...
    lon_resolution:       500.0
    ...                   ...
    period:               1850-1850
    frequency:            nan day(s)
    creation_date:        Created on 2025-12-09T16:24:01Z 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)
    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

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

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

    units :

    meters

    axis :

    X

    standard_name :

    projection_x_coordinate

    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 :

    meters

    axis :

    Y

    standard_name :

    projection_y_coordinate

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

  • 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

  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

  period :

  1850-1850

  frequency :

  nan day(s)

  creation_date :

  Created on 2025-12-09T16:24:01Z 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