Exemple de cas d’utilisation¶

Cet exemple illustre un cas d’utilisation qui couvre les étapes essentielles de la construction d’un modèle hydrologique et de l’analyse de l’impact des changements climatiques :

Identification du bassin versant et de ses principales caractéristiques
- Bassin versant de la rivière Beaurivage dans le sud du Québec, à l’emplacement de la station 023401.
Collecte des données observées
- Données ERA5-Land et données de débit.
Préparation et calage du modèle hydrologique
- GR4JCN émulé par le cadre hydrologique Raven.
Calcul des indicateurs hydrologiques
- Débit moyen estival
- Débit moyen mensuel
- Débit maximum sur 20 et 100 ans
- Débit estival minimal sur 7 jours de récurrence 2 ans
Évaluation de l’impact des changements climatiques
- Simulations CMIP6 post-traitées du jeu de données ESPO-G6-R2

Identification d’un bassin versant et de ses caractéristiques¶

INFO

Pour plus d’informations sur cette section et les options disponibles, consultez le Notebook GIS.

Cette première étape dépend fortement du modèle hydrologique. Puisque nous utiliserons GR4JCN dans notre exemple, nous devons obtenir la superficie du bassin versant, les coordonnées du centroïde et l’altitude. Nous aurons également besoin de la délimitation du bassin versant pour extraire les données météorologiques. Toutes ces informations peuvent être obtenues via la fonction xhydro.gis.watershed_to_raven_hru, qui fait appel à diverses fonctions de ce module.

[1]:

import warnings
warnings.filterwarnings("ignore", category=UserWarning)
warnings.filterwarnings("ignore", category=FutureWarning)

[2]:

from IPython.display import clear_output

import xhydro.gis as xhgis

clear_output(wait=False)

[3]:

# Watershed delineation
coords = (-71.28878, 46.65692)
gdf = xhgis.watershed_to_raven_hru(coords)
gdf

[3]:

	HRU_ID	area	latitude	longitude	elevation	SubId	DowSubId	geometry
0	7120365812	585.585577	46.452161	-71.260464	222.55365	1	-1	POLYGON ((-71.09758 46.40035, -71.09409 46.403...

Puisque xhgis.watershed_delineation extrait le polygone HydroBASINS le plus proche, le bassin versant pourrait ne pas correspondre exactement aux coordonnées demandées. La station de jaugeage 023401 a une aire de drainage associée de 708 km², ce qui diffère de nos résultats. Le débit devra être ajusté en utilisant un facteur d’échelle.

[4]:

gauge_area = 708
scaling_factor = gdf.iloc[0]["area"] / gauge_area
scaling_factor

[4]:

np.float64(0.8270982719703007)

Collecte des données observées¶

[5]:

import geopandas as gpd
import matplotlib.pyplot as plt
import xarray as xr

# For easy access to the specific streamflow data used here
import xdatasets
import xscen

Données météorologiques¶

INFO

Plusieurs librairies peuvent être utilisées pour effectuer ces étapes. Par souci de simplicité, cet exemple utilisera les modules subset et aggregate de la librairie xscen.

Cet exemple utilisera les données journalières ERA5-Land hébergées sur la plateforme PAVICS.

[6]:

# Extraction of ERA5-Land data
meteo_ref = xr.open_dataset(
    "https://pavics.ouranos.ca/twitcher/ows/proxy/thredds/dodsC/datasets/reanalyses/day_ERA5-Land_NAM.ncml",
    engine="netcdf4",
    chunks={"time": 365, "lon": 50, "lat": 50},
)[["pr", "tasmin", "tasmax"]]
meteo_ref

[6]:

<xarray.Dataset> Size: 454GB
Dimensions:  (time: 27790, lat: 801, lon: 1700)
Coordinates:
  * time     (time) datetime64[ns] 222kB 1950-01-01 1950-01-02 ... 2026-01-31
  * lat      (lat) float32 3kB 10.0 10.1 10.2 10.3 10.4 ... 89.7 89.8 89.9 90.0
  * lon      (lon) float32 7kB -179.9 -179.8 -179.7 -179.6 ... -10.2 -10.1 -10.0
Data variables:
    pr       (time, lat, lon) float32 151GB dask.array<chunksize=(365, 50, 50), meta=np.ndarray>
    tasmin   (time, lat, lon) float32 151GB dask.array<chunksize=(365, 50, 50), meta=np.ndarray>
    tasmax   (time, lat, lon) float32 151GB dask.array<chunksize=(365, 50, 50), meta=np.ndarray>
Attributes: (12/30)
    Conventions:               CF-1.9
    cell_methods:              time: mean (interval: 1 day)
    doi:                       https://doi.org/10.24381/cds.e2161bac
    domain:                    NAM
    frequency:                 day
    history:                   [2022-12-25 09:07:39.901698] Converted variabl...
    ...                        ...
    institute_id:              ECMWF
    dataset_id:                ERA5-Land
    abstract:                  ERA5-Land provides hourly high resolution info...
    dataset_description:       https://www.ecmwf.int/en/era5-land
    attribution:               Contains modified Copernicus Climate Change Se...
    citation:                  Muñoz Sabater, J., (2021): ERA5-Land hourly da...

Ce jeu de données couvre l’ensemble du globe et possède plus de 70 ans de données. La première étape consistera donc à sous-échantillonner les données spatialement et temporellement. Pour le sous-échantillonnage spatial, le GeoDataFrame obtenu précédemment peut être utilisé.

[7]:

meteo_ref = meteo_ref.sel(time=slice("1991", "2020"))  # Temporal subsetting
meteo_ref = xscen.spatial.subset(
    meteo_ref, method="shape", tile_buffer=2, shape=gdf
)  # Spatial subsetting, with a buffer of 2 grid cells
meteo_ref

[7]:

<xarray.Dataset> Size: 9MB
Dimensions:  (time: 10958, lat: 8, lon: 8)
Coordinates:
  * time     (time) datetime64[ns] 88kB 1991-01-01 1991-01-02 ... 2020-12-31
  * lat      (lat) float32 32B 46.1 46.2 46.3 46.4 46.5 46.6 46.7 46.8
  * lon      (lon) float32 32B -71.6 -71.5 -71.4 -71.3 -71.2 -71.1 -71.0 -70.9
Data variables:
    pr       (time, lat, lon) float32 3MB dask.array<chunksize=(355, 8, 8), meta=np.ndarray>
    tasmin   (time, lat, lon) float32 3MB dask.array<chunksize=(355, 8, 8), meta=np.ndarray>
    tasmax   (time, lat, lon) float32 3MB dask.array<chunksize=(355, 8, 8), meta=np.ndarray>
Attributes: (12/30)
    Conventions:               CF-1.9
    cell_methods:              time: mean (interval: 1 day)
    doi:                       https://doi.org/10.24381/cds.e2161bac
    domain:                    NAM
    frequency:                 day
    history:                   [2026-03-31 14:35:41] shape spatial subsetting...
    ...                        ...
    institute_id:              ECMWF
    dataset_id:                ERA5-Land
    abstract:                  ERA5-Land provides hourly high resolution info...
    dataset_description:       https://www.ecmwf.int/en/era5-land
    attribution:               Contains modified Copernicus Climate Change Se...
    citation:                  Muñoz Sabater, J., (2021): ERA5-Land hourly da...

[8]:

ax = plt.subplot(1, 1, 1)
meteo_ref.tasmin.isel(time=0).plot(ax=ax)
gdf.plot(ax=ax)

[8]:

<Axes: title={'center': 'time = 1991-01-01'}, xlabel='longitude [degrees_east]', ylabel='latitude [degrees_north]'>

Raven attend des températures en Celsius et des précipitations en millimètres, mais elles sont actuellement dans un format conforme aux conventions CF, soient en Kelvin et en kg m⁻² s⁻¹, respectivement. La fonction xhydro.modelling.format_input peut être utilisée pour préparer les données pour la modélisation avec Raven. Elle effectue la conversion des unités, le renommage des variables et le formatage des coordonnées afin d’assurer la compatibilité avec RavenPy. Dans le cas de données météorologiques maillées—comme dans cet exemple—xHydro utilise des fonctions disponibles dans RavenPy pour attribuer des poids à chaque cellule de la grille en fonction de la portion qui chevauche le bassin versant. Alternativement, les données peuvent être agrégées manuellement avant d’être transmises au modèle.

Par souci de simplification, l’altitude de la grille sera fixée à 450 m. Le calcul de l’altitude des cellules de grille dans ERA5-Land n’est pas toujours trivial et ne fait pas partie du cadre de cet exemple.

[9]:

from pathlib import Path
import tempfile
notebook_folder = Path(tempfile.TemporaryDirectory().name)

import xhydro as xh

# Add altitude data
meteo_ref = meteo_ref.assign_coords(
    {"elevation": xr.ones_like(meteo_ref.pr.isel(time=0).drop_vars("time")) * 450}
)
meteo_ref["elevation"].attrs = {"units": "m"}

meteo_ref, config_meteo_ref = xh.modelling.format_input(
    meteo_ref, model="GR4JCN", save_as=notebook_folder / "_data" / "meteo.nc"
)
meteo_ref

[9]:

<xarray.Dataset> Size: 9MB
Dimensions:    (time: 10958, latitude: 8, longitude: 8)
Coordinates:
  * time       (time) datetime64[ns] 88kB 1991-01-01 1991-01-02 ... 2020-12-31
  * latitude   (latitude) float32 32B 46.1 46.2 46.3 46.4 46.5 46.6 46.7 46.8
  * longitude  (longitude) float32 32B -71.6 -71.5 -71.4 ... -71.1 -71.0 -70.9
    elevation  (latitude, longitude) float32 256B dask.array<chunksize=(8, 8), meta=np.ndarray>
Data variables:
    pr         (time, latitude, longitude) float32 3MB dask.array<chunksize=(355, 8, 8), meta=np.ndarray>
    tasmin     (time, latitude, longitude) float32 3MB dask.array<chunksize=(355, 8, 8), meta=np.ndarray>
    tasmax     (time, latitude, longitude) float32 3MB dask.array<chunksize=(355, 8, 8), meta=np.ndarray>
Attributes: (12/30)
    Conventions:               CF-1.9
    cell_methods:              time: mean (interval: 1 day)
    doi:                       https://doi.org/10.24381/cds.e2161bac
    domain:                    NAM
    frequency:                 day
    history:                   [2026-03-31 14:35:41] shape spatial subsetting...
    ...                        ...
    institute_id:              ECMWF
    dataset_id:                ERA5-Land
    abstract:                  ERA5-Land provides hourly high resolution info...
    dataset_description:       https://www.ecmwf.int/en/era5-land
    attribution:               Contains modified Copernicus Climate Change Se...
    citation:                  Muñoz Sabater, J., (2021): ERA5-Land hourly da...

Cette fonction retourne également de l’information qui sera utilisée plus tard pour instantier le modèle hydrologique:

[10]:

config_meteo_ref

[10]:

{'data_type': ['TEMP_MAX', 'TEMP_MIN', 'PRECIP'],
 'alt_names_meteo': {'TEMP_MAX': 'tasmax',
  'TEMP_MIN': 'tasmin',
  'PRECIP': 'pr'},
 'meteo_file': '/tmp/tmpbnxtefbc/_data/meteo.nc'}

Préparation et calage du modèle hydrologique (xhydro.modelling)¶

INFO

Pour plus d’informations sur cette section et les options disponibles, consultez le Notebook de modélisation hydrologique.

[16]:

import xhydro as xh
from xhydro.modelling.calibration import perform_calibration
from xhydro.modelling.obj_funcs import get_objective_function

La fonction perform_calibration nécessite un argument model_config qui lui permet de construire le modèle hydrologique correspondant. Toutes les informations requises ont été acquises dans les sections précédentes, il ne reste plus qu’à remplir les entrées du modèle RavenPy/GR4JCN.

Pour simplifications, comme l’équivalent en eau de la neige n’est pas actuellement disponible dans la base de données sur PAVICS, « AVG_ANNUAL_SNOW » a été estimé à l’aide de Brown & Brasnett (2010).

[17]:

# Model configuration
model_config = {
    "model_name": "GR4JCN",
    "workdir": notebook_folder / "model",
    "overwrite": True,
    "parameters": [0.529, -3.396, 407.29, 1.072, 16.9, 0.947],
    "hru": gdf,
    "start_date": "1991-01-01",
    "end_date": "2020-12-31",
    "rain_snow_fraction": "RAINSNOW_DINGMAN",
    "evaporation": "PET_HARGREAVES_1985",
    "global_parameter": {"AVG_ANNUAL_SNOW": 100.00},
    **config_meteo_ref,  # Reuse information gathered earlier
}

# Parameter bounds for GR4JCN
bounds_low = [0.01, -15.0, 10.0, 0.0, 1.0, 0.0]
bounds_high = [2.5, 10.0, 700.0, 7.0, 30.0, 1.0]

[18]:

# Calibration / validation period
mask_calib = xr.where(qobs.time.dt.year <= 2010, 1, 0).values
mask_valid = xr.where(qobs.time.dt.year > 2010, 1, 0).values

# Model calibration
best_parameters, best_simulation, best_objfun = perform_calibration(
    model_config,
    "kge",
    qobs=notebook_folder / "_data" / "qobs.nc",
    bounds_low=bounds_low,
    bounds_high=bounds_high,
    evaluations=8,
    algorithm="DDS",
    mask=mask_calib,
    sampler_kwargs=dict(trials=1),
)

Initializing the  Dynamically Dimensioned Search (DDS) algorithm  with  8  repetitions
The objective function will be maximized
Starting the DDS algorithm with 8 repetitions...
Finding best starting point for trial 1 using 5 random samples.
1 of 8, maximal objective function=0.129327, time remaining: 00:00:10
Initialize database...
['csv', 'hdf5', 'ram', 'sql', 'custom', 'noData']
2 of 8, maximal objective function=0.129327, time remaining: 00:00:11
3 of 8, maximal objective function=0.293959, time remaining: 00:00:09
4 of 8, maximal objective function=0.293959, time remaining: 00:00:08
5 of 8, maximal objective function=0.293959, time remaining: 00:00:05
6 of 8, maximal objective function=0.293959, time remaining: 00:00:03
7 of 8, maximal objective function=0.390197, time remaining: 00:00:00
8 of 8, maximal objective function=0.393656, time remaining: 23:59:57
Best solution found has obj function value of 0.3936562927318237 at 5



*** Final SPOTPY summary ***
Total Duration: 24.88 seconds
Total Repetitions: 8
Maximal objective value: 0.393656
Corresponding parameter setting:
param0: 1.13242
param1: -2.46413
param2: 43.5229
param3: 1.33209
param4: 24.153
param5: 0.0766439
******************************

Best parameter set:
param0=1.132418270728551, param1=-2.464133107201744, param2=43.522924628974806, param3=1.3320909887768697, param4=24.153016559360484, param5=0.07664389863678882
Run number 7 has the highest objectivefunction with: 0.3937

Pour réduire les temps de calcul dans cet exemple, seulement 10 étapes ont été utilisées pour la fonction de calage, ce qui est bien en dessous du nombre recommandé. Les paramètres ci-dessous ont été obtenus en exécutant le code ci-dessus avec 150 évaluations.

[19]:

# Replace the results with parameters obtained using 150 evaluations
best_parameters = [
    0.3580270511815579,
    -2.187141388684563,
    24.012067980309702,
    0.000781,
    1.9330212374187332,
    0.5491789347783598,
]
model_config["parameters"] = best_parameters

best_simulation = xh.modelling.hydrological_model(model_config).run()

Le vrai KGE doit être calculé à partir d’une période de validation, en utilisant get_objective_function.

[20]:

get_objective_function(
    qobs=qobs.q,
    qsim=best_simulation,
    obj_func="kge",
    mask=mask_valid,
).values

[20]:

array(0.68497379)

[21]:

ax = plt.figure(figsize=(10, 5))
qobs.q.plot(color="k", linewidth=3)
best_simulation.q.plot(color="r")

[21]:

[<matplotlib.lines.Line2D at 0x7f5f930e23c0>]

Impacts des changements climatiques¶

INFO

Pour plus d’informations sur cette section et les options disponibles, consultez le Notebook d’analyse de l’impact des changements climatiques.

Cet exemple gardera l’analyse de l’impact des changements climatiques relativement simple.

Calculez la différence entre les périodes future et de référence en utilisant xhydro.cc.compute_deltas.
Utilisez ces différences pour calculer les statistiques d’ensemble en utilisant xhydro.cc.ensemble_stats : percentiles de l’ensemble et accord entre les modèles climatiques.

[38]:

# Differences
deltas = xh.cc.compute_deltas(
    ds_indicators, reference_horizon="1991-2020", kind="%", rename_variables=False
).isel(horizon=-1)

# Save the results
deltas.to_netcdf(notebook_folder / "_data" / f"deltas_sim0.nc")

deltas.squeeze()

Il existe de nombreuses façons de créer l’ensemble lui-même. Si vous utilisez un dictionnaire de jeux de données, la clé sera utilisée pour nommer chaque élément de la nouvelle dimension realization. Cela peut être très utile lorsque l’on effectue des analyses plus détaillées ou lorsque l’on souhaite pondérer les différents modèles en fonction, par exemple, du nombre de simulations disponibles. Dans notre cas, puisque nous souhaitons seulement calculer des statistiques d’ensemble, nous pouvons simplifier et fournir simplement une liste.

[39]:

import pooch

# Acquire deltas for the other 13 simulations
from xhydro.testing.helpers import (  # In-house function to access xhydro-testdata
    deveraux,
)

deltas_files = deveraux().fetch("use_case/deltas.zip", processor=pooch.Unzip())
deltas_files = xclim.ensembles.create_ensemble(deltas_files)

# Fix variable names and combine with the file we just created
deltas_files = deltas_files.rename(
    {"streamflow_max_annual": "q_max_annual", "streamflow7_min_summer": "q7_min_summer"}
)
deltas_sim0 = xr.open_dataset(
    notebook_folder / "_data" / f"deltas_sim0.nc"
).assign_coords({"realization": 13})
deltas_files = xr.concat([deltas_files, deltas_sim0], dim="realization")
clear_output(wait=False)

[40]:

# Statistics to compute
statistics = {
    "ensemble_percentiles": {"values": [10, 25, 50, 75, 90], "split": False},
    "robustness_fractions": {"test": None},
}

ens_stats = xh.cc.ensemble_stats(deltas_files, statistics)

ens_stats

[40]:

<xarray.Dataset> Size: 2kB
Dimensions:                         (percentiles: 5, month: 12, return_period: 2)
Coordinates:
  * percentiles                     (percentiles) int64 40B 10 25 50 75 90
  * month                           (month) <U3 144B 'JAN' 'FEB' ... 'NOV' 'DEC'
  * return_period                   (return_period) int64 16B 20 100
    p_quantile                      (return_period) float64 16B 0.95 0.99
    basin_name                      <U7 28B 'sub_001'
    horizon                         <U9 36B '2070-2099'
    subbasin_id                     <U1 4B '1'
    elevation                       float32 4B 222.6
    drainage_area                   float64 8B 585.6
    centroid_longitude              float64 8B -71.26
    centroid_latitude               float64 8B 46.45
Data variables: (12/32)
    qmoy_summer                     (percentiles) float64 40B -14.0 ... 19.14
    qmoy_monthly                    (month, percentiles) float64 480B dask.array<chunksize=(12, 5), meta=np.ndarray>
    q_max_annual                    (return_period, percentiles) float64 80B dask.array<chunksize=(2, 5), meta=np.ndarray>
    q7_min_summer                   (percentiles) float64 40B -28.54 ... 11.68
    qmoy_summer_changed             float64 8B 1.0
    qmoy_summer_positive            float64 8B 0.5714
    ...                              ...
    q7_min_summer_positive          float64 8B 0.2143
    q7_min_summer_changed_positive  float64 8B 0.2143
    q7_min_summer_negative          float64 8B 0.7857
    q7_min_summer_changed_negative  float64 8B 0.7857
    q7_min_summer_valid             float64 8B 1.0
    q7_min_summer_agree             float64 8B 0.7857
Attributes:
    cat:xrfreq:            fx
    cat:frequency:         fx
    cat:processing_level:  ensemble
    cat:variable:          qmoy_monthly
    ensemble_size:         4

[41]:

# Recreate the boxplots based on the computed percentiles
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(13, 5), sharey=True)

ax = plt.subplot(1, 3, 1)
for i, rp in enumerate(ens_stats.return_period.values):
    stats = [
        {
            "label": rp,
            "med": ens_stats.q_max_annual.sel(percentiles=50, return_period=rp).values,
            "q1": ens_stats.q_max_annual.sel(percentiles=25, return_period=rp).values,
            "q3": ens_stats.q_max_annual.sel(percentiles=75, return_period=rp).values,
            "whislo": ens_stats.q_max_annual.sel(
                percentiles=10, return_period=rp
            ).values,
            "whishi": ens_stats.q_max_annual.sel(
                percentiles=90, return_period=rp
            ).values,
        }
    ]

    ax.bxp(stats, showfliers=False, positions=[i], widths=0.5)
ax.set_title("Maximum annual streamflow")
plt.xlabel("Return period")
plt.ylabel("Difference Fut-Hist (%)")

ax = plt.subplot(1, 3, 2)
for i, rp in enumerate(ens_stats.return_period.values):
    stats = [
        {
            "label": rp,
            "med": ens_stats.q7_min_summer.sel(percentiles=50).values,
            "q1": ens_stats.q7_min_summer.sel(percentiles=25).values,
            "q3": ens_stats.q7_min_summer.sel(percentiles=75).values,
            "whislo": ens_stats.q7_min_summer.sel(percentiles=10).values,
            "whishi": ens_stats.q7_min_summer.sel(percentiles=90).values,
        }
    ]

    ax.bxp(stats, showfliers=False, positions=[i], widths=0.5)
ax.set_title("Minimum summer streamflow (7-day avg)")
plt.xlabel("")

ax = plt.subplot(1, 3, 3)
stats = [
    {
        "label": "",
        "med": ens_stats.qmoy_summer.sel(percentiles=50).values,
        "q1": ens_stats.qmoy_summer.sel(percentiles=25).values,
        "q3": ens_stats.qmoy_summer.sel(percentiles=75).values,
        "whislo": ens_stats.qmoy_summer.sel(percentiles=10).values,
        "whishi": ens_stats.qmoy_summer.sel(percentiles=90).values,
    }
]

ax.bxp(stats, showfliers=False, positions=[i], widths=0.25)
ax.set_title("Mean summer flow")

plt.show()

[42]:

print(
    f"Fraction of simulations with a positive change (maximum streamflow): {ens_stats.q_max_annual_positive.values}"
)
print(
    f"Fraction of simulations with a positive change (minimum summer streamflow): {ens_stats.q7_min_summer_positive.values}"
)
print(
    f"Fraction of simulations with a positive change (mean summer streamflow): {ens_stats.qmoy_summer_positive.values}"
)

Fraction of simulations with a positive change (maximum streamflow): [0.5        0.64285714]
Fraction of simulations with a positive change (minimum summer streamflow): 0.21428571428571427
Fraction of simulations with a positive change (mean summer streamflow): 0.5714285714285714

Exemple de cas d’utilisation¶

Identification d’un bassin versant et de ses caractéristiques¶

Collecte des données observées¶

Données météorologiques¶

Données hydrométriques¶

Préparation et calage du modèle hydrologique (xhydro.modelling)¶

Calcul des indicateurs hydroclimatologiques¶

Indicateurs non fréquenciels¶

Analyse fréquencielle¶

Simulations et indicateurs de débit futur¶

Données météorologiques futures¶

Données de débit et indicateurs futurs¶

Impacts des changements climatiques¶