9. Extreme value analysis¶

This module provides a wrapper for the Extremes.jl package, enabling its integration with Xarray. The Extremes.jl package allows the analysis of extreme values, providing the following features:

Methods for block maxima and threshold exceedances (e.g., genextreme, gumbel_r, genpareto).
Parameter estimation using probability-weighted moments, maximum likelihood, and Bayesian procedures.
Stationary and non-stationary models.
Return level estimation.

For more information about the Extremes.jl package, see: https://doi.org/10.18637/jss.v109.i06 and https://github.com/jojal5/Extremes.jl

[1]:

import holoviews as hv
import hvplot.xarray  # noqa
import numpy as np
import pandas as pd
import panel as pn
import pooch

from xhydro.extreme_value_analysis.parameterestimation import fit, return_level
from xhydro.testing.helpers import deveraux

Detected IPython. Loading juliacall extension. See https://juliapy.github.io/PythonCall.jl/stable/compat/#IPython

9.1. Data¶

This example uses GFDL ESM4 annual total precipitation data from 1955 to 2100 for 97 virtual stations in Quebec.

[2]:

eva_data = deveraux().fetch("extreme_value_analysis/NOAA_GFDL_ESM4.zip", pooch.Unzip())

df = pd.read_csv(eva_data[0], parse_dates=[0])[
    ["time", "station_num", "station_name", "total_precip"]
]
ds = df.to_xarray()
ds = ds.set_coords(["time", "station_num", "station_name"])
ds = ds.set_index(index=["station_num", "time"])
ds = ds.unstack("index")

Downloading file 'extreme_value_analysis/NOAA_GFDL_ESM4.zip' from 'https://raw.githubusercontent.com/hydrologie/xhydro-testdata/v2025.1.14/data/extreme_value_analysis/NOAA_GFDL_ESM4.zip' to 'C:\Users\ospin\AppData\Local\xhydro-testdata\xhydro-testdata\Cache\v2025.1.14'.
Unzipping contents of 'C:\Users\ospin\AppData\Local\xhydro-testdata\xhydro-testdata\Cache\v2025.1.14\extreme_value_analysis\NOAA_GFDL_ESM4.zip' to 'C:\Users\ospin\AppData\Local\xhydro-testdata\xhydro-testdata\Cache\v2025.1.14\extreme_value_analysis\NOAA_GFDL_ESM4.zip.unzip'

[3]:

station_sel = pn.widgets.Select(
    name="Select Station", options=ds.coords["station_num"].values.tolist()
)


def plot_station_data(station_num):
    station_data = ds.total_precip.sel(station_num=station_num)
    plot = station_data.hvplot.line(x="time", title=f"Virtual station {station_num}")
    return plot


interactive_plot = pn.bind(plot_station_data, station_num=station_sel)
pn.Column(station_sel, interactive_plot).servable()

[3]:

9.6. Comparison of the return level using the stationary and non-stationary model¶

[9]:

station_selector = pn.widgets.Select(
    name="Select Station",
    options=rtnlv_non_stationary.coords["station_num"].values.tolist(),
)


def plot_station_data(station_num):

    non_stationary = rtnlv_non_stationary.sel(station_num=station_num)
    stationary = rtnlv_stationary.sel(station_num=station_num)

    line_plot_sta = hv.Curve(
        (non_stationary.return_level.values, stationary.total_precip.values[0]),
        label="Stationary",
    ).opts(color="blue")
    area_plot = hv.Area(
        (
            list(non_stationary.return_level.values),
            stationary.total_precip_lower.values[0],
            stationary.total_precip_upper.values[0],
        ),
        vdims=["y", "y2"],
    ).opts(alpha=0.3, color="blue", line_width=0)

    area_plot_non = non_stationary.hvplot.area(
        x="return_level",
        y="total_precip_lower",
        y2="total_precip_upper",
        alpha=0.3,
        color="red",
        line_width=0,
    )
    line_plot_non = non_stationary["total_precip"].hvplot.line(
        x="return_level", label="Non-stationary", color="red"
    )

    plot = area_plot * line_plot_sta * area_plot_non * line_plot_non

    plot = plot.opts(
        title=f"Virtual Station {station_num}",
        xlabel="Return Level (Time)",
        ylabel="Total Precipitation",
    )

    return plot


interactive_plot = pn.bind(plot_station_data, station_num=station_selector)


pn.Column(station_selector, interactive_plot).servable()

[9]:

9.7. Working with dask.array chunks¶

Currently, the extreme_value_analysis function is not thread-safe. Therefore, it is recommended to use the dask.scheduler= »processes » option to ensures that tasks are executed in separate Python processes.

[10]:

ds_c = ds.chunk({"time": -1, "station_num": 1})

fit_stationary_c = fit(
    ds_c,
    dist="genextreme",
    method="ml",
    variables=["total_precip"],
    confidence_level=0.95,
)
fit_stationary_c

C:\Users\ospin\OneDrive\Documents\Github\xhydro\src\xhydro\extreme_value_analysis\parameterestimation.py:215: UserWarning: Dataset contains chunks. It is recommended to use scheduler='processes' to compute the results.

[10]:

<xarray.Dataset> Size: 8kB
Dimensions:             (station_num: 97, dparams: 3)
Coordinates:
  * station_num         (station_num) int64 776B 1001 1004 1008 ... 15403 15414
  * dparams             (dparams) <U5 60B 'shape' 'loc' 'scale'
Data variables:
    total_precip        (station_num, dparams) float64 2kB dask.array<chunksize=(1, 3), meta=np.ndarray>
    total_precip_lower  (station_num, dparams) float64 2kB dask.array<chunksize=(1, 3), meta=np.ndarray>
    total_precip_upper  (station_num, dparams) float64 2kB dask.array<chunksize=(1, 3), meta=np.ndarray>
Attributes:
    long_name:         genextreme parameters
    dist:              genextreme
    method:            Maximum likelihood
    confidence_level:  0.95

[11]:

fit_stationary_c.compute(scheduler="processes")

Detected IPython. Loading juliacall extension. See https://juliapy.github.io/PythonCall.jl/stable/compat/#IPython
Detected IPython. Loading juliacall extension. See https://juliapy.github.io/PythonCall.jl/stable/compat/#IPython
Detected IPython. Loading juliacall extension. See https://juliapy.github.io/PythonCall.jl/stable/compat/#IPython
Detected IPython. Loading juliacall extension. See https://juliapy.github.io/PythonCall.jl/stable/compat/#IPython
Detected IPython. Loading juliacall extension. See https://juliapy.github.io/PythonCall.jl/stable/compat/#IPython
Detected IPython. Loading juliacall extension. See https://juliapy.github.io/PythonCall.jl/stable/compat/#IPython
Detected IPython. Loading juliacall extension. See https://juliapy.github.io/PythonCall.jl/stable/compat/#IPython
Detected IPython. Loading juliacall extension. See https://juliapy.github.io/PythonCall.jl/stable/compat/#IPython
Detected IPython. Loading juliacall extension. See https://juliapy.github.io/PythonCall.jl/stable/compat/#IPython
Detected IPython. Loading juliacall extension. See https://juliapy.github.io/PythonCall.jl/stable/compat/#IPython
Detected IPython. Loading juliacall extension. See https://juliapy.github.io/PythonCall.jl/stable/compat/#IPython
Detected IPython. Loading juliacall extension. See https://juliapy.github.io/PythonCall.jl/stable/compat/#IPython
Detected IPython. Loading juliacall extension. See https://juliapy.github.io/PythonCall.jl/stable/compat/#IPython
Detected IPython. Loading juliacall extension. See https://juliapy.github.io/PythonCall.jl/stable/compat/#IPython
Detected IPython. Loading juliacall extension. See https://juliapy.github.io/PythonCall.jl/stable/compat/#IPython
Detected IPython. Loading juliacall extension. See https://juliapy.github.io/PythonCall.jl/stable/compat/#IPython
Detected IPython. Loading juliacall extension. See https://juliapy.github.io/PythonCall.jl/stable/compat/#IPython

9. Extreme value analysis¶

9.1. Data¶

9.2. Parameter estimation for stationary model¶

9.3. Return level estimation for stationary model¶

9.4. Parameter estimation for non-stationary model¶

9.5. Return level estimation for non-stationary model¶

9.6. Comparison of the return level using the stationary and non-stationary model¶

9.7. Working with dask.array chunks¶