Skip to content

Extremes

Block-extrema (block maxima / minima for a generalized extreme-value workflow) and peaks-over-threshold (POT) diagnostics: threshold selection, exceedances, and Poisson-process counts and statistics.

Functions

These pure functions back the operators above; each takes xr.DataArray/xr.Dataset and a dim: argument.

block_maxima 

block_maxima(da: DataArray, block_size: int = 365, time: str = 'time', boundary: str = 'trim', side: str = 'left') -> DataArray

Coarsen along time and take the maximum per block.

Parameters:

Name Type Description Default
da DataArray

Input DataArray with a time dimension.

 required
block_size int

Number of time steps per block.

 365
time str

Name of the time dimension.

 'time'
boundary str

Passed to .coarsen ("trim", "pad").

 'trim'
side str

Passed to .coarsen ("left" or "right").

 'left'

Returns:

Type Description
DataArray

DataArray of block maxima.
Source code in packages/xrtoolz-ops/src/xrops/extremes/_src/extremes.py 
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
def block_maxima(
    da: xr.DataArray,
    block_size: int = 365,
    time: str = "time",
    boundary: str = "trim",
    side: str = "left",
) -> xr.DataArray:
    """Coarsen along ``time`` and take the maximum per block.

    Args:
        da: Input DataArray with a ``time`` dimension.
        block_size: Number of time steps per block.
        time: Name of the time dimension.
        boundary: Passed to ``.coarsen`` (``"trim"``, ``"pad"``).
        side: Passed to ``.coarsen`` (``"left"`` or ``"right"``).

    Returns:
        DataArray of block maxima.
    """
    return da.coarsen({time: block_size}, boundary=boundary, side=side).max()

block_minima 

block_minima(da: DataArray, block_size: int = 365, time: str = 'time', boundary: str = 'trim', side: str = 'left') -> DataArray

Coarsen along time and take the minimum per block.

Source code in packages/xrtoolz-ops/src/xrops/extremes/_src/extremes.py 
39
40
41
42
43
44
45
46
47
def block_minima(
    da: xr.DataArray,
    block_size: int = 365,
    time: str = "time",
    boundary: str = "trim",
    side: str = "left",
) -> xr.DataArray:
    """Coarsen along ``time`` and take the minimum per block."""
    return da.coarsen({time: block_size}, boundary=boundary, side=side).min()

pot_threshold 

pot_threshold(da: DataArray, quantile: float = 0.98, time: str = 'time') -> DataArray

Threshold for peaks-over-threshold: the quantile of da.

Source code in packages/xrtoolz-ops/src/xrops/extremes/_src/extremes.py 
50
51
52
53
54
55
56
def pot_threshold(
    da: xr.DataArray,
    quantile: float = 0.98,
    time: str = "time",
) -> xr.DataArray:
    """Threshold for peaks-over-threshold: the ``quantile`` of ``da``."""
    return da.quantile(quantile, dim=time)

pot_exceedances 

pot_exceedances(da: DataArray, quantile: float = 0.98, decluster_freq: int | None = None, time: str = 'time') -> DataArray

Return values above the quantile threshold, optionally declustered.

When decluster_freq is set, values are declustered by taking the block maximum over that many steps; this converts the peaks series into independent exceedances (standard POT practice).

Parameters:

Name Type Description Default
da DataArray

Input DataArray.

 required
quantile float

Quantile (in [0, 1]) used to compute the threshold.

 0.98
decluster_freq int | None

If given, apply a block-maximum declustering filter with this block size and drop resulting NaNs.

 None
time str

Name of the time dimension.

 'time'

Returns:

Type Description
DataArray

DataArray of exceedances. If decluster_freq is None the
DataArray

returned array has the same shape as da with sub-threshold
DataArray

points replaced by NaN.
Source code in packages/xrtoolz-ops/src/xrops/extremes/_src/extremes.py 
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
def pot_exceedances(
    da: xr.DataArray,
    quantile: float = 0.98,
    decluster_freq: int | None = None,
    time: str = "time",
) -> xr.DataArray:
    """Return values above the ``quantile`` threshold, optionally declustered.

    When ``decluster_freq`` is set, values are declustered by taking the
    block maximum over that many steps; this converts the peaks series
    into independent exceedances (standard POT practice).

    Args:
        da: Input DataArray.
        quantile: Quantile (in [0, 1]) used to compute the threshold.
        decluster_freq: If given, apply a block-maximum declustering
            filter with this block size and drop resulting NaNs.
        time: Name of the time dimension.

    Returns:
        DataArray of exceedances. If ``decluster_freq is None`` the
        returned array has the same shape as ``da`` with sub-threshold
        points replaced by NaN.
    """
    threshold = pot_threshold(da, quantile=quantile, time=time)
    above = da.where(da >= threshold)
    if decluster_freq is None:
        return above
    return block_maxima(above, block_size=decluster_freq, time=time).dropna(time)

pp_counts 

pp_counts(da: DataArray, quantile: float = 0.98, block_size: int = 5, time: str = 'time', boundary: str = 'trim', side: str = 'left') -> DataArray

Count exceedances of a high quantile within each time block.

Parameters:

Name Type Description Default
da DataArray

Input DataArray.

 required
quantile float

Quantile used to build the threshold.

 0.98
block_size int

Size of the time block in which to count.

 5
time str

Name of the time dimension.

 'time'
boundary, side

Passed to :func:xarray.DataArray.coarsen.

 required

Returns:

Type Description
DataArray

DataArray with one count per block along the (coarsened) time
DataArray

dimension.
Source code in packages/xrtoolz-ops/src/xrops/extremes/_src/extremes.py 
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
def pp_counts(
    da: xr.DataArray,
    quantile: float = 0.98,
    block_size: int = 5,
    time: str = "time",
    boundary: str = "trim",
    side: str = "left",
) -> xr.DataArray:
    """Count exceedances of a high quantile within each time block.

    Args:
        da: Input DataArray.
        quantile: Quantile used to build the threshold.
        block_size: Size of the time block in which to count.
        time: Name of the time dimension.
        boundary, side: Passed to :func:`xarray.DataArray.coarsen`.

    Returns:
        DataArray with one count per block along the (coarsened) time
        dimension.
    """
    return _pp_apply(
        da,
        quantile=quantile,
        block_size=block_size,
        time=time,
        boundary=boundary,
        side=side,
        reducer=lambda vals, thr: int(np.sum(vals > thr)),
    )

pp_stats 

pp_stats(da: DataArray, quantile: float = 0.98, block_size: int = 5, statistic: Callable[[ndarray], float] = mean, time: str = 'time', boundary: str = 'trim', side: str = 'left') -> DataArray

Summarize exceedances within each time block.

Applies statistic to the subset of values in each block that exceed the quantile threshold.

Parameters:

Name Type Description Default
da DataArray

Input DataArray.

 required
quantile float

Quantile used to build the threshold.

 0.98
block_size int

Size of the time block.

 5
statistic Callable[[ndarray], float]

Reduction to apply to the exceedances (default numpy.mean).

 mean
time str

Name of the time dimension.

 'time'
boundary, side

Passed to :func:xarray.DataArray.coarsen.

 required

Returns:

Type Description
DataArray

DataArray of per-block summary statistics.
Source code in packages/xrtoolz-ops/src/xrops/extremes/_src/extremes.py 
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
def pp_stats(
    da: xr.DataArray,
    quantile: float = 0.98,
    block_size: int = 5,
    statistic: Callable[[np.ndarray], float] = np.mean,
    time: str = "time",
    boundary: str = "trim",
    side: str = "left",
) -> xr.DataArray:
    """Summarize exceedances within each time block.

    Applies ``statistic`` to the subset of values in each block that
    exceed the ``quantile`` threshold.

    Args:
        da: Input DataArray.
        quantile: Quantile used to build the threshold.
        block_size: Size of the time block.
        statistic: Reduction to apply to the exceedances (default
            ``numpy.mean``).
        time: Name of the time dimension.
        boundary, side: Passed to :func:`xarray.DataArray.coarsen`.

    Returns:
        DataArray of per-block summary statistics.
    """
    return _pp_apply(
        da,
        quantile=quantile,
        block_size=block_size,
        time=time,
        boundary=boundary,
        side=side,
        reducer=lambda vals, thr: (
            float(statistic(vals[vals > thr])) if np.any(vals > thr) else np.nan
        ),
    )