6. Analyses fréquentielles régionales¶
[1]:
import matplotlib.pyplot as plt
import xdatasets
from lmoments3.distr import KappaGen
from sklearn.cluster import HDBSCAN, OPTICS, AgglomerativeClustering
import xhydro as xh
import xhydro.frequency_analysis as xhfa
import xhydro.gis as xhgis
ERROR 1: PROJ: proj_create_from_database: Open of /home/docs/checkouts/readthedocs.org/user_builds/xhydro-fr/conda/latest/share/proj failed
Redefining 'percent' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
Redefining '%' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
Redefining 'year' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
Redefining 'yr' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
Redefining 'C' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
Redefining 'd' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
Redefining 'h' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
Redefining 'degrees_north' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
Redefining 'degrees_east' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
Redefining 'degrees' (<class 'pint.delegates.txt_defparser.plain.UnitDefinition'>)
Redefining '[speed]' (<class 'pint.delegates.txt_defparser.plain.DerivedDimensionDefinition'>)
Ce Notebook montrera comment utiliser la librairie xHydro
pour réaliser des analyses fréquentielles régionales sur un jeu de données de débit. Comme les étapes initiales pour l’analyse fréquentielle régionale sont les mêmes que pour l’analyse fréquentielle locale, les utilisateurs sont invités à consulter le Notebook Analyses fréquentielles locales pour un aperçu.
Dans cet exemple, nous utiliserons le même jeu de données de stations hydrométriques couvrant une partie du sud du Québec, assurant la continuité avec l’analyse précédente tout en s’étendant à une échelle régionale. Cependant, comme les analyses régionales nécessitent l’accès à plusieurs sources de données, nous allons extraire les débits pour 15 stations.
[2]:
ds = (
xdatasets.Query(
**{
"datasets": {
"deh": {
"id": ["02*"],
"regulated": ["Natural"],
"variables": ["streamflow"],
}
},
"time": {"start": "1970-01-01", "minimum_duration": (30 * 365, "d")},
}
)
.data.squeeze()
.load()
)
# This dataset lacks some attributes, so let's add them.
ds["id"].attrs["cf_role"] = "timeseries_id"
ds["streamflow"].attrs = {
"long_name": "Streamflow",
"units": "m3 s-1",
"standard_name": "water_volume_transport_in_river_channel",
"cell_methods": "time: mean",
}
# Clean some of the coordinates that are not needed for this example
ds = ds.drop_vars([c for c in ds.coords if c not in ["id", "time", "name"]])
timeargs = {
"annual": {},
}
ds_4fa = xh.indicators.get_yearly_op(
ds, op="max", timeargs=timeargs, missing="pct", missing_options={"tolerance": 0.15}
)
ds_4fa
[2]:
<xarray.Dataset> Size: 4kB Dimensions: (id: 15, time: 56) Coordinates: * id (id) object 120B '020404' '020602' ... '024013' name (id) object 120B 'York' 'Dartmouth' ... 'Bécancour' * time (time) datetime64[ns] 448B 1970-01-01 ... 2025-01-01 Data variables: streamflow_max_annual (id, time) float32 3kB nan nan nan ... nan nan nan Attributes: cat:xrfreq: YS-JAN cat:frequency: yr cat:processing_level: indicators cat:variable: ('streamflow_max_annual',) cat:id:
6.1. Variables explicatives¶
Dans les analyses fréquentielles régionales, les variables explicatives sont utilisées pour aider à expliquer la variation spatiale des extrêmes hydrologiques à travers différents lieux. Ces variables peuvent inclure des facteurs tels que la superficie du bassin versant, l’altitude, les précipitations, l’utilisation des terres et le type de sol, entre autres. En intégrant ces variables explicatives, nous pouvons prendre en compte l’influence des caractéristiques géographiques et environnementales sur les événements extrêmes, ce qui permet d’obtenir des prédictions régionales plus précises.
6.1.1. a) Extraction des caractéristiques des bassins versants à l’aide de xhydro.gis
¶
Dans cet exemple, nous utiliserons les propriétés des bassins versants comme variables explicatives principales. Cependant, d’autres variables, telles que les moyennes climatiques ou les données d’utilisation des terres, peuvent également être intégrées en fonction des exigences de l’analyse. Pour les étapes détaillées sur comment extraire et travailler avec ces données, consultez le Notebook Opérations SIG.
[3]:
gdf = xdatasets.Query(
**{
"datasets": {
"deh_polygons": {
"id": ["02*"],
"regulated": ["Natural"],
"variables": ["streamflow"],
}
},
"time": {"start": "1970-01-01", "minimum_duration": (30 * 365, "d")},
}
).data.reset_index()
dswp = xhgis.watershed_properties(
gdf[["Station", "geometry"]], unique_id="Station", output_format="xarray"
)
cent = dswp["centroid"].to_numpy()
lon = [ele[0] for ele in cent]
lat = [ele[1] for ele in cent]
dswp = dswp.assign(lon=("Station", lon))
dswp = dswp.assign(lat=("Station", lat))
dswp = dswp.drop("centroid")
dswp
[3]:
<xarray.Dataset> Size: 720B Dimensions: (Station: 15) Coordinates: * Station (Station) object 120B '020404' '020602' ... '024007' '024013' Data variables: area (Station) float64 120B 6.641e+08 6.27e+08 ... 2.331e+09 2.318e+08 perimeter (Station) float64 120B 2.496e+05 2.067e+05 ... 1.151e+05 gravelius (Station) float64 120B 2.732 2.329 2.51 ... 2.429 2.906 2.132 lon (Station) float64 120B -65.24 -64.91 -66.19 ... -71.6 -71.33 lat (Station) float64 120B 48.9 49.04 48.96 ... 46.1 46.19 46.08
6.1.2. b) Analyse en composantes principales¶
Après avoir acquis les variables explicatives, l’étape suivante consiste à les traiter en utilisant l’Analyse en Composantes Principales (ACP) pour réduire la dimensionnalité du jeu de données. L’ACP permet de simplifier le jeu de données en transformant les variables originales en un ensemble réduit de composantes non corrélées, tout en conservant la plupart des variations dans les données. Cela est réalisé en utilisant la fonction xhydro.frequency_analysis.regional.fit_pca
.
[4]:
help(xhfa.regional.fit_pca)
Help on function fit_pca in module xhydro.frequency_analysis.regional:
fit_pca(ds: xarray.core.dataset.Dataset, **kwargs) -> tuple
Perform Principal Component Analysis (PCA) on the input dataset.
This function scales the input data, applies PCA transformation, and returns
the transformed data along with the PCA object.
Parameters
----------
ds : xr.Dataset
Input dataset to perform PCA on.
\*\*kwargs : dict
Additional keyword arguments to pass to the PCA constructor.
Returns
-------
tuple: A tuple containing:
- data_pca (xr.DataArray): PCA-transformed data with 'Station' and 'components' as coordinates.
- obj_pca (sklearn.decomposition.PCA): Fitted PCA object.
Notes
-----
- The input data is scaled before PCA is applied.
- The number of components in the output depends on the n_components parameter passed to PCA.
[5]:
data_pca, pca = xhfa.regional.fit_pca(dswp, n_components=3)
data_pca
[5]:
<xarray.DataArray (Station: 15, components: 3)> Size: 360B array([[ 4.72371072e-01, 2.79138991e+00, 6.39636458e-01], [-7.78947174e-01, 2.80342902e+00, -7.03864289e-01], [ 6.42722349e-02, 2.44442638e+00, -1.70232009e-01], [-4.62904770e-01, 3.38966959e-01, 2.11537163e-01], [ 2.77366893e-01, 1.29617141e-01, -1.70186429e-03], [-2.66789043e+00, -7.45264001e-01, -3.07127663e-01], [ 1.50697107e+00, -4.68425447e-01, 2.93481446e-02], [ 1.06459845e+00, -6.33396135e-01, 1.21710218e+00], [-4.48977316e-01, -7.45172165e-01, -5.82713040e-01], [-2.67587558e-01, -9.69452822e-01, -8.49911487e-01], [-1.14083539e+00, -7.60570200e-01, 1.14297648e+00], [-5.34449889e-01, -1.00569987e+00, 8.72775507e-01], [ 3.69627958e-01, -1.10358882e+00, -4.21783113e-01], [ 4.55838155e+00, -8.44355876e-01, -5.52199238e-01], [-2.01199670e+00, -1.23190407e+00, -5.23843228e-01]]) Coordinates: * Station (Station) object 120B '020404' '020602' ... '024007' '024013' * components (components) int64 24B 0 1 2 Attributes: long_name: Fitted Scaled Data description: Fitted scaled data with StandardScaler and PCA from sk... fitted_variables: ['area', 'perimeter', 'gravelius', 'lon', 'lat']
Les résultats montrent que la corrélation entre les composantes est proche de 0, ce qui suggère que les trois premières composantes sont suffisamment indépendantes. Cela indique que ces composantes peuvent être utilisées efficacement pour le reste de notre analyse, car elles capturent la majorité de la variation des variables explicatives, sans chevauchement ou multicolinéarité significative.
[6]:
data_pca.to_dataframe(name="value").reset_index().pivot(
index="Station", columns="components"
).corr()
[6]:
value | ||||
---|---|---|---|---|
components | 0 | 1 | 2 | |
components | ||||
value | 0 | 1.000000e+00 | 1.582751e-16 | 1.692547e-16 |
1 | 1.582751e-16 | 1.000000e+00 | -4.270722e-16 | |
2 | 1.692547e-16 | -4.270722e-16 | 1.000000e+00 |
6.1.3. c) Clustering (regroupement)¶
Les résultats de l’ACP peuvent être utilisés pour regrouper les stations en clusters en fonction de leurs similarités dans les composantes principales. Le clustering aide à identifier des régions avec des caractéristiques similaires, permettant des analyses fréquentielles régionales plus ciblées et précises. Cette étape est réalisée en utilisant xhydro.frequency_analysis.regional.get_group_from_fit
, qui prend en charge les méthodes de clustering de sklearn.cluster
. Dans cet exemple, nous utiliserons AgglomerativeClustering
pour former 3 clusters en fonction des résultats de l’ACP.
[7]:
Help on function get_group_from_fit in module xhydro.frequency_analysis.regional:
get_group_from_fit(model: collections.abc.Callable, param: dict, sample: xarray.core.dataset.Dataset | xarray.core.dataarray.DataArray) -> list
Get indices of groups from a fit using the specified model and parameters.
Parameters
----------
model : callable
Model class or instance with a fit method.
param : dict
Parameters for the model.
sample : xr.Dataset or xr.DataArray
Data sample to fit the model.
Returns
-------
list :
List of indices for each non-excluded group.
[8]:
groups = xhfa.regional.get_group_from_fit(
AgglomerativeClustering, {"n_clusters": 3}, data_pca
)
groups
[8]:
[array(['022507', '022704', '023002', '023303', '023401', '023422',
'023428', '023432', '023701', '024003', '024013'], dtype=object),
array(['020404', '020602', '021407'], dtype=object),
array(['024007'], dtype=object)]
[9]:
ax = plt.subplot(1, 1, 1)
gdf[gdf["Station"].isin(groups[0])].plot(ax=ax, color="red")
gdf[gdf["Station"].isin(groups[1])].plot(ax=ax, color="green")
gdf[gdf["Station"].isin(groups[2])].plot(ax=ax, color="blue")
[9]:
<Axes: >

6.2. Analyse fréquentielle régionale¶
Hosking et Wallis ont introduit une méthode pour l’analyse fréquentielle régionale basée sur les L-moments, qui sont particulièrement utiles pour l’analyse des événements extrêmes à travers différentes régions. Voici un aperçu des principaux concepts :
L-Moments : Les L-moments sont des mesures statistiques dérivées des combinaisons linéaires des statistiques d’ordre. Ils offrent une alternative plus robuste aux moments traditionnels (comme la moyenne et la variance), en particulier en présence de valeurs aberrantes ou lorsque l’on traite de petits échantillons. Les L-moments sont largement utilisés dans les études hydrologiques pour estimer les caractéristiques des valeurs extrêmes, car ils fournissent des estimations plus stables et fiables.
Analyse fréquentielle régionale : Cette technique consiste à regrouper les données de plusieurs sites ou régions pour estimer la distribution de fréquence des événements extrêmes, tels que les inondations ou les sécheresses. L’approche permet de déterminer les paramètres statistiques communs à différents endroits, offrant ainsi une perspective plus généralisée et régionale sur le comportement des événements extrêmes. La méthode de Hosking et Wallis se concentre sur l’ajustement d’une distribution fréquentielle régionale et l’évaluation de sa capacité à capturer les caractéristiques des données.
L-Moments régionaux : Il s’agit des L-moments calculés à partir des données regroupées de plusieurs sites au sein d’une région. En appliquant les L-moments régionaux, nous pouvons estimer les paramètres des distributions fréquentielles régionales et mieux comprendre l’occurrence des événements extrêmes dans différentes parties de la région. Cette méthode améliore la précision et la robustesse des prédictions des événements extrêmes, en particulier lorsque des données de plusieurs sites sont disponibles.
Commençons par calculer les L-moments pour chaque station en utilisant xhydro.frequency_analysis.regional.calc_moments
.
[10]:
Help on function calc_moments in module xhydro.frequency_analysis.regional:
calc_moments(ds: xarray.core.dataset.Dataset) -> xarray.core.dataset.Dataset
Calculate L-moments for multiple stations.
Parameters
----------
ds : xr.Dataset
A vector of stations, where each element is an array-like object containing the data for which to calculate L-moments.
Returns
-------
xr.Dataset
L-moment dataset with a new lmom dimension.
Notes
-----
NaN values in each station are removed before calculating L-moments.
The function uses the `moment_l` function to calculate L-moments for each individual stations.
Equations are based on Hosking, J. R. M., & Wallis, J. R. (1997). Regional frequency analysis (p. 240).
[11]:
ds_moments = xhfa.regional.calc_moments(ds_4fa)
ds_moments
[11]:
<xarray.Dataset> Size: 1kB Dimensions: (id: 15, lmom: 6) Coordinates: * id (id) object 120B '020404' '020602' ... '024013' name (id) object 120B 'York' 'Dartmouth' ... 'Bécancour' * lmom (lmom) <U4 96B 'l1' 'l2' 'l3' 'tau' 'tau3' 'tau4' Data variables: streamflow_max_annual (id, lmom) float64 720B 145.1 31.03 ... 0.0697 0.1307 Attributes: cat:xrfreq: YS-JAN cat:frequency: yr cat:processing_level: indicators cat:variable: ('streamflow_max_annual',) cat:id:
Pour réaliser une analyse fréquentielle régionale, nous devons remodeler nos jeux de données de maxima annuels et nos L-moments en fonction des regroupements identifiés par l’algorithme de clustering. Comme il n’existe pas de convention de nommage standardisée pour cette nouvelle dimension, xHydro
utilise le nom group_id
pour les regroupements de clusters. Ce processus de remodelage est géré par la fonction xhydro.frequency_analysis.regional.group_ds
, qui organise les données en fonction des attributions de clusters et les prépare pour une analyse plus approfondie.
[12]:
help(xhfa.regional.group_ds)
Help on function group_ds in module xhydro.frequency_analysis.regional:
group_ds(ds: xarray.core.dataset.Dataset, *, groups: list) -> xarray.core.dataset.Dataset
Group a dataset by a list of groups.
Parameters
----------
ds : xr.Dataset
The input dataset to be grouped.
groups : list
A list of groups to be used for grouping the dataset.
Returns
-------
xr.Dataset
A new dataset with the grouped data.
[13]:
ds_groups = xhfa.regional.group_ds(ds_4fa, groups=groups)
ds_moments_groups = xhfa.regional.group_ds(ds_moments, groups=groups)
Nous devons maintenant calculer deux scores importants pour évaluer l’analyse fréquentielle régionale : le H-Score et le Z-Score. Le H-Score mesure l’homogénéité des données à travers différents sites ou régions par rapport au modèle régional. Il permet de déterminer dans quelle mesure les données de plusieurs sites s’alignent avec la distribution régionale. Le Z-Score quantifie l’écart entre les valeurs observées et attendues, standardisé par leur variabilité. Il aide à déterminer dans quelle mesure la distribution théorique (basée sur le modèle régional) correspond aux données observées.
H-Score (Score d’homogénéité)
- H < 1 : HomogèneCela indique que les données des différents sites sont assez similaires et s’alignent bien avec le modèle régional, ce qui suggère que le modèle régional est adapté pour l’ensemble de la région.
- 1 ≤ H < 2 : Peut-être homogèneCette plage suggère une certaine hétérogénéité, ce qui signifie que les données peuvent encore s’adapter généralement au modèle régional, mais qu’il existe des variations que le modèle ne capture pas entièrement.
- H ≥ 2 : HétérogèneUn score aussi élevé indique des différences significatives entre les sites ou les régions, suggérant que le modèle régional pourrait ne pas être adapté à toutes les données. La région peut être trop diverse ou nécessiter des ajustements du modèle.
Z-Score (Qualité de l’ajustement)
Z-Score faible : Un Z-Score faible suggère un bon ajustement entre le modèle et les données observées. En général, un Z-Score absolu inférieur à 1,64 indique que le modèle est approprié et que l’ajustement est statistiquement acceptable.
Z-Score élevé : Un Z-Score élevé suggère des écarts significatifs entre les valeurs observées et attendues. Un Z-Score absolu supérieur à 1,64 implique que le modèle ne correspond peut-être pas bien aux données et que des ajustements supplémentaires ou un modèle différent pourraient être nécessaires.
La fonction xhydro.frequency_analysis.regional.calc_h_z
est utilisée pour calculer le H-Score et le Z-Score pour l’analyse fréquentielle régionale.
AVERTISSEMENT
L’argument kap
dans la fonction xhydro.frequency_analysis.regional.calc_h_z
attend un objet de distribution Kappa3 généré à l’aide de lmoments3.distr.KappaGen()
. En raison de problèmes de licence, la librairie lmoments3
ne peut pas être incluse dans les exigences de xHydro
et doit être installée indépendamment.
[14]:
help(xhfa.regional.calc_h_z)
Help on function calc_h_z in module xhydro.frequency_analysis.regional:
calc_h_z(ds_groups: xarray.core.dataset.Dataset, ds_moments_groups: xarray.core.dataset.Dataset, *, kap: object, seed: int | None = None) -> xarray.core.dataset.Dataset
Calculate heterogeneity measure H and Z-score for regional frequency analysis.
Parameters
----------
ds_groups : xr.Dataset
Dataset containing grouped data.
ds_moments_groups : xr.Dataset
Dataset containing L-moments for grouped data.
kap : scipy.stats.kappa3
Kappa3 distribution object.
seed : int, optional
Random seed for reproducibility. Default is None.
Returns
-------
xr.Dataset
Dataset containing calculated H values and Z-scores for each group.
Notes
-----
This function applies the heterogeneity measure and Z-score calculations
to grouped data for regional frequency analysis. It uses L-moments and
the Kappa3 distribution in the process.
This function does not support lazy evaluation.
Equations are based on Hosking, J. R. M., & Wallis, J. R. (1997). Regional frequency analysis (p. 240).
[15]:
ds_hz = xhfa.regional.calc_h_z(ds_groups, ds_moments_groups, kap=KappaGen())
ds_hz
[15]:
<xarray.Dataset> Size: 96B Dimensions: (group_id: 3, crit: 2) Coordinates: * group_id (group_id) int64 24B 0 1 2 lmom <U4 16B 'l1' * crit (crit) <U1 8B 'H' 'Z' Data variables: streamflow_max_annual (crit, group_id) float64 48B 0.6387 0.6584 ... nan
Nous pouvons filtrer les résultats pour inclure uniquement les données ayant des scores H et Z inférieurs aux seuils donnés en utilisant la fonction xhydro.frequency_analysis.regional.mask_h_z
. Par défaut, les seuils sont définis à 1 pour H et 1,64 pour Z. Cependant, ces seuils peuvent être personnalisés en fonction des exigences spécifiques de l’analyse fréquentielle.
[16]:
help(xhfa.regional.mask_h_z)
Help on function mask_h_z in module xhydro.frequency_analysis.regional:
mask_h_z(ds: xarray.core.dataset.Dataset, *, thresh_h: float | None = 1, thresh_z: float | None = 1.64) -> xarray.core.dataarray.DataArray
Create a boolean mask based on heterogeneity measure H and Z-score thresholds.
Parameters
----------
ds : xr.Dataset
Dataset containing H and Z values for each group.
thresh_h : float, optional
Threshold for the heterogeneity measure H. Default is 1.
thresh_z : float, optional
Threshold for the absolute Z-score. Default is 1.64.
Returns
-------
xr.DataArray
Boolean mask where True indicates groups that meet both threshold criteria.
[17]:
mask = xhfa.regional.mask_h_z(ds_hz)
ds_groups_h1 = ds_groups.where(mask).load()
ds_moments_groups_h1 = ds_moments_groups.where(mask).load()
La dernière étape de l’analyse fréquentielle régionale est similaire à l’analyse fréquentielle locale, mais au lieu d’utiliser xhydro.frequency_analysis.local.parametric_quantiles
, nous utilisons xhydro.frequency_analysis.regional.calculate_rp_from_afr
pour calculer les périodes de retour. De plus, pour éviter de réaliser l’analyse sur des régions trop petites, l’argument remove_small_regions
peut être utilisé pour exclure les régions ayant moins d’un nombre spécifié de stations. Par défaut, ce seuil est fixé à 5.
[18]:
Help on function calculate_rp_from_afr in module xhydro.frequency_analysis.regional:
calculate_rp_from_afr(ds_groups: xarray.core.dataset.Dataset, ds_moments_groups: xarray.core.dataset.Dataset, *, rp: <built-in function array>, l1: xarray.core.dataarray.DataArray | None = None) -> xarray.core.dataarray.DataArray
Calculate return periods from Annual Flow Regime (AFR) analysis.
Parameters
----------
ds_groups : xr.Dataset
Dataset containing grouped flow data.
ds_moments_groups : xr.Dataset
Dataset containing L-moments for grouped data.
rp : array-like
Return periods to calculate.
l1 : xr.DataArray, optional
First L-moment (location) values. L-moment can be specified for ungauged catchments.
If None, values are taken from ds_moments_groups.
Returns
-------
xr.DataArray
Calculated return periods for each group and specified return period.
Notes
-----
This function calculates return periods using the Annual Flow Regime method.
If l1 is not provided, it uses the first L-moment from ds_moments_groups.
The function internally calls calculate_ic_from_AFR to compute the flood index.
Equations are based on Hosking, J. R. M., & Wallis, J. R. (1997). Regional frequency analysis (p. 240).
[19]:
Help on function remove_small_regions in module xhydro.frequency_analysis.regional:
remove_small_regions(ds: xarray.core.dataset.Dataset, *, thresh: int = 5) -> xarray.core.dataset.Dataset
Remove regions from the dataset that have fewer than the threshold number of stations.
Parameters
----------
ds : xr.Dataset
The input dataset containing regions and stations.
thresh : int, optional
The minimum number of stations required for a region to be kept. Default is 5.
Returns
-------
xr.Dataset
The dataset with small regions removed.
[20]:
q_t = xhfa.regional.calculate_rp_from_afr(
ds_groups_h1, ds_moments_groups_h1, rp=[2, 20, 100]
)
q_t
[20]:
<xarray.Dataset> Size: 2kB Dimensions: (group_id: 3, return_period: 3, id: 14) Coordinates: * group_id (group_id) int64 24B 0 1 2 lmom <U4 16B 'l1' * return_period (return_period) int64 24B 2 20 100 * id (id) object 112B '020404' '020602' ... '024013' name (group_id, id) object 336B nan nan nan ... nan nan Data variables: streamflow_max_annual (group_id, return_period, id) float64 1kB nan ... nan
[21]:
q_t = xhfa.regional.remove_small_regions(q_t)
q_t
[21]:
<xarray.Dataset> Size: 608B Dimensions: (group_id: 1, return_period: 3, id: 14) Coordinates: * group_id (group_id) int64 8B 0 lmom <U4 16B 'l1' * return_period (return_period) int64 24B 2 20 100 * id (id) object 112B '020404' '020602' ... '024013' name (group_id, id) object 112B nan nan ... 'Bécancour' Data variables: streamflow_max_annual (group_id, return_period, id) float64 336B nan ......
Affichons les résultats pour une station, en comparant les analyses locales et régionales.
[22]:
q_reg = q_t.sel(id="023401").dropna(dim="group_id", how="all")
reg = q_reg.streamflow_max_annual.squeeze()
[23]:
params_loc = xhfa.local.fit(ds_4fa.sel(id="023401"))
q_loc = xhfa.local.parametric_quantiles(params_loc, t=[2, 20, 100])
loc = q_loc.sel(scipy_dist="genextreme").streamflow_max_annual
[24]:
fig = plt.figure(figsize=(15, 4))
plt.plot(reg.return_period.values, reg.values, "red", label="Regional analysis")
plt.plot(loc.return_period.values, loc.values, "black", label="Local analysis")
plt.xscale("log")
plt.grid(visible=True)
plt.xlabel("Return period (years)")
plt.ylabel("Streamflow (m$^3$/s)")
plt.legend()
[24]:
<matplotlib.legend.Legend at 0x74ecbc23bad0>

6.3. Incertitudes¶
Les incertitudes sont un aspect important de l’analyse fréquentielle et doivent être prises en compte lors de l’interprétation des résultats. Ces incertitudes proviennent souvent de la qualité des données, du choix de la distribution et de l’estimation des paramètres. Bien que les visualisations puissent fournir des informations sur l’ajustement du modèle, il est essentiel de quantifier et de prendre en compte les incertitudes, telles que les intervalles de confiance pour les estimations des paramètres, afin d’assurer des conclusions robustes.
6.3.1. a) Bootstrap des observations¶
Une méthode pour quantifier les incertitudes consiste à effectuer un bootstrap des observations. Dans cet exemple, nous effectuerons le bootstrap un petit nombre de fois pour illustrer le processus, bien qu’en pratique, un nombre plus élevé d’itérations (par exemple, 5000) soit recommandé pour obtenir des estimations plus fiables. Le bootstrap consiste à rééchantillonner les données observées avec remplacement pour générer plusieurs ensembles de données, qui peuvent ensuite être utilisés pour évaluer la variabilité des paramètres et des résultats du modèle.
Cela peut être réalisé en appelant xhydro.frequency_analysis.uncertainties.bootstrap_obs
. Une différence importante avec les analyses locales est qu’aucune étape d’ajustement n’est nécessaire pour les analyses fréquentielles régionales, ce qui rend cette fonction beaucoup plus rapide. Les résultats sont obtenus en utilisant d’abord xhydro.frequency_analysis.uncertainties.calc_moments_iter
pour calculer les L-moments pour tous les échantillons de bootstrap. Ensuite, xhydro.frequency_analysis.uncertainties.calc_q_iter
est utilisé pour calculer les périodes de retour à partir de ces L-moments.
[25]:
Help on function calc_moments_iter in module xhydro.frequency_analysis.uncertainties:
calc_moments_iter(ds_samples: xarray.core.dataset.Dataset) -> xarray.core.dataset.Dataset
Calculate L-moments for each bootstrap sample.
Parameters
----------
ds_samples : xarray.Dataset
The bootstrap samples.
Returns
-------
xarray.Dataset
L-moments for each bootstrap sample.
[26]:
Help on function calc_q_iter in module xhydro.frequency_analysis.uncertainties:
calc_q_iter(bv: str, var: str, ds_groups: xarray.core.dataset.Dataset, ds_moments_iter: xarray.core.dataset.Dataset, *, return_periods: <built-in function array>, small_regions_threshold: int | None = 5, l1: xarray.core.dataarray.DataArray | None = None) -> xarray.core.dataarray.DataArray
Calculate quantiles for each bootstrap sample and group.
Parameters
----------
bv : str
The basin identifier or 'all' to proceed on all bv (needed for ungauged).
var : str
The variable name.
ds_groups : xarray.Dataset
The grouped data.
ds_moments_iter : xarray.Dataset
The L-moments for each bootstrap sample.
return_periods : array-like
The return periods to calculate quantiles for.
small_regions_threshold : int, optional
The threshold for removing small regions. Default is 5.
l1 : xr.DataArray, optional
First L-moment (location) values. L-moment can be specified for ungauged catchments.
If `None`, values are taken from ds_moments_iter.
Returns
-------
xarray.DataArray
Quantiles for each bootstrap sample and group.
[27]:
ds_reg_boot = xhfa.uncertainties.bootstrap_obs(ds_4fa, n_samples=35)
ds_moments_iter = xhfa.uncertainties.calc_moments_iter(ds_reg_boot).load()
ds_moments_iter
[27]:
<xarray.Dataset> Size: 26kB Dimensions: (samples: 35, id: 15, lmom: 6) Coordinates: * id (id) object 120B '020404' '020602' ... '024013' name (id) object 120B 'York' 'Dartmouth' ... 'Bécancour' * samples (samples) int64 280B 0 1 2 3 4 5 ... 30 31 32 33 34 * lmom (lmom) <U4 96B 'l1' 'l2' 'l3' 'tau' 'tau3' 'tau4' Data variables: streamflow_max_annual (samples, id, lmom) float64 25kB 142.9 ... 0.08028
[28]:
q_reg_boot = xhfa.uncertainties.calc_q_iter(
bv="023401",
var="streamflow_max_annual",
ds_groups=ds_groups_h1,
ds_moments_iter=ds_moments_iter,
return_periods=[2, 20, 100],
)
q_reg_boot
[28]:
<xarray.Dataset> Size: 10kB Dimensions: (return_period: 3, id: 11, samples: 35) Coordinates: lmom <U4 16B 'l1' * return_period (return_period) int64 24B 2 20 100 * id (id) object 88B '022507' '022704' ... '024013' name (id) object 88B 'Du Loup' 'Ouelle' ... 'Bécancour' * samples (samples) object 280B MultiIndex * group_id (samples) int64 280B 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 * obs_samples (samples) int64 280B 0 1 2 3 4 5 ... 30 31 32 33 34 Data variables: streamflow_max_annual (return_period, id, samples) float64 9kB 94.42 ......
6.3.2. b) Utilisation de plusieurs régions¶
Une autre approche pour estimer l’incertitude consiste à considérer plusieurs régions pour un bassin versant d’intérêt. Cela peut être réalisé en appliquant différentes méthodes de clustering ou en effectuant une procédure de jackknife sur la liste des stations. Comme ce processus peut être intensif en termes de calcul, nous présentons ici un exemple simplifié pour illustrer le potentiel de cette méthode. L’objectif principal est de démontrer comment la prise en compte de plusieurs régions peut aider à évaluer l’incertitude dans l’analyse fréquentielle régionale.
Tout d’abord, nous utiliserons xhydro.frequency_analysis.uncertainties.generate_combinations
sur les résultats de notre ACP pour générer des listes de stations. Cette fonction crée des combinaisons qui excluent un certain nombre de stations à chaque fois, nous permettant ainsi d’évaluer comment l’exclusion de différentes stations affecte les regroupements.
[29]:
Help on function generate_combinations in module xhydro.frequency_analysis.uncertainties:
generate_combinations(da: xarray.core.dataarray.DataArray, *, n: int) -> list
Generate combinations of indices omitting up to N indices.
Parameters
----------
da : xarray.DataArray
Input DataArray.
n : int
Number of indices to omit in each combination.
Returns
-------
list
List of all combinations.
[30]:
combinations_list = xhfa.uncertainties.generate_combinations(data_pca, n=2)
print(
f"This generated {len(combinations_list)} lists of stations, varying from {len(combinations_list[0])} to {len(combinations_list[-1])} stations."
)
This generated 121 lists of stations, varying from 15 to 13 stations.
Ensuite, nous combinerons ces listes avec trois méthodes de clustering et, pour chaque méthode, nous essaierons de modifier certains des paramètres.
[31]:
clust_options = {
AgglomerativeClustering: {"arg_name": "n_clusters", "range": range(2, 12)},
HDBSCAN: {"arg_name": "min_cluster_size", "range": range(6, 7)},
OPTICS: {"arg_name": "min_samples", "range": range(4, 5)},
}
[32]:
groups = []
for model in [AgglomerativeClustering, HDBSCAN, OPTICS]:
for p in clust_options[model]["range"]:
d_param = {}
d_param[clust_options[model]["arg_name"]] = p
for combination in combinations_list:
# Extract data for the current combination
data_com = data_pca.sel(Station=list(combination))
# Get groups from the fit and add to the list
groups = groups + xhfa.regional.get_group_from_fit(model, d_param, data_com)
unique_groups = [list(x) for x in {tuple(x) for x in groups}]
Les étapes suivantes suivent une approche similaire à celles précédentes, mais maintenant avec plusieurs régions.
[33]:
ds_groups_multiregion = xhfa.regional.group_ds(ds_4fa, groups=unique_groups)
ds_moments_groups_multiregion = xhfa.regional.group_ds(ds_moments, groups=unique_groups)
ds_h_z_multiregion = xhfa.regional.calc_h_z(
ds_groups_multiregion, ds_moments_groups_multiregion, kap=KappaGen()
)
mask_multiregion = xhfa.regional.mask_h_z(ds_h_z_multiregion)
ds_groups_h1_multiregion = ds_groups_multiregion.where(mask_multiregion).load()
ds_moments_groups_h1_multiregion = ds_moments_groups_multiregion.where(
mask_multiregion
).load()
q_t_multiregion = xhfa.regional.calculate_rp_from_afr(
ds_groups_h1_multiregion, ds_moments_groups_h1_multiregion, rp=[2, 20, 100]
)
q_t_multiregion = xhfa.regional.remove_small_regions(q_t_multiregion)
q = q_t_multiregion.sel(id="023401").dropna(dim="group_id", how="all")
q
[33]:
<xarray.Dataset> Size: 9kB Dimensions: (group_id: 219, return_period: 3) Coordinates: * group_id (group_id) int64 2kB 2 4 5 8 13 ... 423 424 425 428 lmom <U4 16B 'l1' * return_period (return_period) int64 24B 2 20 100 id <U6 24B '023401' name (group_id) object 2kB 'Beaurivage' ... 'Beaurivage' Data variables: streamflow_max_annual (group_id, return_period) float64 5kB 189.6 ... 418.2
6.3.3. c) Combinaison du bootstrap avec plusieurs régions¶
Il est possible de combiner les deux méthodes—bootstrap et plusieurs régions—en intégrant les éléments suivants :
Les L-moments bootstrapés de
xhydro.frequency_analysis.uncertainties.calc_moments_iter
.Les résultats masqués de
xhydro.frequency_analysis.regional.group_ds
qui ont été filtrés en fonction des scores H et Z.
Lorsque nous combinons ces méthodes, xhydro.frequency_analysis.uncertainties.calc_q_iter
vérifiera dans combien de valeurs distinctes de group_id
la station est présente et empilera les données en conséquence avec les échantillons bootstrapés. Par exemple, si nous avons 35 bootstrap et environ 220 régions, cela générera un total d’environ 7 700 échantillons pour une analyse ultérieure.
[34]:
q_reg_multiregion_boot = xhfa.uncertainties.calc_q_iter(
"023401",
"streamflow_max_annual",
ds_groups_h1_multiregion,
ds_moments_iter,
return_periods=[2, 20, 100],
)
q_reg_multiregion_boot
[34]:
<xarray.Dataset> Size: 4MB Dimensions: (return_period: 3, id: 15, samples: 7665) Coordinates: lmom <U4 16B 'l1' * return_period (return_period) int64 24B 2 20 100 * id (id) object 120B '020404' '020602' ... '024013' name (id, samples) object 920kB nan nan ... 'Bécancour' * samples (samples) object 61kB MultiIndex * group_id (samples) int64 61kB 2 2 2 2 2 ... 428 428 428 428 * obs_samples (samples) int64 61kB 0 1 2 3 4 5 ... 30 31 32 33 34 Data variables: streamflow_max_annual (return_period, id, samples) float64 3MB nan ... 1...
6.3.4. d) Comparaison¶
Illustrons la différence entre les approches.
[35]:
fig, ax = plt.subplots(1, 3)
fig.set_figheight(4)
fig.set_figwidth(15)
# Subset the data
q_reg_boots = q_reg_boot.streamflow_max_annual.sel(id="023401")
q_t_multiregions = q_t_multiregion.streamflow_max_annual.sel(id="023401")
q_t_multiregion_boots = q_reg_multiregion_boot.streamflow_max_annual.sel(id="023401")
def _make_plot(data, dim, label):
# Original fit
ax.plot(
loc.return_period.values,
loc,
"black",
label="Local frequency analysis",
)
ax.plot(
reg.return_period.values,
reg,
"red",
label="Original regional frequency analysis",
)
ax.plot(
data.return_period.values,
data.quantile(0.5, dim),
"green",
label=label,
)
data_05 = data.quantile(0.05, dim)
data_95 = data.quantile(0.95, dim)
ax.fill_between(
data.return_period.values, data_05, data_95, alpha=0.2, color="green"
)
plt.xscale("log")
plt.grid(visible=True)
plt.xlabel("Return period (years)")
plt.ylabel("Streamflow (m$^3$/s)")
plt.ylim([150, 500])
ax.legend()
ax = plt.subplot(1, 3, 1)
_make_plot(q_reg_boots, "samples", "Bootstrapped observations")
ax = plt.subplot(1, 3, 2)
_make_plot(q_t_multiregions, "group_id", "Multiple regions")
ax = plt.subplot(1, 3, 3)
_make_plot(q_t_multiregion_boots, "samples", "Combined")
