1. Opérations SIG¶

Les opérations de Système d’Information Géographique (SIG) sont cruciales pour l’analyse hydrologique. Cette section illustre comment utiliser la librairie xHydro pour effectuer des tâches SIG clés, y compris la délimitation des limites de bassins versants et l’extraction de variables physiographiques, climatiques et géographiques essentielles à l’échelle du bassin versant.

[1]:

import leafmap
import numpy as np
import pandas as pd
import xarray as xr
import xclim
import xdatasets

import xhydro as xh
import xhydro.gis as xhgis

1.1. Délimitation du bassin versant¶

Actuellement, la délimitation des bassins versants utilise HydroBASINS de HydroSHEDS <https://www.hydrosheds.org/products/hydrobasins>`__ (hybas_na_lev01-12_v1c) et est compatible avec toute localisation en Amérique du Nord. Le processus consiste à identifier tous les sous-bassins en amont 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, permettant de sélectionner les sorties ou de visualiser les limites du bassin versant résultant. Bien que l’utilisation de la carte ne soit pas obligatoire pour effectuer les calculs, elle améliore considérablement la visualisation et la compréhension du processus de délimitation du bassin versant.

[2]:

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

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

Une première option consiste à fournir une liste de coordonnées. Dans l’exemple ci-dessous, nous sélectionnons deux points, chacun représentant l’exutoire pour les bassins versants des rivières Montmorency et Beaurivage, respectivement.

[3]:

lng_lat = [
    (-71.15069, 46.89431),  # Montmorency river watershed
    (-71.28512, 46.66141),  # Beaurivage river watershed
]

1.1.2. b) À partir de marqueurs sur une carte¶

Plutôt que d’utiliser une liste, une approche plus interactive permet de sélectionner directement les sorties à partir de la carte existante m, en utilisant le bouton Draw a marker situé à gauche de la carte. L’image ci-dessous montre comment sélectionner les points en faisant glisser des marqueurs vers les emplacements souhaités sur la carte.

Après avoir sélectionné des points en utilisant l’approche a) ou b), ou une combinaison des deux, nous pouvons initier le calcul de délimitation du bassin versant. Cela se fait en utilisant la fonction xhydro.gis.watershed_delineation.

[4]:

help(xhgis.watershed_delineation)

Help on function watershed_delineation in module xhydro.gis:

watershed_delineation(*, coordinates: 'list[tuple] | tuple | None' = None, map: 'leafmap.Map | None' = None) -> 'gpd.GeoDataFrame'
    Calculate watershed delineation from pour point.

    Watershed delineation can be computed at any location in North America using HydroBASINS (hybas_na_lev01-12_v1c).
    The process involves assessing all upstream sub-basins from a specified pour point and consolidating them into a
    unified watershed.

    Parameters
    ----------
    coordinates : list of tuple, tuple, optional
        Geographic coordinates (longitude, latitude) for the location where watershed delineation will be conducted.
    map : leafmap.Map, optional
        If the function receives both a map and coordinates as inputs, it will generate and display watershed
        boundaries on the map. Additionally, any markers present on the map will be transformed into
        corresponding watershed boundaries for each marker.

    Returns
    -------
    gpd.GeoDataFrame
        GeoDataFrame containing the watershed boundaries.

    Warnings
    --------
    This function relies on an Amazon S3-hosted dataset to delineate watersheds.

[5]:

gdf = xhgis.watershed_delineation(coordinates=lng_lat, map=m)
gdf

[5]:

	HYBAS_ID	Upstream Area (sq. km).	geometry	category	color
0	7120034480	1157.6	POLYGON ((-71.15208 46.8838, -71.15612 46.8869...	5	#081d58
1	7120365812	693.5	POLYGON ((-71.09758 46.40035, -71.09409 46.403...	1	#ffffd9

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

[6]:

m.zoom_to_gdf(gdf)

1.1.3. c) À partir de xdatasets ¶

Bien que la délimitation automatique des limites de bassins versants soit un outil précieux, les utilisateurs sont encouragés à utiliser les limites officielles des bassins versants lorsqu’elles sont disponibles, plutôt que d’en générer de nouvelles. La librairie xdatasets par exemple, héberge quelques limites officielles qui peuvent être accessibles en utilisant la fonction xdatasets.Query. À ce jour, les sources de bassins versants suivantes sont disponibles :

Source	Nom du jeu de données
(Canada) Archives nationales des données hydrologiques : HYDAT	hydat_polygons
(Québec) Ministère provincial de l’Environnement (MELCCFP/DPEH)	deh_polygons
(Québec) Hydro-Québec	hq_polygons

[7]:

xdatasets.Query(
    **{
        "datasets": {
            "deh_polygons": {
                "id": ["031502", "0421*"],
                "regulated": ["Natural"],
            }
        }
    }
).data.reset_index()

[7]:

	Station	Superficie	geometry
0	031502	15.708960	POLYGON ((-72.50126 46.21216, -72.50086 46.213...
1	042103	579.479614	POLYGON ((-78.49014 46.64514, -78.4901 46.6450...

1.2. Extraction des propriétés du bassin versant¶

Une fois les limites du bassin versant obtenues, nous pouvons extraire des propriétés telles que des informations géographiques, la classification de l’utulisation des sols et des données climatologiques.

1.2.1. a) Propriétés géographiques du bassin versant¶

Tout d’abord, xhydro.gis.watershed_properties peut être utilisé pour extraire les propriétés géographiques du bassin versant, y compris le périmètre, la surface totale, le coefficient de Gravelius et le centre du bassin. Il est important de noter que cette fonction retourne toutes les colonnes présentes dans l’argument gpd.GeoDataFrame fourni.

[8]:

help(xhgis.watershed_properties)

Help on function watershed_properties in module xhydro.gis:

watershed_properties(gdf: 'gpd.GeoDataFrame', unique_id: 'str | None' = None, projected_crs: 'int' = 6622, output_format: 'str' = 'geopandas') -> 'gpd.GeoDataFrame | xr.Dataset'
    Watershed properties extracted from a gpd.GeoDataFrame.

    The calculated properties are :
    - area
    - perimeter
    - gravelius
    - centroid coordinates

    Parameters
    ----------
    gdf : gpd.GeoDataFrame
        GeoDataFrame containing watershed polygons. Must have a defined .crs attribute.
    unique_id : str, optional
        The column name in the GeoDataFrame that serves as a unique identifier.
    projected_crs : int
        The projected coordinate reference system (crs) to utilize for calculations, such as determining watershed area.
    output_format : str
        One of either `xarray` (or `xr.Dataset`) or `geopandas` (or `gpd.GeoDataFrame`).

    Returns
    -------
    gpd.GeoDataFrame or xr.Dataset
        Output dataset containing the watershed properties.

[9]:

xhgis.watershed_properties(gdf)

[9]:

	HYBAS_ID	Upstream Area (sq. km).	category	color	area	perimeter	gravelius	centroid
0	7120034480	1157.6	5	#081d58	1.155329e+09	262827.855265	2.181291	(-71.0919799473914, 47.257838767692064)
1	7120365812	693.5	1	#ffffd9	5.845762e+08	158155.350574	1.845263	(-71.26046365140014, 46.45216102567744)

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

[10]:

xhgis.watershed_properties(gdf, unique_id="HYBAS_ID", output_format="xarray")

1.2.2. b) Propriétés de surface¶

Nous pouvons utiliser xhydro.gis.surface_properties pour extraire les propriétés de surface pour le même gpd.GeoDataFrame, telles que la pente et l’aspect. Par défaut, ces propriétés sont calculées en utilisant le modèle d’élévation numérique GLO-90 de Copernicus au 22 avril 2021. Cependant, la source et la date peuvent être modifiées à l’aide des arguments de la fonction.

[11]:

help(xhgis.surface_properties)

Help on function surface_properties in module xhydro.gis:

surface_properties(gdf: 'gpd.GeoDataFrame', unique_id: 'str | None' = None, projected_crs: 'int' = 6622, output_format: 'str' = 'geopandas', operation: 'str' = 'mean', dataset_date: 'str' = '2021-04-22', collection: 'str' = 'cop-dem-glo-90') -> 'gpd.GeoDataFrame | xr.Dataset'
    Surface properties for watersheds.

    Surface properties are calculated using Copernicus's GLO-90 Digital Elevation Model.
    By default, the dataset has a geographic coordinate system (EPSG: 4326) and this function expects a projected crs for more accurate results.

    The calculated properties are :
    - elevation (meters)
    - slope (degrees)
    - aspect ratio (degrees)

    Parameters
    ----------
    gdf : gpd.GeoDataFrame
        GeoDataFrame containing watershed polygons. Must have a defined .crs attribute.
    unique_id : str, optional
        The column name in the GeoDataFrame that serves as a unique identifier.
    projected_crs : int
        The projected coordinate reference system (crs) to utilize for calculations, such as determining watershed area.
    output_format : str
        One of either `xarray` (or `xr.Dataset`) or `geopandas` (or `gpd.GeoDataFrame`).
    operation : str
        Aggregation statistics such as `mean` or `sum`.
    dataset_date : str
        Date (%Y-%m-%d) for which to select the imagery from the dataset. Date must be available.
    collection : str
        Collection name from the Planetary Computer STAC Catalog. Default is `cop-dem-glo-90`.

    Returns
    -------
    gpd.GeoDataFrame or xr.Dataset
        Output dataset containing the surface properties.

    Warnings
    --------
    This function relies on the Microsoft Planetary Computer's STAC Catalog to retrieve the Digital Elevation Model (DEM) data.

[12]:

xhgis.surface_properties(gdf)

[12]:

	elevation	slope	aspect	platform	proj:epsg	proj:shape	epsg	gsd	time	band	spatial_ref
geometry
0	724.319763	5.564893	186.539612	TanDEM-X	4326	{1200}	4326	90	2021-04-22	data	0
1	224.041672	1.458935	210.367889	TanDEM-X	4326	{1200}	4326	90	2021-04-22	data	0

Encore une fois, pour plus de commodité, nous pouvons afficher les résultats au format xarray.Dataset :

[13]:

xhgis.surface_properties(gdf, unique_id="HYBAS_ID", output_format="xarray")

1.2.4. d) Indicateurs climatiques¶

L’extraction des indicateurs climatiques est l’étape la plus complexe, car elle nécessite l’accès à un jeu de données météorologiques, puis la subdivision et la moyenne des données sur les différents bassins versants contenus dans l’objet GeoDataFrame. Ces étapes dépassent le cadre de xHydro, et les utilisateurs devront recourir à d’autres bibliothèques pour cette tâche. Une approche, décrite dans l”exemple de cas d’utilisation, implique l’utilisation d’une combinaison de xscen et xESMF. Une autre approche, démontrée ici, utilise xdatasets. En effet, xdatasets permet l’extraction de tous les pixels d’un jeu de données à l’intérieur d’un bassin versant tout en tenant compte du poids de chaque pixel intersectant ce dernier.

Pour cet exemple, nous utiliserons les données de la réanalyse ERA5-Land couvrant la période de 1981 à 2010.

[17]:

datasets = {
    "era5_land_reanalysis": {"variables": ["t2m", "tp"]},
}
space = {
    "clip": "polygon",  # bbox, point or polygon
    "averaging": True,
    "geometry": gdf,
    "unique_id": "HYBAS_ID",
}
time = {
    "start": "1981-01-01",
    "end": "1982-12-31",
    "timezone": "America/Montreal",
    "timestep": "D",
    "aggregation": {"tp": np.nansum, "t2m": np.nanmean},
}

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

Temporal operations: processing tp with era5_land_reanalysis: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:00<00:00,  5.91it/s]

Parce que les étapes suivantes utilisent xclim en arrière-plan, le jeu de données doit être conforme aux normes CF. Au minimum, le xarray.DataArray utilisé doit respecter les principes suivants :

Le jeu de données doit inclure une dimension time, généralement à une fréquence quotidienne, sans pas de temps manquant (bien que les NaNs soient supportés). Si vos données dévient de ce format, une prudence particulière doit être prise lors de l’interprétation des résultats.
S’il y a une dimension spatiale 1D, comme HYBAS_ID dans l’exemple ci-dessous, elle doit avoir un attribut cf_role avec la valeur timeseries_id.
La variable doit inclure au moins un attribut units. Bien que non obligatoire, d’autres attributs tels que long_name et cell_methods sont attendus par xclim, et des avertissements seront émis si ils sont manquants.
Les noms des variables doivent correspondre à ceux supportés par xclim.

Le code suivant formatera les données ERA5-Land obtenues depuis xdatasets pour ajouter les métadonnées manquantes et changer les noms des variables et les unités :

[18]:

ds = ds.rename({"t2m_nanmean": "tas", "tp_nansum": "pr"})
ds["tas"] = xclim.core.units.convert_units_to(ds["tas"], "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")
ds

Les indicateurs climatiques peuvent être calculés de différentes manières, et pour la plupart des tâches simples, utiliser directement xclim est toujours une option viable. xclim propose une liste étendue d” indicateurs disponibles dans un cadre convivial et flexible.

[19]:

xclim.atmos.tg_max(ds.tas)

Pour des tâches plus complexes, ou lorsqu’il s’agit de calculer plusieurs indicateurs simultanément, xhydro.indicators.compute_indicators est souvent la méthode privilégiée. Cette fonction permet aux utilisateurs de créer et de passer plusieurs indicateurs à la fois, soit en fournissant une liste d’indicateurs personnalisés créés avec xclim.core.indicator.Indicator.from_dict (voir la boîte INFO ci-dessous), soit en référant un fichier YAML (voir la documentation de xscen).

INFO

Les indicateurs personnalisés dans xHydro sont construits en suivant le formatage YAML requis par xclim.

Un indicateur personnalisé construit avec xclim.core.indicator.Indicator.from_dict nécessitera ces éléments :

« data » : Un dictionnaire avec les informations suivantes :
- « base » : L“« ID YAML » obtenu depuis cette page.
- « input » : Un dictionnaire reliant l’entrée xclim par défaut au nom de votre variable. Nécessaire uniquement si elle diffère. Dans le lien ci-dessus, ce sont les chaînes suivant « Uses: ».
- « parameters » : Un dictionnaire contenant tous les autres paramètres pour un indicateur donné. Dans le lien ci-dessus, la manière la plus simple d’y accéder est de cliquer sur le lien dans le coin supérieur droit de la boîte décrivant un indicateur donné.
- Des entrées supplémentaires peuvent être utilisées ici, comme décrit dans la documentation xclim sous « identifier ».
« identifier » : Un nom personnalisé pour votre indicateur. Ce sera le nom retourné dans les résultats.
« module »: Nécessaire, mais peut être n’importe quoi. Pour éviter un écrasement accidentel des indicateurs xclim, il est préférable d’utiliser quelque chose de différent de : [« atmos », « land », « generic »].

[20]:

help(xh.indicators.compute_indicators)

Help on function compute_indicators in module xscen.indicators:

compute_indicators(ds: xarray.core.dataset.Dataset, indicators: str | os.PathLike | collections.abc.Sequence[xclim.core.indicator.Indicator] | collections.abc.Sequence[tuple[str, xclim.core.indicator.Indicator]] | module, *, periods: list[str] | list[list[str]] | None = None, restrict_years: bool = True, to_level: str | None = 'indicators', rechunk_input: bool = False) -> dict
    Calculate variables and indicators based on a YAML call to xclim.

    The function cuts the output to be the same years as the inputs.
    Hence, if an indicator creates a timestep outside the original year range (e.g. the first DJF for QS-DEC),
    it will not appear in the output.

    Parameters
    ----------
    ds : xr.Dataset
        Dataset to use for the indicators.
    indicators : str | os.PathLike | Sequence[Indicator] | Sequence[tuple[str, Indicator]] | ModuleType
        Path to a YAML file that instructs on how to calculate missing variables.
        Can also be only the "stem", if translations and custom indices are implemented.
        Can be the indicator module directly, or a sequence of indicators or a sequence of
        tuples (indicator name, indicator) as returned by `iter_indicators()`.
    periods : list of str or list of lists of str, optional
        Either [start, end] or list of [start, end] of continuous periods over which to compute the indicators.
        This is needed when the time axis of ds contains some jumps in time.
        If None, the dataset will be considered continuous.
    restrict_years : bool
        If True, cut the time axis to be within the same years as the input.
        This is mostly useful for frequencies that do not start in January, such as QS-DEC.
        In that instance, `xclim` would start on previous_year-12-01 (DJF), with a NaN.
        `restrict_years` will cut that first timestep. This should have no effect on YS and MS indicators.
    to_level : str, optional
        The processing level to assign to the output.
        If None, the processing level of the inputs is preserved.
    rechunk_input : bool
        If True, the dataset is rechunked with :py:func:`flox.xarray.rechunk_for_blockwise`
        according to the resampling frequency of the indicator. Each rechunking is done
        once per frequency with :py:func:`xscen.utils.rechunk_for_resample`.

    Returns
    -------
    dict
        Dictionary (keys = timedeltas) with indicators separated by temporal resolution.

    See Also
    --------
    xclim.indicators, xclim.core.indicator.build_indicator_module_from_yaml

[21]:

indicators = [
    # Maximum summer temperature
    xclim.core.indicator.Indicator.from_dict(
        data={
            "base": "tg_max",
            "parameters": {"indexer": {"month": [6, 7, 8]}},
        },
        identifier="tg_max_summer",
        module="example",
    ),
    # Number of days with more than 5 mm of precipitation
    xclim.core.indicator.Indicator.from_dict(
        data={
            "base": "wetdays",
            "parameters": {"thresh": "5 mm d-1"},
        },
        identifier="wet5mm",
        module="example",
    ),
]

xh.indicators.compute_indicators(ds, indicators)["YS-JAN"]

Enfin, xhydro.indicators.get_yearly_op, également basé sur la librairie xclim, fournit une méthode flexible pour obtenir les valeurs annuelles de statistiques spécifiques, telles que les maxima annuels ou saisonniers. Cela peut être particulièrement utile, par exemple, lors de l’extraction des données brutes nécessaires pour les analyses fréquentielles.

[22]:

help(xh.indicators.get_yearly_op)

Help on function get_yearly_op in module xhydro.indicators.generic:

get_yearly_op(ds, op, *, input_var: str = 'streamflow', window: int = 1, timeargs: dict | None = None, missing: str = 'skip', missing_options: dict | None = None, interpolate_na: bool = False) -> xarray.core.dataset.Dataset
    Compute yearly operations on a variable.

    Parameters
    ----------
    ds : xr.Dataset
        Dataset containing the variable to compute the operation on.
    op : str
        Operation to compute. One of ["max", "min", "mean", "sum"].
    input_var : str
        Name of the input variable. Defaults to "streamflow".
    window : int
        Size of the rolling window. A "mean" operation is performed on the rolling window before the call to xclim.
        This parameter cannot be used with the "sum" operation.
    timeargs : dict, optional
        Dictionary of time arguments for the operation.
        Keys are the name of the period that will be added to the results (e.g. "winter", "summer", "annual").
        Values are up to two dictionaries, with both being optional.
        The first is {'freq': str}, where str is a frequency supported by xarray (e.g. "YS", "YS-JAN", "YS-DEC").
        It needs to be a yearly frequency. Defaults to "YS-JAN".
        The second is an indexer as supported by :py:func:`xclim.core.calendar.select_time`.
        Defaults to {}, which means the whole year.
        See :py:func:`xclim.core.calendar.select_time` for more information.
        Examples: {"winter": {"freq": "YS-DEC", "date_bounds": ["12-01", "02-28"]}}, {"jan": {"freq": "YS", "month": 1}}, {"annual": {}}.
    missing : str
        How to handle missing values. One of "skip", "any", "at_least_n", "pct", "wmo".
        See :py:func:`xclim.core.missing` for more information.
    missing_options : dict, optional
        Dictionary of options for the missing values' method. See :py:func:`xclim.core.missing` for more information.
    interpolate_na : bool
        Whether to interpolate missing values before computing the operation. Only used with the "sum" operation.
        Defaults to False.

    Returns
    -------
    xr.Dataset
        Dataset containing the computed operations, with one variable per indexer.
        The name of the variable follows the pattern `{input_var}{window}_{op}_{indexer}`.

    Notes
    -----
    If you want to perform a frequency analysis on a frequency that is finer than annual, simply use multiple timeargs
    (e.g. 1 per month) to create multiple distinct variables.

L’argument timeargs repose sur des indexeurs compatibles avec xclim.core.calendar.select_time. Quatre types d’indexeurs sont actuellement acceptés :

month, suivi d’une séquence de numéros de mois.
season, suivi de un ou plusieurs parmi 'DJF', 'MAM', 'JJA', et 'SON'.
doy_bounds, suivi d’une séquence représentant les bornes inclusives de la période à considérer ('début', 'fin').
date_bounds, similaire à doy_bounds, mais utilisant un format mois-jour ('%m-%d').

Ensuite, nous spécifions les opérations à effectuer pour chaque variable. Les opérations prises en charge incluent "max", "min", "mean", et "sum".

[23]:

timeargs = {
    "january": {"month": [1]},
    "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"],
    "pr": ["sum"],
}

La combinaison de timeargs et operations par le produit cartésien permet la génération efficace d’un ensemble complet d’indicateurs climatiques.

[24]:

ds_climatology = xr.merge(
    [
        xh.indicators.get_yearly_op(ds, input_var=variable, op=op, timeargs=timeargs)
        for (variable, ops) in operations.items()
        for op in ops
    ]
)
ds_climatology

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

[25]:

pd.set_option("display.max_rows", 100)
ds_climatology.mean("time").to_dataframe().T

[25]:

HYBAS_ID	7120034480	7120365812
spatial_agg	polygon	polygon
timestep	D	D
source	era5_land_reanalysis	era5_land_reanalysis
tas_max_january	-5.979976	-1.652047
tas_max_spring	17.686927	20.677712
tas_max_summer_fall	20.757984	24.204149
tas_max_year	20.757984	24.204149
tas_mean_january	-18.967689	-16.441736
tas_mean_spring	0.424023	4.184447
tas_mean_summer_fall	9.293055	12.396911
tas_mean_year	1.134716	4.487463
pr_sum_january	61.612786	57.209695
pr_sum_spring	467.706516	393.650887
pr_sum_summer_fall	614.891367	639.558307
pr_sum_year	1351.15288	1290.998159