1. Module SIG

Les opérations SIG font partie intégrante des processus hydrologiques. Cette page montre comment utiliser xhydro pour effectuer des manipulations SIG telles que la délimitation des limites des bassins versants et l’extraction de variables physiographiques, climatologiques et géographiques à l’échelle du bassin versant.

[1]:
import leafmap
import pandas as pd
import xarray as xr
import xclim
import xdatasets as xd
from IPython.core.display import HTML, display

import xhydro.gis as xhgis
from xhydro.indicators import get_yearly_op

1.1. Délimitation du bassin versant

Actuellement, la délimitation des bassins versants utilise HydroBASINS (hybas_na_lev01-12_v1c) et peut fonctionner n’importe où en Amérique du Nord. Le processus consiste à évaluer tous les sous-bassins en amont à partir d’un exutoire spécifié et à les consolider en un bassin versant unifié. La librairie leafmap est utilisée pour générer des cartes interactives. Cette carte sert à sélectionner les exutoires ou à visualiser les limites des bassins versants qui en résultent. Bien que l’utilisation de la carte ne soit pas indispensable pour effectuer les calculs, elle s’avère utile à des fins de visualisation.

[2]:
m = leafmap.Map(center=(48.63, -74.71), zoom=5, basemap="USGS Hydrography")
m
[2]:

1.1.1. a) À partir d’une liste de coordonnées

Dans ce scénario, nous sélectionnons deux points d’écoulement, chacun représentant respectivement l’exutoire des bassins versants du Lac Saint-Jean et de la rivière des Outaouais.

[3]:
lng_lat = [
    (-69.81971, 48.14720),  # Lac Saint-Jean watershed
    (-74.393438, 45.572442),  # Ottawa river watershed
]

1.1.2. b) À partir de marqueurs sur une carte

Au lieu d’utiliser une liste, une approche plus interactive consiste à sélectionner directement les points à partir de la carte m existante. L’image suivante illustre le processus de sélection des points d’écoulement en faisant glisser des marqueurs vers les emplacements souhaités sur la carte.

test

The next cell is only useful for the documentation as it simulates a user selecting an outlet from the map m. You should instead remove this code and interact with the map in object m as shown above by positioning markers at sites of interest

[4]:
m.draw_features = [
    {
        "type": "Feature",
        "properties": {},
        "geometry": {"type": "Point", "coordinates": [-73.118597, 46.042467]},
    }
]  # Richelieu watershed

Après avoir sélectionné les points en utilisant l’approche a) ou b), ou une combinaison des deux, nous pouvons lancer le calcul de délimitation du bassin versant.

[5]:
gdf = xhgis.watershed_delineation(coordinates=lng_lat, map=m)
gdf
[5]:
HYBAS_ID Upstream Area (sq. km). geometry category color
0 7120034330 87595.8 POLYGON ((-74.37864 48.88141, -74.37452 48.886... 1 #ffffd9
1 7120034330 87595.8 POLYGON ((-74.33808 49.28658, -74.34388 49.288... 1 #ffffd9
2 7120034330 87595.8 POLYGON ((-74.38747 49.38815, -74.38578 49.386... 1 #ffffd9
3 7120398781 144026.8 POLYGON ((-80.07991 46.7786, -80.08529 46.7827... 2 #c7e9b4
4 7120398781 144026.8 POLYGON ((-78.59583 45.475, -78.59675 45.46602... 2 #c7e9b4
5 7120382860 23717.7 POLYGON ((-73.77437 43.36757, -73.77557 43.388... 1 #ffffd9

Les résultats sont stockés dans un objet GeoPandas gpd.GeoDataFrame (gdf), ce qui nous permet d’enregistrer nos polygones dans divers formats courants tels qu’un ESRI Shapefile ou GeoJSON. Si une carte m est présente, les polygones y seront automatiquement ajoutés. Si vous souhaitez visualiser la carte, tapez simplement m dans la cellule de code pour la restituer. Si l’affichage direct de la carte n’est pas compatible avec votre interpréteur de notebooks, vous pouvez utiliser le code suivant pour extraire le code HTML de la carte et la tracer :

[6]:
m.zoom_to_gdf(gdf)

1.1.3. c) Depuis xdatasets

La délimitation automatique des limites des bassins versants est un outil précieux dans la boîte à outils, mais les utilisateurs sont encouragés à utiliser les limites officielles des bassins versants si elles existent déjà, au lieu d’en créer de nouvelles. Cette fonctionnalité récupère une liste de bassins à partir des ensembles de données pris en charge par xdatasets, et sur demande, xdatasets fournit un gpd.GeoDataFrame contenant les limites précalculées de ces bassins. Les sources de bassin versant suivantes sont disponibles actuellement :

Source

Nom du jeu de données

DEH

deh_polygons

HYDAT

hydat_polygons

HQ

hq_polygons

[7]:
gdf = xd.Query(
    **{
        "datasets": {
            "deh_polygons": {
                "id": ["031*", "0421*"],
                "regulated": ["Natural"],
            }
        }
    }
).data.reset_index()

gdf
[7]:
Station Superficie geometry
0 031501 21.868620 POLYGON ((-72.47379 46.2334, -72.46888 46.2287...
1 031502 15.708960 POLYGON ((-72.50126 46.21216, -72.50086 46.213...
2 042103 579.479614 POLYGON ((-78.49014 46.64514, -78.4901 46.6450...

1.2. Extraire les propriétés du bassin versant

Après avoir obtenu les limites de nos bassins versants, nous pouvons extraire des propriétés telles que des informations géographiques, la classification de l’utilisation des terres et des données climatologiques des bassins versants délimités.

1.2.1. a) Propriétés géographiques des bassins versants

Dans un premier temps, nous extrayons les propriétés géographiques du bassin versant, notamment le périmètre, la superficie totale, le coefficient de Gravelius et le centroïde du bassin. Il est important de noter que cette fonction renvoie toutes les colonnes présentes dans l’argument gpd.GeoDataFrame fourni.

[8]:
Station Superficie area perimeter gravelius centroid
0 031501 21.868620 2.186862e+07 27186.996845 1.640007 (-72.48631199105834, 46.22277542928622)
1 031502 15.708960 1.570896e+07 20263.293021 1.442220 (-72.47966677792694, 46.21359517038631)
2 042103 579.479614 5.794796e+08 283765.058390 3.325331 (-78.37036445281987, 46.48287117609677)

Pour plus de commodité, nous pouvons également récupérer les mêmes résultats sous la forme d’un xarray.Dataset :

[9]:
xhgis.watershed_properties(
    gdf[["Station", "geometry"]], unique_id="Station", output_format="xarray"
)
[9]:
<xarray.Dataset> Size: 120B
Dimensions:    (Station: 3)
Coordinates:
  * Station    (Station) object 24B '031501' '031502' '042103'
Data variables:
    area       (Station) float64 24B 2.187e+07 1.571e+07 5.795e+08
    perimeter  (Station) float64 24B 2.719e+04 2.026e+04 2.838e+05
    gravelius  (Station) float64 24B 1.64 1.442 3.325
    centroid   (Station) object 24B (-72.48631199105834, 46.22277542928622) ....

We can also extract the surface properties for the same gpd.GeoDataFrame :

[10]:
elevation slope aspect proj:shape band platform epsg time gsd proj:epsg spatial_ref
geometry
0 33.573009 0.324613 239.025970 {1200} data TanDEM-X 4326 2021-04-22 90 4326 0
1 51.393295 0.518484 242.431335 {1200} data TanDEM-X 4326 2021-04-22 90 4326 0
2 358.549866 2.500644 178.557648 {1200} data TanDEM-X 4326 2021-04-22 90 4326 0

Again, for convenience, we can output the results in xarray.Dataset format :

[11]:
xhgis.surface_properties(
    gdf[["Station", "geometry"]], output_format="xarray", unique_id="Station"
)
[11]:
<xarray.Dataset> Size: 180B
Dimensions:      (Station: 3)
Coordinates:
    proj:shape   object 8B {1200}
    band         <U4 16B 'data'
    platform     <U8 32B 'TanDEM-X'
    epsg         int64 8B 4326
    time         datetime64[ns] 8B 2021-04-22
    gsd          int64 8B 90
    proj:epsg    int64 8B 4326
    spatial_ref  int64 8B 0
    geometry     (Station) object 24B POLYGON ((-306224.9316606918 257197.438...
  * Station      (Station) object 24B '031501' '031502' '042103'
Data variables:
    elevation    (Station) float32 12B 33.57 51.39 358.5
    slope        (Station) float32 12B 0.3246 0.5185 2.501
    aspect       (Station) float32 12B 239.0 242.4 178.6
Attributes:
    spec:        RasterSpec(epsg=4326, bounds=(-79.00083333333333, 46.0, -72....
    resolution:  0.0008333333333333334
    _FillValue:  nan

1.2.2. b) Classification de l’utilisation des terres

La classification de l’utilisation des terres est alimentée par le catalogue STAC de Planetary Computer. Il utilise par défaut l’ensemble de données « 10m Annual Land Use Land Cover » (io-lulc-9-class), mais d’autres collections peuvent être spécifiées en utilisant l’argument collection.

[12]:
df = xhgis.land_use_classification(gdf, unique_id="Station")
df
[12]:
pct_built_area pct_crops pct_trees pct_rangeland pct_water pct_snow/ice pct_bare_ground pct_flooded_vegetation
Station
031501 0.015321 0.724102 0.255548 0.005029 0.00000 0.000000e+00 0.000000 0.000000
031502 0.009526 0.670792 0.312680 0.007002 0.00000 0.000000e+00 0.000000 0.000000
042103 0.000013 0.000000 0.890441 0.024052 0.08537 3.444143e-07 0.000011 0.000113
[13]:
ax = xhgis.land_use_plot(gdf, unique_id="Station", idx=2)
../_images/notebooks_gis_29_0.png

1.2.3. c) Indicateurs climatiques

L’étape d’extraction des indicateurs climatiques est la plus complexe. En effet, pour ce faire, l’accès à un jeu de données météorologiques pour les différents bassins versants au sein de notre objet gdf est requis. Heureusement, xdatasets facilite précisément de telles opérations. En effet, xdatasets permet d’extraire d’un jeu de données sur grille tous les pixels contenus dans un bassin versant tout en respectant la proportion du bassin versant coupant chaque pixel. Par la suite, la fonction get_yearly_op, construite sur la bibliothèque xclim , offre une flexibilité impressionnante dans la définition d’indicateurs adaptés aux besoins de l’utilisateur.

Pour lancer le processus, nous utilisons les données de réanalyse ERA5-Land couvrant la période de 1981 à 2010 comme jeu de données climatologiques.

[14]:
datasets = {
    "era5_land_reanalysis": {"variables": ["t2m", "tp", "sd"]},
}
space = {
    "clip": "polygon",  # bbox, point or polygon
    "averaging": True,
    "geometry": gdf.iloc[0:2],  # select the 2 first polygons
    "unique_id": "Station",
}
time = {
    "start": "1981-01-01",
    "end": "2010-12-31",
    "timezone": "America/Montreal",
}

ds = xd.Query(datasets=datasets, space=space, time=time).data.squeeze()

Étant donné que les prochaines étapes utilisent xclim sous le capot, le jeu de données doit être conforme aux normes CF. Au minimum, le xarray.DataArray utilisé doit suivre ces principes :

  • Le jeu de données a besoin d’une dimension temporelle, généralement à une fréquence quotidienne sans pas de temps manquants (les NaN sont pris en charge). Si vos données diffèrent de celles-ci, vous devrez être très prudent sur les résultats fournis.

  • S’il existe une dimension spatiale, telle que ``Station``dans l’exemple ci-dessous, elle a besoin d’un attribut cf_role avec timeseries_id comme valeur.

  • La variable aura au moins besoin d’un attribut units, bien que d’autres attributs tels que long_name et cell_methods soient également attendus par xclim et des avertissements seront générés s’ils sont manquants .

  • Bien que cela ne soit pas nécessaire pour get_yearly_op, les noms de variables doivent être parmi ceux pris en charge ici pour une compatibilité maximale.

Le code suivant ajoute les attributs manquants :

[15]:
ds = ds.rename({"t2m": "tas", "tp": "pr", "sd": "snd"})
ds["tas"] = xclim.core.units.convert_units_to(
    ds["tas"], "degC"
)  # Convert from degK to degC
ds["tas"].attrs.update({"cell_methods": "time: mean"})

ds["pr"].attrs.update({"units": "m d-1", "cell_methods": "time: mean within days"})
ds["pr"] = xclim.core.units.convert_units_to(
    ds["pr"], "mm d-1"
)  # Convert from m/d to mm/d

ds["snd"].attrs.update({"units": "m", "cell_methods": "time: mean within days"})
ds["snd"] = xclim.core.units.convert_units_to(ds["snd"], "mm")  # Convert from m to mm
ds
[15]:
<xarray.Dataset> Size: 15MB
Dimensions:  (Station: 2, time: 262968)
Coordinates:
  * time     (time) datetime64[ns] 2MB 1981-01-01 ... 2010-12-31T23:00:00
  * Station  (Station) object 16B '031501' '031502'
    source   <U20 80B 'era5_land_reanalysis'
Data variables:
    tas      (Station, time) float64 4MB -23.29 -23.28 -23.49 ... 2.734 2.454
    pr       (Station, time) float64 4MB 0.0 0.0 0.0 ... 0.01582 0.3079 1.186
    snd      (Station, time) float64 4MB 55.07 55.07 55.07 ... 35.57 35.27 35.01

Dans la deuxième étape, nous pouvons définir les saisons à l’aide d’indexeurs compatibles avec xclim.core.calendar.select_time. Il existe actuellement quatre types d’indexeurs acceptés :

  • month, suivi d’une séquence de numéros de mois.

  • season, suivi d’un ou plusieurs de 'DJF', 'MAM', 'JJA' et 'SON'.

  • doy_bounds, suivi d’une séquence représentant les limites inclusives de la période à considérer ('start', 'end').

  • date_bounds, qui est le même que ci-dessus, mais en utilisant un format mois-jour ('%m-%d').

Ensuite, nous précisons les opérations que nous avons l’intention de calculer pour chaque variable. Les opérations prises en charge incluent \max\, \min\, \mean\ et \sum\.

[16]:
timeargs = {
    "01": {"month": [1]},
    "02": {"month": [2]},
    "03": {"month": [3]},
    "04": {"month": [4]},
    "05": {"month": [5]},
    "06": {"month": [6]},
    "07": {"month": [7]},
    "08": {"month": [8]},
    "09": {"month": [9]},
    "10": {"month": [10]},
    "11": {"month": [11]},
    "12": {"month": [12]},
    "spring": {"date_bounds": ["02-11", "06-19"]},
    "summer_fall": {"date_bounds": ["06-20", "11-19"]},
    "year": {"date_bounds": ["01-01", "12-31"]},
}

operations = {
    "tas": ["max", "mean", "min"],
    "pr": ["sum"],
    "snd": ["mean"],
}

La combinaison des timeargs et des operations à travers le produit cartésien produit rapidement une vaste gamme d’indicateurs climatiques.

[17]:
ds_climatology = xr.merge(
    [
        get_yearly_op(ds, input_var=variable, op=op, timeargs=timeargs)
        for (variable, ops) in operations.items()
        for op in ops
    ]
)
ds_climatology
[17]:
<xarray.Dataset> Size: 36kB
Dimensions:               (Station: 2, time: 30)
Coordinates:
  * Station               (Station) object 16B '031501' '031502'
    source                <U20 80B 'era5_land_reanalysis'
  * time                  (time) datetime64[ns] 240B 1981-01-01 ... 2010-01-01
Data variables: (12/75)
    tas_max_01            (Station, time) float64 480B 5.1 3.616 ... 9.685
    tas_max_02            (Station, time) float64 480B 13.23 1.781 ... 5.404
    tas_max_03            (Station, time) float64 480B 11.61 10.72 ... 11.53
    tas_max_04            (Station, time) float64 480B 19.4 20.66 ... 25.72 25.3
    tas_max_05            (Station, time) float64 480B 24.39 28.17 ... 31.67
    tas_max_06            (Station, time) float64 480B 29.32 28.33 ... 28.63
    ...                    ...
    snd_mean_10           (Station, time) float64 480B 0.3611 ... 0.0389
    snd_mean_11           (Station, time) float64 480B 4.29 2.369 ... 2.396
    snd_mean_12           (Station, time) float64 480B 44.07 5.757 ... 29.87
    snd_mean_spring       (Station, time) float64 480B 6.949 92.08 ... 16.39
    snd_mean_summer_fall  (Station, time) float64 480B 0.1345 0.1227 ... 0.03939
    snd_mean_year         (Station, time) float64 480B 14.36 47.45 ... 15.62
Attributes:
    cat:variable:          ('tas_max_01',)
    cat:xrfreq:            YS-JAN
    cat:frequency:         yr
    cat:processing_level:  indicators
    cat:id:                

Les mêmes données peuvent également être visualisées sous la forme d’un pd.DataFrame :

[18]:
pd.set_option("display.max_rows", 100)
ds_climatology.mean("time").to_dataframe().T
[18]:
Station 031501 031502
source era5_land_reanalysis era5_land_reanalysis
tas_max_01 5.79328 5.83979
tas_max_02 5.709544 5.74949
tas_max_03 12.451726 12.531728
tas_max_04 21.224634 21.254693
tas_max_05 26.618743 26.62262
tas_max_06 29.915254 29.914234
tas_max_07 30.538177 30.541674
tas_max_08 29.401676 29.396607
tas_max_09 26.409743 26.420074
tas_max_10 20.433635 20.457859
tas_max_11 14.840299 14.884756
tas_max_12 7.110035 7.155206
tas_max_spring 29.27284 29.27054
tas_max_summer_fall 31.171837 31.172982
tas_max_year 31.373761 31.372966
tas_mean_01 -11.098958 -11.070966
tas_mean_02 -8.575376 -8.555364
tas_mean_03 -3.011887 -2.994205
tas_mean_04 5.168779 5.193515
tas_mean_05 12.849464 12.856895
tas_mean_06 18.246934 18.250927
tas_mean_07 20.60195 20.605273
tas_mean_08 19.556622 19.560292
tas_mean_09 15.100446 15.109986
tas_mean_10 8.060541 8.074052
tas_mean_11 1.330069 1.348995
tas_mean_12 -6.698234 -6.672344
tas_mean_spring 5.031998 5.047103
tas_mean_summer_fall 14.466165 14.474516
tas_mean_year 6.028106 6.042794
tas_min_01 -28.778872 -28.809396
tas_min_02 -25.183904 -25.231202
tas_min_03 -20.468419 -20.494124
tas_min_04 -7.592577 -7.580951
tas_min_05 1.642597 1.65338
tas_min_06 8.113888 8.126952
tas_min_07 12.12059 12.114125
tas_min_08 10.297272 10.30083
tas_min_09 4.89994 4.900005
tas_min_10 -1.616498 -1.610984
tas_min_11 -11.544179 -11.551885
tas_min_12 -24.749929 -24.769652
tas_min_spring -24.09364 -24.113668
tas_min_summer_fall -6.582494 -6.571843
tas_min_year -30.032823 -30.080185
pr_sum_01 79.061732 79.090602
pr_sum_02 66.888931 66.917058
pr_sum_03 76.444733 76.478452
pr_sum_04 93.790183 93.890757
pr_sum_05 99.64739 99.727621
pr_sum_06 110.389847 110.467514
pr_sum_07 118.354004 118.343317
pr_sum_08 108.421897 108.597886
pr_sum_09 105.164435 105.113248
pr_sum_10 104.317676 104.318942
pr_sum_11 105.9594 105.997777
pr_sum_12 92.276419 92.311662
pr_sum_spring 376.994477 377.308115
pr_sum_summer_fall 552.196939 552.293122
pr_sum_year 1160.716646 1161.254838
snd_mean_01 67.3871 66.592959
snd_mean_02 95.084661 93.688919
snd_mean_03 95.804241 93.601279
snd_mean_04 20.813106 19.620342
snd_mean_05 0.02517 0.021946
snd_mean_06 -0.0 -0.0
snd_mean_07 -0.0 -0.0
snd_mean_08 -0.0 -0.0
snd_mean_09 0.000001 0.000001
snd_mean_10 0.218364 0.213499
snd_mean_11 4.655122 4.610496
snd_mean_12 30.770399 30.495989
snd_mean_spring 41.788584 40.769412
snd_mean_summer_fall 0.319896 0.315362
snd_mean_year 25.923512 25.435699