4. Frequency analysis module - Local analysis

[1]:

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

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

4.1. Extraction et préparation des données

Dans cet exemple, nous effectuerons une analyse fréquentielle en utilisant des séries chronologiques historiques provenant de différents sites. Nous commençons par obtenir un ensemble de données comprenant des informations hydrologiques. Ici, nous utilisons la librairie xdatasets pour obtenir des données hydrologiques du Ministère de l’Environnement, de la Lutte contre les changements climatiques, de la Faune et des Parcs au Québec, Canada. Plus précisément, notre requête se concentre sur les stations dont l’ID commence par 020, qui possèdent un régime d’écoulement naturel et qui se limitent à des données sur le débit des cours d’eau.

Les utilisateurs peuvent préférer générer leur propre xarray.DataArray en utilisant leur propre jeu de données. Au minimum, le xarray.DataArray utilisé pour l’analyse fréquentielle doit suivre les principes suivants :

  • L’ensemble de données a besoin d’une dimension time.

  • S’il y a une dimension spatiale, comme id dans l’exemple ci-dessous, il faut un attribut cf_role avec timeseries_id comme valeur.

  • La variable aura au moins besoin d’un attribut units, bien que d’autres attributs tels que long_name et cell_methods soient également attendus par xclim (qui est appelé à différents moments de l’analyse fréquentielle) et des avertissements seront générés s’ils sont manquants.

[2]:

ds = (
    xd.Query(
        **{
            "datasets": {
                "deh": {
                    "id": ["020*"],
                    "regulated": ["Natural"],
                    "variables": ["streamflow"],
                }
            },
            "time": {"start": "1970-01-01", "minimum_duration": (15 * 365, "d")},
        }
    )
    .data.squeeze()
    .load()
)

# This dataset lacks some of the aforementioned attributes, so we need to 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",
}

# For the purpose of this example, we keep only 2 stations
ds = ds.isel(id=slice(3, 5))
ds

[2]:

<xarray.Dataset> Size: 328kB
Dimensions:        (id: 2, time: 20454)
Coordinates: (12/15)
    drainage_area  (id) float32 8B 626.0 1.2e+03
    end_date       (id) datetime64[ns] 16B 2024-06-02 1996-08-13
  * id             (id) object 16B '020602' '020802'
    latitude       (id) float32 8B 48.98 49.2
    longitude      (id) float32 8B -64.7 -65.29
    name           (id) object 16B 'Dartmouth' 'Madeleine'
    ...             ...
    spatial_agg    <U9 36B 'watershed'
    start_date     (id) datetime64[ns] 16B 1970-10-01 1970-01-01
  * time           (time) datetime64[ns] 164kB 1970-01-01 ... 2025-12-31
    time_agg       <U4 16B 'mean'
    timestep       <U1 4B 'D'
    variable       <U10 40B 'streamflow'
Data variables:
    streamflow     (id, time) float32 164kB nan nan nan nan ... nan nan nan nan
[3]:

ds.streamflow.dropna("time", how="all").hvplot(
    x="time", grid=True, widget_location="bottom", groupby="id"
)

[3]:

4.2. Personnalisation des paramètres d’analyse

4.2.1. a) Définition des saisons

Nous pouvons définir des saisons en utilisant des indexeurs compatibles avec xclim.core.calendar.select_time. Il y a actuellement quatre types d’indexeurs acceptés :

  • month, suivi d’une séquence de numéros de mois.

  • season, suivi d’un ou plusieurs des éléments suivants : "DJF", "MAM", "JJA" et "SON".

  • doy_bounds, suivi d’une séquence représentant les limites inclusives de la période à considérer (début, fin).

  • date_bounds, qui est le même que ci-dessus, mais en utilisant le format mois-jour ('%m-%d').

Pour obtenir les maxima de blocs par xhydro.indicators.get_yearly_op, les indexeurs doivent être regroupés dans un dictionnaire, la clé étant l’étiquette à donner à la période de l’année demandée. Une deuxième clé peut être utilisée pour donner des instructions sur la fréquence de rééchantillonnage, par exemple pour envelopper l’année pour l’hiver.

[4]:

# Some examples
timeargs = {
    "spring": {"date_bounds": ["02-11", "06-19"]},
    "summer": {"doy_bounds": [152, 243]},
    "fall": {"month": [9, 10, 11]},
    "winter": {
        "season": ["DJF"],
        "freq": "YS-DEC",
    },  # To correctly wrap around the year, we need to specify the resampling frequency.
    "august": {"month": [8]},
    "annual": {},
}

4.2.2. b) Obtenir des maxima de bloc

Après avoir sélectionné chaque saison souhaitée, nous pouvons extraire les séries temporelles de maxima de blocs de chaque station en utilisant xhydro.indicators.get_yearly_op. Les principaux arguments sont :

  • op : l’opération à calculer. Parmi "max", "min", "mean", "sum".

  • input_var : le nom de la variable. La valeur par défaut est "streamflow".

  • window : la taille de la fenêtre roulante. Un "mean" est effectué sur la fenêtre mobile avant l’opération op.

  • timeargs : comme défini précédemment. Laisser à None pour obtenir les maxima annuels.

  • missing et missing_options : pour définir les tolérances pour les données manquantes. Voir cette page pour plus d’informations.

  • interpolate_na : s’il faut interpoler les données manquantes avant l’opération op. Uniquement utilisé pour sum.

La fonction renvoie un xarray.Dataset avec 1 variable par indexeur.

[5]:

# Here, we hide years with more than 15% of missing data.
ds_4fa = xh.indicators.get_yearly_op(
    ds, op="max", timeargs=timeargs, missing="pct", missing_options={"tolerance": 0.15}
)

ds_4fa

[5]:

<xarray.Dataset> Size: 4kB
Dimensions:                (id: 2, time: 56)
Coordinates: (12/15)
    drainage_area          (id) float32 8B 626.0 1.2e+03
    end_date               (id) datetime64[ns] 16B 2024-06-02 1996-08-13
  * id                     (id) object 16B '020602' '020802'
    latitude               (id) float32 8B 48.98 49.2
    longitude              (id) float32 8B -64.7 -65.29
    name                   (id) object 16B 'Dartmouth' 'Madeleine'
    ...                     ...
    spatial_agg            <U9 36B 'watershed'
    start_date             (id) datetime64[ns] 16B 1970-10-01 1970-01-01
    time_agg               <U4 16B 'mean'
    timestep               <U1 4B 'D'
    variable               <U10 40B 'streamflow'
  * time                   (time) datetime64[ns] 448B 1970-01-01 ... 2025-01-01
Data variables:
    streamflow_max_spring  (id, time) float32 448B nan 241.0 223.0 ... nan nan
    streamflow_max_summer  (id, time) float32 448B nan 27.2 123.0 ... nan nan
    streamflow_max_fall    (id, time) float32 448B nan 11.6 18.2 ... nan nan nan
    streamflow_max_august  (id, time) float32 448B nan 12.1 21.9 ... nan nan nan
    streamflow_max_annual  (id, time) float32 448B nan 241.0 223.0 ... nan nan
    streamflow_max_winter  (id, time) float32 448B 8.5 5.52 3.48 ... nan nan nan
Attributes:
    cat:frequency:         yr
    cat:processing_level:  indicators
    cat:id:                
[6]:

ds_4fa.streamflow_max_summer.dropna("time", how="all").hvplot(
    x="time", grid=True, widget_location="bottom", groupby="id"
)

[6]:

4.2.3. c) Utilisation de saisons personnalisées par année ou par station

L’utilisation de plages de dates individualisées pour chaque année ou chaque bassin versant n’est pas explicitement supportée, les utilisateurs doivent donc masquer leurs données avant d’appeler get_yearly_op. Notez que dans ce cas, missing` doit être ajusté en conséquence.

[7]:

# Create a mask beforehand

nyears = np.unique(ds.time.dt.year).size
dom_start = xr.DataArray(
    np.random.randint(1, 30, size=(nyears,)).astype("str"),
    dims=("year"),
    coords={"year": np.unique(ds.time.dt.year)},
)
dom_end = xr.DataArray(
    np.random.randint(1, 30, size=(nyears,)).astype("str"),
    dims=("year"),
    coords={"year": np.unique(ds.time.dt.year)},
)

mask = xr.zeros_like(ds["streamflow"])
for y in dom_start.year.values:
    # Random mask of dates per year, between April and June.
    mask.loc[
        {
            "time": slice(
                str(y) + "-04-" + str(dom_start.sel(year=y).item()),
                str(y) + "-06-" + str(dom_end.sel(year=y).item()),
            )
        }
    ] = 1

[8]:

mask.hvplot(x="time", grid=True, widget_location="bottom", groupby="id")

[8]:

[9]:

# The name of the indexer will be used to identify the variable created here
timeargs_custom = {"custom": {}}

# We use where() to mask the data that we want to ignore
masked = ds.where(mask == 1)
# Since we masked almost all the year, our tolerance for missing data should be changed accordingly
missing = "at_least_n"
missing_options = {"n": 45}

# We use xr.merge() to combine the results with the previous dataset.
ds_4fa = xr.merge(
    [
        ds_4fa,
        xh.indicators.get_yearly_op(
            masked,
            op="max",
            timeargs=timeargs_custom,
            missing=missing,
            missing_options=missing_options,
        ),
    ]
)

[10]:

ds_4fa.streamflow_max_custom.dropna("time", how="all").hvplot(
    x="time", grid=True, widget_location="bottom", groupby="id"
)

[10]:

4.2.4. d) Calcul des volumes

L’analyse de fréquence peut également être effectuée sur des volumes, en utilisant un flux de travail similaire. La principale différence est que si nous partons d’un débit, nous devrons d’abord le convertir en volume en utilisant xhydro.indicators.compute_volume. De plus, si nécessaire, get_yearly_op a un argument interpolate_na qui peut être utilisé pour interpoler les données manquantes avant la somme.

[11]:

# Get a daily volume from a daily streamflow
ds["volume"] = xh.indicators.compute_volume(ds["streamflow"], out_units="hm3")

# We'll take slightly different indexers
timeargs_vol = {"spring": {"date_bounds": ["04-30", "06-15"]}, "annual": {}}

# The operation that we want here is the sum, not the max.
ds_4fa = xr.merge(
    [
        ds_4fa,
        xh.indicators.get_yearly_op(
            ds,
            op="sum",
            input_var="volume",
            timeargs=timeargs_vol,
            missing="pct",
            missing_options={"tolerance": 0.15},
            interpolate_na=True,
        ),
    ]
)

[12]:

ds_4fa.volume_sum_spring.dropna("time", how="all").hvplot(
    x="time", grid=True, widget_location="bottom", groupby="id"
)

[12]:

4.3. Analyse fréquentielle locale

Une fois que nous avons nos maximums annuels (ou volumes/minimums), la première étape d’une analyse de fréquence locale est d’appeler xhfa.local.fit pour obtenir les paramètres de distribution. Les options sont les suivantes :

  • distributions : liste de distributions SciPy. La valeur par défaut est ["expon", "gamma", "genextreme", "genpareto", "gumbel_r", "pearson3", "weibull_min"].

  • min_years : nombre minimum d’années nécessaires pour ajuster les données.

  • method : méthode d’ajustement. La valeur par défaut est le maximum de vraisemblance.

[13]:

# To speed up the Notebook, we'll only perform the analysis on a subset of variables
params = xhfa.local.fit(
    ds_4fa[["streamflow_max_spring", "volume_sum_spring"]], min_years=15
)

params

[13]:

<xarray.Dataset> Size: 1kB
Dimensions:                (scipy_dist: 4, id: 2, dparams: 4)
Coordinates: (12/17)
    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'
    drainage_area          (id) float32 8B dask.array<chunksize=(2,), meta=np.ndarray>
    end_date               (id) datetime64[ns] 16B dask.array<chunksize=(2,), meta=np.ndarray>
    ...                     ...
    source                 <U102 408B 'Ministère de l’Environnement, de la Lu...
    spatial_agg            <U9 36B 'watershed'
    start_date             (id) datetime64[ns] 16B dask.array<chunksize=(2,), meta=np.ndarray>
    time_agg               <U4 16B 'mean'
    timestep               <U1 4B 'D'
    variable               <U10 40B 'streamflow'
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 d’information tels que l’AIC, le BIC et l’AICC sont utiles pour déterminer quelle distribution statistique est la mieux adaptée à un lieu donné. Ces trois critères peuvent être calculés en utilisant xhfa.local.criteria.

[14]:

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

criteria

[14]:

<xarray.Dataset> Size: 1kB
Dimensions:                (id: 2, scipy_dist: 4, criterion: 3)
Coordinates: (12/17)
    drainage_area          (id) float32 8B 626.0 1.2e+03
    end_date               (id) datetime64[ns] 16B 2024-06-02 1996-08-13
  * id                     (id) object 16B '020602' '020802'
    latitude               (id) float32 8B 48.98 49.2
    longitude              (id) float32 8B -64.7 -65.29
    name                   (id) object 16B 'Dartmouth' 'Madeleine'
    ...                     ...
    time_agg               <U4 16B 'mean'
    timestep               <U1 4B 'D'
    variable               <U10 40B 'streamflow'
    horizon                <U9 36B '1970-2025'
  * scipy_dist             (scipy_dist) <U10 160B 'genextreme' ... 'expon'
  * criterion              (criterion) <U4 48B 'aic' 'bic' 'aicc'
Data variables:
    volume_sum_spring      (scipy_dist, id, criterion) float64 192B dask.array<chunksize=(1, 2, 3), meta=np.ndarray>
    streamflow_max_spring  (scipy_dist, id, criterion) float64 192B dask.array<chunksize=(1, 2, 3), meta=np.ndarray>
Attributes:
    cat:frequency:         yr
    cat:processing_level:  indicators
    cat:id:

Enfin, les périodes de retour peuvent être obtenues en utilisant xhfa.local.parametric_quantiles. Les options sont les suivantes :

  • t : la (les) période(s) de retour en années.

  • mode : si la période de retour est la probabilité de dépassement (max) ou de non-dépassement (min). La valeur par défaut est "max".

[15]:

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

rp.load()

[15]:

<xarray.Dataset> Size: 1kB
Dimensions:                (id: 2, scipy_dist: 4, return_period: 2)
Coordinates: (12/18)
    horizon                <U9 36B '1970-2025'
  * id                     (id) object 16B '020602' '020802'
  * scipy_dist             (scipy_dist) <U10 160B 'genextreme' ... 'expon'
    drainage_area          (id) float32 8B 626.0 1.2e+03
    end_date               (id) datetime64[ns] 16B 2024-06-02 1996-08-13
    latitude               (id) float32 8B 48.98 49.2
    ...                     ...
    start_date             (id) datetime64[ns] 16B 1970-10-01 1970-01-01
    time_agg               <U4 16B 'mean'
    timestep               <U1 4B 'D'
    variable               <U10 40B 'streamflow'
  * return_period          (return_period) int64 16B 20 100
    p_quantile             (return_period) float64 16B 0.95 0.99
Data variables:
    streamflow_max_spring  (scipy_dist, return_period, id) float64 128B 303.0...
    volume_sum_spring      (scipy_dist, return_period, id) float64 128B 407.2...
Attributes:
    cat:frequency:         yr
    cat:processing_level:  indicators
    cat:id:

Dans une prochaine version, le traçage des figures sera géré par une fonction appropriée. Pour l’instant, nous allons montrer un exemple dans ce Notebook utilisant des utilitaires préliminaires.

xhfa.local._prepare_plots génère les points de données nécessaires pour tracer les résultats de l’analyse fréquentielle. Si log=True, il renverra des valeurs x espacées de logarithme entre xmin et xmax.

[16]:

data = xhfa.local._prepare_plots(params, xmin=1, xmax=1000, npoints=50, log=True)
data.load()

[16]:

<xarray.Dataset> Size: 8kB
Dimensions:                (id: 2, scipy_dist: 4, return_period: 50)
Coordinates: (12/18)
    horizon                <U9 36B '1970-2025'
  * id                     (id) object 16B '020602' '020802'
  * scipy_dist             (scipy_dist) <U10 160B 'genextreme' ... 'expon'
    drainage_area          (id) float32 8B 626.0 1.2e+03
    end_date               (id) datetime64[ns] 16B 2024-06-02 1996-08-13
    latitude               (id) float32 8B 48.98 49.2
    ...                     ...
    start_date             (id) datetime64[ns] 16B 1970-10-01 1970-01-01
    time_agg               <U4 16B 'mean'
    timestep               <U1 4B 'D'
    variable               <U10 40B 'streamflow'
  * return_period          (return_period) float64 400B 1.0 1.151 ... 1e+03
    p_quantile             (return_period) float64 400B 0.0 0.1315 ... 0.999
Data variables:
    streamflow_max_spring  (scipy_dist, return_period, id) float64 3kB -inf ....
    volume_sum_spring      (scipy_dist, return_period, id) float64 3kB -inf ....
Attributes:
    cat:frequency:         yr
    cat:processing_level:  indicators
    cat:id: