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 l’aire de drainage, 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. Ces informations proviendront de HydroBASINS via le module xhydro.gis.

Délimitation du bassin versant et caractéristiques de base : xhydro.gis.watershed_delineation et xhydro.gis.watershed_properties
Altitude : xhydro.gis.surface_properties

[1]:

# Workaround for determining the notebook folder within a running notebook
# This cell is not visible when the documentation is built.

from __future__ import annotations

try:
    from _finder import _find_current_folder

    notebook_folder = _find_current_folder()
except ImportError:
    from pathlib import Path

    notebook_folder = Path().cwd()

(notebook_folder / "_data").mkdir(exist_ok=True)

[2]:

import leafmap
import numpy as np
import pandas as pd

import xhydro.gis as xhgis

[3]:

# Watershed delineation
coords = (-71.28878, 46.65692)
m = leafmap.Map(center=(coords[1], coords[0]), zoom=8)
gdf = xhgis.watershed_delineation(coordinates=coords, map=m)
gdf

[3]:

	HYBAS_ID	Upstream Area (sq. km).	geometry	category	color
0	7120365812	693.5	POLYGON ((-71.09758 46.40035, -71.09409 46.403...	1	#ffffd9

Afficher m produira cette carte :

[4]:

# Watershed properties
gdf_wat = xhgis.watershed_properties(gdf)
gdf_wat

[4]:

	HYBAS_ID	Upstream Area (sq. km).	category	color	area	perimeter	gravelius	centroid
0	7120365812	693.5	1	#ffffd9	5.845762e+08	158155.350574	1.845263	(-71.26046365140014, 46.45216102567744)

[5]:

# Surface properties
gdf_surf = xhgis.surface_properties(gdf)
gdf_surf

[5]:

	elevation	slope	aspect	time	id	proj:transform	platform	epsg	proj:shape	proj:epsg	title	band	gsd	spatial_ref
geometry
0	224.201508	1.416482	214.740387	2021-04-22	Copernicus_DSM_COG_30_N46_00_W072_00_DEM	{0.0, -0.0008333333333333334, 47.0004166666666...	TanDEM-X	4326	{1200}	4326	N46_00_W072_00	data	90	0

Nous conserverons les informations nécessaires et les formaterons pour GR4JCN.

[6]:

drainage_area = np.array([gdf_wat["Upstream Area (sq. km)."].iloc[0]])
latitude = np.array([gdf_wat.centroid[0][1]])
longitude = np.array([gdf_wat.centroid[0][0]])
elevation = np.array([gdf_surf["elevation"].iloc[0]])

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.

[7]:

gauge_area = 708
scaling_factor = drainage_area / gauge_area
scaling_factor

[7]:

array([0.97951977])

Collecte des données observées¶

[8]:

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.

[9]:

# 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

[9]:

<xarray.Dataset> Size: 429GB
Dimensions:  (time: 26298, lat: 800, lon: 1700)
Coordinates:
  * lat      (lat) float32 3kB 10.0 10.1 10.2 10.3 10.4 ... 89.6 89.7 89.8 89.9
  * lon      (lon) float32 7kB -179.9 -179.8 -179.7 -179.6 ... -10.2 -10.1 -10.0
  * time     (time) datetime64[ns] 210kB 1950-01-01 1950-01-02 ... 2021-12-31
Data variables:
    pr       (time, lat, lon) float32 143GB dask.array<chunksize=(365, 50, 50), meta=np.ndarray>
    tasmin   (time, lat, lon) float32 143GB dask.array<chunksize=(365, 50, 50), meta=np.ndarray>
    tasmax   (time, lat, lon) float32 143GB dask.array<chunksize=(365, 50, 50), meta=np.ndarray>
Attributes: (12/26)
    Conventions:          CF-1.8
    cell_methods:         time: mean (interval: 1 day)
    data_specs_version:   00.00.07
    domain:               NAM
    format:               netcdf
    frequency:            day
    ...                   ...
    dataset_description:  https://www.ecmwf.int/en/era5-land
    license_type:         permissive
    license:              Please acknowledge the use of ERA5-Land as stated i...
    attribution:          Contains modified Copernicus Climate Change Service...
    citation:             Muñoz Sabater, J., (2021): ERA5-Land hourly data fr...
    doi:                  https://doi.org/10.24381/cds.e2161bac

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é.

[10]:

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 = meteo_ref.assign_coords({"crs": meteo_ref.crs})
meteo_ref

[10]:

<xarray.Dataset> Size: 9MB
Dimensions:  (time: 10958, lat: 8, lon: 8)
Coordinates:
  * 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
  * time     (time) datetime64[ns] 88kB 1991-01-01 1991-01-02 ... 2020-12-31
    crs      int64 8B 1
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/27)
    Conventions:          CF-1.8
    cell_methods:         time: mean (interval: 1 day)
    data_specs_version:   00.00.07
    domain:               NAM
    format:               netcdf
    frequency:            day
    ...                   ...
    license_type:         permissive
    license:              Please acknowledge the use of ERA5-Land as stated i...
    attribution:          Contains modified Copernicus Climate Change Service...
    citation:             Muñoz Sabater, J., (2021): ERA5-Land hourly data fr...
    doi:                  https://doi.org/10.24381/cds.e2161bac
    crs:                  EPSG:4326

[11]:

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

[11]:

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

Puisque GR4JCN est un modèle global, nous devons moyenner les données météorologiques sur le bassin versant. Cela sera accompli en utilisant xscen.spatial_mean. Il existe plusieurs méthodes, mais pour les fichiers shapefile, la plus précise est « xesmf », qui utilise le SpatialAverager de xESMF pour suivre précisément la représentation de chaque point de grille sur le polygone. Cela peut être très lent pour des polygones détaillés avec de nombreux sommets, mais l’argument « simplify_tolerance » aide à accélérer ce processus.

[12]:

meteo_ref = xscen.spatial_mean(
    meteo_ref,
    method="xesmf",
    region={"method": "shape", "shape": gdf},
    simplify_tolerance=0.1,
    kwargs={"skipna": True},
)

# Some housekeeping to associate the HYBAS_ID dimension (coming from the shapefile) to other coordinates added from the GeoDataFrame.
for c in meteo_ref.coords:
    if len(meteo_ref[c].dims) == 0:
        if c not in ["lon", "lat"]:
            meteo_ref = meteo_ref.drop_vars([c])
        else:
            meteo_ref[c] = meteo_ref[c].expand_dims("HYBAS_ID")
meteo_ref

[12]:

<xarray.Dataset> Size: 351kB
Dimensions:   (time: 10958, HYBAS_ID: 1)
Coordinates:
  * time      (time) datetime64[ns] 88kB 1991-01-01 1991-01-02 ... 2020-12-31
    lon       (HYBAS_ID) float64 8B -71.28
    lat       (HYBAS_ID) float64 8B 46.46
  * HYBAS_ID  (HYBAS_ID) int64 8B 7120365812
Data variables:
    pr        (time) float64 88kB dask.array<chunksize=(355,), meta=np.ndarray>
    tasmin    (time) float64 88kB dask.array<chunksize=(355,), meta=np.ndarray>
    tasmax    (time) float64 88kB dask.array<chunksize=(355,), meta=np.ndarray>
Attributes: (12/28)
    Conventions:          CF-1.8
    cell_methods:         time: mean (interval: 1 day)
    data_specs_version:   00.00.07
    domain:               NAM
    format:               netcdf
    frequency:            day
    ...                   ...
    license:              Please acknowledge the use of ERA5-Land as stated i...
    attribution:          Contains modified Copernicus Climate Change Service...
    citation:             Muñoz Sabater, J., (2021): ERA5-Land hourly data fr...
    doi:                  https://doi.org/10.24381/cds.e2161bac
    crs:                  EPSG:4326
    regrid_method:        conservative

Raven attend des températures en Celsius et des précipitations en millimètres, mais elles sont actuellement conforme aux conventions CF, en Kelvin et en kg m-2 s-1, respectivement. Encore une fois, xscen peut être utilisé pour cette opération.

[13]:

meteo_ref = xscen.utils.change_units(
    meteo_ref, {"tasmax": "degC", "tasmin": "degC", "pr": "mm"}
)

[14]:

# Save
meteo_ref.to_netcdf(notebook_folder / "_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.

[20]:

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 simplifier, l’altitude de la station météorologique sera considérée comme l’altitude moyenne du bassin versant réel. Le calcul de l’altitude des points de grille dans ERA5-Land n’est pas toujours trivial et ne fait pas partie de l’exemple. De même, l’équivalent en eau de neige n’est actuellement pas disponible dans la base de données de PAVICS, donc « AVG_ANNUAL_SNOW » a été estimé approximativement à l’aide de Brown & Brasnett (2010).

[21]:

# Model configuration
model_config = {
    "model_name": "GR4JCN",
    "parameters": [0.529, -3.396, 407.29, 1.072, 16.9, 0.947],
    "drainage_area": drainage_area,
    "elevation": elevation,
    "latitude": latitude,
    "longitude": longitude,
    "start_date": "1991-01-01",
    "end_date": "2020-12-31",
    "meteo_file": notebook_folder / "_data" / "meteo.nc",
    "data_type": ["TEMP_MAX", "TEMP_MIN", "PRECIP"],
    "alt_names_meteo": {"TEMP_MIN": "tasmin", "TEMP_MAX": "tasmax", "PRECIP": "pr"},
    "meteo_station_properties": {
        "ALL": {"elevation": elevation, "latitude": latitude, "longitude": longitude}
    },
    "rain_snow_fraction": "RAINSNOW_DINGMAN",
    "evaporation": "PET_HARGREAVES_1985",
    "global_parameter": {"AVG_ANNUAL_SNOW": 100.00},
}

# 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]

[22]:

# 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 algotrithm with 8 repetitions...
Finding best starting point for trial 1 using 5 random samples.
1 of 8, maximal objective function=0.164749, time remaining: 00:00:30
Initialize database...
['csv', 'hdf5', 'ram', 'sql', 'custom', 'noData']
2 of 8, maximal objective function=0.164749, time remaining: 00:00:30
3 of 8, maximal objective function=0.239819, time remaining: 00:00:25
4 of 8, maximal objective function=0.466832, time remaining: 00:00:20
5 of 8, maximal objective function=0.466832, time remaining: 00:00:14
6 of 8, maximal objective function=0.466832, time remaining: 00:00:07
7 of 8, maximal objective function=0.481349, time remaining: 00:00:00
8 of 8, maximal objective function=0.481349, time remaining: 23:59:53
Best solution found has obj function value of 0.48134887704825247 at 5



*** Final SPOTPY summary ***
Total Duration: 63.2 seconds
Total Repetitions: 8
Maximal objective value: 0.481349
Corresponding parameter setting:
param0: 0.46099
param1: -9.77586
param2: 216.456
param3: 0.384076
param4: 13.2655
param5: 0.563789
******************************

Best parameter set:
param0=0.4609896022995832, param1=-9.775860334349964, param2=216.455502549914, param3=0.38407572808259516, param4=13.265462690482325, param5=0.5637885058813862
Run number 6 has the highest objectivefunction with: 0.4813

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.

[23]:

# 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.

[24]:

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

[24]:

array(0.69988518)

[25]:

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

[25]:

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

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.

[42]:

# 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.

[43]:

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.extend([notebook_folder / "_data" / f"deltas_sim0.nc"])

[44]:

# 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

[44]:

<xarray.Dataset> Size: 2kB
Dimensions:                                  (percentiles: 5, month: 12,
                                              return_period: 2)
Coordinates:
    basin_name                               <U7 28B 'sub_001'
    horizon                                  <U9 36B '2070-2099'
  * percentiles                              (percentiles) int64 40B 10 ... 90
  * month                                    (month) <U3 144B 'JAN' ... 'DEC'
  * return_period                            (return_period) int64 16B 20 100
    p_quantile                               (return_period) float64 16B 0.95...
Data variables: (12/32)
    qmoy_summer                              (percentiles) float64 40B -14.0 ...
    qmoy_monthly                             (month, percentiles) float64 480B dask.array<chunksize=(12, 5), meta=np.ndarray>
    streamflow_max_annual                    (return_period, percentiles) float64 80B dask.array<chunksize=(2, 5), meta=np.ndarray>
    streamflow7_min_summer                   (percentiles) float64 40B -28.54...
    qmoy_summer_changed                      float64 8B 1.0
    qmoy_summer_positive                     float64 8B 0.5714
    ...                                       ...
    streamflow7_min_summer_positive          float64 8B 0.2143
    streamflow7_min_summer_changed_positive  float64 8B 0.2143
    streamflow7_min_summer_negative          float64 8B 0.7857
    streamflow7_min_summer_changed_negative  float64 8B 0.7857
    streamflow7_min_summer_valid             float64 8B 1.0
    streamflow7_min_summer_agree             float64 8B 0.7857
Attributes:
    cat:xrfreq:            fx
    cat:frequency:         fx
    cat:processing_level:  ensemble
    cat:variable:          qmoy_monthly
    cat:id:
    ensemble_size:         14

[45]:

# 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.streamflow_max_annual.sel(
                percentiles=50, return_period=rp
            ).values,
            "q1": ens_stats.streamflow_max_annual.sel(
                percentiles=25, return_period=rp
            ).values,
            "q3": ens_stats.streamflow_max_annual.sel(
                percentiles=75, return_period=rp
            ).values,
            "whislo": ens_stats.streamflow_max_annual.sel(
                percentiles=10, return_period=rp
            ).values,
            "whishi": ens_stats.streamflow_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.streamflow7_min_summer.sel(percentiles=50).values,
            "q1": ens_stats.streamflow7_min_summer.sel(percentiles=25).values,
            "q3": ens_stats.streamflow7_min_summer.sel(percentiles=75).values,
            "whislo": ens_stats.streamflow7_min_summer.sel(percentiles=10).values,
            "whishi": ens_stats.streamflow7_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()

[46]:

print(
    f"Fraction of simulations with a positive change (maximum streamflow): {ens_stats.streamflow_max_annual_positive.values}"
)
print(
    f"Fraction of simulations with a positive change (minimum summer streamflow): {ens_stats.streamflow7_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¶