5. Analyses fréquentielles locales¶

[1]:

# Basic imports
import hvplot.xarray  # noqa
import numpy as np
import xarray as xr
import xdatasets

import xhydro as xh
import xhydro.frequency_analysis as xhfa

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'>)

5.3. Analyse fréquentielle locale¶

Après avoir extrait les données brutes, telles que les maximums ou minimums annuels, l’analyse fréquentielle locale est réalisée en trois étapes :

Utiliser xhydro.frequency_analysis.local.fit pour déterminer le meilleur ensemble de paramètres pour un nombre donné de distributions statistiques.
(Facultatif) Utiliser xhydro.frequency_analysis.local.criteria pour calculer les critères de qualité d’ajustement et évaluer à quel point chaque distribution statistique correspond bien aux données.
Utiliser xhydro.frequency_analysis.local.parametric_quantiles pour calculer les niveaux de retour basés sur les paramètres ajustés.

Pour accélérer le Notebook, nous effectuerons l’analyse uniquement sur un sous-ensemble de variables.

[15]:

help(xhfa.local.fit)

Help on function fit in module xhydro.frequency_analysis.local:

fit(ds, *, distributions: list[str] | None = None, min_years: int | None = None, method: str = 'ML', periods: list[str] | list[list[str]] | None = None) -> xarray.core.dataset.Dataset
    Fit multiple distributions to data.

    Parameters
    ----------
    ds : xr.Dataset
        Dataset containing the data to fit. All variables will be fitted.
    distributions : list of str, optional
        List of distribution names as defined in `scipy.stats`. See https://docs.scipy.org/doc/scipy/reference/stats.html#continuous-distributions.
        Defaults to ["genextreme", "pearson3", "gumbel_r", "expon"].
    min_years : int, optional
        Minimum number of years required for a distribution to be fitted.
    method : str
        Fitting method. Defaults to "ML" (maximum likelihood).
    periods : list of str or list of list of str, optional
        Either [start, end] or list of [start, end] of periods to be considered.
        If multiple periods are given, the output will have a `horizon` dimension.
        If None, all data is used.

    Returns
    -------
    xr.Dataset
        Dataset containing the parameters of the fitted distributions, with a new dimension `scipy_dist` containing the distribution names.

    Notes
    -----
    In order to combine the parameters of multiple distributions, the size of the `dparams` dimension is set to the
    maximum number of unique parameters between the distributions.

La fonction fit permet d’ajuster simultanément plusieurs distributions statistiques, telles que ["genextreme", "pearson3", "gumbel_r", "expon"]. Comme différentes distributions ont des ensembles de paramètres variables (et parfois des conventions de nommage différentes), xHydro gère cette complexité en utilisant une dimension dparams et en remplissant avec des valeurs NaN là où nécessaire. Lorsque les résultats interagissent avec SciPy, comme la fonction parametric_quantiles, seuls les paramètres pertinents pour chaque distribution sont passés. Les distributions sélectionnées sont stockées dans une nouvelle dimension scipy_dist.

[16]:

params = xhfa.local.fit(
    ds_4fa[["streamflow_max_spring", "volume_sum_spring"]], min_years=15
)

params

[16]:

<xarray.Dataset> Size: 820B
Dimensions:                (scipy_dist: 4, id: 2, dparams: 4)
Coordinates:
    horizon                <U9 36B '1970-2025'
  * id                     (id) object 16B '020602' '020802'
  * dparams                (dparams) <U5 80B 'c' 'skew' 'loc' 'scale'
  * scipy_dist             (scipy_dist) <U10 160B 'genextreme' ... 'expon'
    name                   (id) object 16B dask.array<chunksize=(2,), meta=np.ndarray>
Data variables:
    streamflow_max_spring  (scipy_dist, id, dparams) float64 256B dask.array<chunksize=(1, 2, 4), meta=np.ndarray>
    volume_sum_spring      (scipy_dist, id, dparams) float64 256B dask.array<chunksize=(1, 2, 4), meta=np.ndarray>
Attributes:
    cat:frequency:         yr
    cat:processing_level:  indicators
    cat:id:

Les critères tels que l’AIC (Critère d’Information d’Akaike), le BIC (Critère Bayésien d’Information), et l’AICC (AIC corrigé) sont des outils précieux pour comparer l’ajustement de différents modèles statistiques. Ces critères équilibrent la qualité d’ajustement d’un modèle avec sa complexité, aidant à éviter l’overfitting. L’AIC et l’AICC sont particulièrement utiles pour comparer des modèles avec des nombres de paramètres différents, tandis que le BIC a tendance à pénaliser plus sévèrement la complexité, ce qui le rend plus conservateur. Des valeurs plus faibles de ces critères indiquent de meilleures performances du modèle, l’AICC étant particulièrement utile pour les petits échantillons. En utilisant ces métriques, nous pouvons déterminer objectivement le modèle le plus approprié pour nos données.

Ces trois critères peuvent être consultés en utilisant xhydro.frequency_analysis.local.criteria. Les résultats sont ajoutés à une nouvelle dimension criterion. Dans cet exemple, l’AIC, le BIC et l’AICC donnent tous une faible indication que la distribution Generalized Extreme Value (GEV) est la meilleure correspondance pour les données, bien que la distribution de Gumbel puisse également être un choix approprié. En revanche, le Pearson III n’a pas convergé et la distribution exponentielle a été rejetée en fonction de ces critères, ce qui suggère qu’elles ne correspondent pas correctement aux données.

[17]:

help(xhfa.local.criteria)

Help on function criteria in module xhydro.frequency_analysis.local:

criteria(ds: xarray.core.dataset.Dataset, p: xarray.core.dataset.Dataset) -> xarray.core.dataset.Dataset
    Compute information criteria (AIC, BIC, AICC) from fitted distributions, using the log-likelihood.

    Parameters
    ----------
    ds : xr.Dataset
        Dataset containing the yearly data that was fitted.
    p : xr.Dataset
        Dataset containing the parameters of the fitted distributions.
        Must have a dimension `dparams` containing the parameter names and a dimension `scipy_dist` containing the distribution names.

    Returns
    -------
    xr.Dataset
        Dataset containing the information criteria for the distributions.

[18]:

criteria = xhfa.local.criteria(
    ds_4fa[["streamflow_max_spring", "volume_sum_spring"]], params
)

criteria["streamflow_max_spring"].isel(id=0).compute()

Enfin, les périodes de retour peuvent être obtenues en utilisant xhfa.local.parametric_quantiles.

[19]:

help(xhfa.local.parametric_quantiles)

Help on function parametric_quantiles in module xhydro.frequency_analysis.local:

parametric_quantiles(p: xarray.core.dataset.Dataset, t: float | list[float], mode: str = 'max') -> xarray.core.dataset.Dataset
    Compute quantiles from fitted distributions.

    Parameters
    ----------
    p : xr.Dataset
        Dataset containing the parameters of the fitted distributions.
        Must have a dimension `dparams` containing the parameter names and a dimension `scipy_dist` containing the distribution names.
    t : float or list of float
        Return period(s) in years.
    mode : {'max', 'min'}
        Whether the return period is the probability of exceedance (max) or non-exceedance (min).

    Returns
    -------
    xr.Dataset
        Dataset containing the quantiles of the distributions.

[20]:

rp = xhfa.local.parametric_quantiles(params, t=[2, 20, 100])

rp.load()

Dans une future version, l’affichage des graphiques sera géré par une fonction dédiée. Pour l’instant, nous démontrons le processus en utilisant les utilitaires préliminaires dans ce Notebook.

La fonction xhfa.local._prepare_plots génère les points de données nécessaires pour visualiser les résultats de l’analyse fréquentielle. Si log=True, elle renverra des valeurs x espacées logarithmiquement entre xmin et xmax. En parallèle, xhfa.local._get_plotting_positions calcule les positions des graphiques pour toutes les variables dans le jeu de données. Elle accepte les arguments alpha et beta. Pour des valeurs typiques, consultez la documentation SciPy. Par défaut, (0.4, 0.4) est utilisé, ce qui correspond à la méthode des quantiles non biaisée (Cunnane).

[21]:

data = xhfa.local._prepare_plots(params, xmin=1, xmax=1000, npoints=50, log=True)
pp = xhfa.local._get_plotting_positions(ds_4fa[["streamflow_max_spring"]])

[22]:

# Let's plot the observations
p1 = data.streamflow_max_spring.hvplot(
    x="return_period", by="scipy_dist", grid=True, groupby=["id"], logx=True
)
data.streamflow_max_spring.hvplot(
    x="return_period", by="scipy_dist", grid=True, groupby=["id"], logx=True
)

# Let's now plot the distributions
p2 = pp.hvplot.scatter(
    x="streamflow_max_spring_pp",
    y="streamflow_max_spring",
    grid=True,
    groupby=["id"],
    logx=True,
)
pp.hvplot.scatter(
    x="streamflow_max_spring_pp",
    y="streamflow_max_spring",
    grid=True,
    groupby=["id"],
    logx=True,
)

# And now combining the plots
p1 * p2

[22]:

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

Pour gérer l’intensité computationnelle, nous nous concentrerons sur un seul bassin versant et limiterons l’analyse aux deux distributions qui semblaient les mieux correspondre aux données.

[23]:

ds_4fa = ds_4fa.sel(id="020602")[["streamflow_max_spring"]]
params = params.sel(id="020602", scipy_dist=["genextreme", "gumbel_r"])[
    ["streamflow_max_spring"]
]

5.4.2. b) Bootstrap des distributions¶

Dans cette approche, plutôt que de rééchantillonner directement les observations, nous rééchantillonnons les distributions ajustées pour estimer l’incertitude. Cette méthode nous permet d’évaluer la variabilité des paramètres des distributions ajustées. Comme dans l’exemple précédent, nous effectuerons un nombre minimal d’itérations de bootstrap pour réduire la charge computationnelle, mais en pratique, un nombre plus élevé d’itérations fournirait des estimations plus robustes de l’incertitude.

Cela peut être accompli en appelant xhydro.frequency_analysis.uncertainties.bootstrap_dist. Contrairement à bootstrap_obs, cette méthode ne prend pas en charge l’évaluation paresseuse et requiert une fonction spécifique pour l’étape d’ajustement : xhydro.frequency_analysis.uncertainties.fit_boot_dist.

[27]:

help(xhfa.uncertainties.bootstrap_dist)

Help on function bootstrap_dist in module xhydro.frequency_analysis.uncertainties:

bootstrap_dist(ds_obs: xarray.core.dataset.Dataset, ds_params: xarray.core.dataset.Dataset, *, n_samples: int) -> xarray.core.dataset.Dataset
    Generate bootstrap samples from a fitted distribution.

    Parameters
    ----------
    ds_obs : xarray.Dataset
        The observed data.
    ds_params : xarray.Dataset
        The fitted distribution parameters.
    n_samples : int
        The number of bootstrap samples to generate.

    Returns
    -------
    xarray.Dataset
        Bootstrap samples with dimensions [samples, time].

    Notes
    -----
    This function does not support lazy evaluation.

[28]:

help(xhfa.uncertainties.fit_boot_dist)

Help on function fit_boot_dist in module xhydro.frequency_analysis.uncertainties:

fit_boot_dist(ds: xarray.core.dataset.Dataset) -> xarray.core.dataset.Dataset
    Fit distributions to bootstrap samples.

    Parameters
    ----------
    ds : xarray.Dataset
        The bootstrap samples to fit.

    Returns
    -------
    xarray.Dataset
        Fitted distribution parameters for each bootstrap sample.

[29]:

tmp_values = xhfa.uncertainties.bootstrap_dist(ds_4fa, params, n_samples=35)
params_boot_dist = xhfa.uncertainties.fit_boot_dist(tmp_values)
params_boot_dist

[29]:

<xarray.Dataset> Size: 2kB
Dimensions:                (scipy_dist: 2, samples: 35, dparams: 3)
Coordinates:
  * samples                (samples) int64 280B 0 1 2 3 4 5 ... 30 31 32 33 34
  * dparams                (dparams) <U5 60B 'c' 'loc' 'scale'
  * scipy_dist             (scipy_dist) <U10 80B 'genextreme' 'gumbel_r'
    horizon                <U9 36B '1970-2025'
    id                     <U6 24B '020602'
    name                   object 8B 'Dartmouth'
Data variables:
    streamflow_max_spring  (scipy_dist, samples, dparams) float64 2kB dask.array<chunksize=(1, 35, 2), meta=np.ndarray>

[30]:

rp_boot_dist = xhfa.local.parametric_quantiles(
    params_boot_dist.load(), t=[2, 20, 100]
)  # Lazy computing is not supported
rp_boot_dist

5.4.3. c) Comparaison¶

Illustrons la différence entre les deux approches.

[31]:

import matplotlib.pyplot as plt

fig, ax = plt.subplots()
fig.set_figheight(4)
fig.set_figwidth(15)

# Subset the data
rp_orig = rp.streamflow_max_spring.sel(id="020602", scipy_dist="genextreme")
boot_obs = rp_boot_obs.streamflow_max_spring.sel(scipy_dist="genextreme")
boot_dist = rp_boot_dist.streamflow_max_spring.sel(scipy_dist="genextreme")

# Original fit
ax.plot(
    rp_orig.return_period.values,
    rp_orig,
    "black",
    label="Original fit",
)

ax.plot(
    boot_obs.return_period.values,
    boot_obs.quantile(0.5, "samples"),
    "red",
    label="Bootstrapped observations",
)
boot_obs_05 = boot_obs.quantile(0.05, "samples")
boot_obs_95 = boot_obs.quantile(0.95, "samples")
ax.fill_between(
    boot_obs.return_period.values, boot_obs_05, boot_obs_95, alpha=0.2, color="red"
)

ax.plot(
    boot_dist.return_period.values,
    boot_dist.quantile(0.5, "samples"),
    "blue",
    label="Bootstrapped distribution",
)
boot_dist_05 = boot_dist.quantile(0.05, "samples")
boot_dist_95 = boot_dist.quantile(0.95, "samples")
ax.fill_between(
    boot_dist.return_period.values, boot_dist_05, boot_dist_95, alpha=0.2, color="blue"
)

plt.xscale("log")
plt.grid(visible=True)
plt.xlabel("Return period (years)")
plt.ylabel("Streamflow (m$^3$/s)")
ax.legend()

[31]:

<matplotlib.legend.Legend at 0x72905c698200>

../_images/notebooks_local_frequency_analysis_48_1.png

5. Analyses fréquentielles locales¶

5.1. Extraction et préparation des données¶

5.2. Acquisition des maxima par blocs¶

5.2.1. a) Définir les saisons¶

5.2.2. b) Définir les critères pour les valeurs manquantes¶

5.2.3. c) Exemple simple¶

5.2.4. d) Exemple avancé : Utiliser des saisons personnalisées par année ou par station¶

5.2.5. e) Variable alternative : Calcul des volumes¶

5.3. Analyse fréquentielle locale¶

5.4. Incertitudes¶

5.4.1. a) Bootstrap des observations¶

5.4.2. b) Bootstrap des distributions¶

5.4.3. c) Comparaison¶