Extraction

xarray-native extremes extraction — the preprocessing that turns a raw spatiotemporal series into the sample an EVT model actually expects. Inputs are xr.DataArrays; array-valued results preserve coordinates, dimensions, and metadata, while fully reduced results (constant_threshold, quantile_threshold reduced over all dimensions, estimate_extremal_index on 1-D data) come back as plain floats.

The two routes through this module mirror the two EVT paradigms: block maxima feed a GEV, and threshold exceedances feed a GPD.

Block maxima

Take the maximum over fixed blocks (usually years) and model those with a GEV. Simple and robust, but it discards all but one observation per block — r_largest_block_maxima recovers some of that lost information by keeping the top \(r\) order statistics, and sliding_block_maxima trades independence for a larger sample.

Extraction utilities for extreme values on xarray data.

Block maxima, threshold selection, and declustering — the three standard preprocessing steps for extreme value analysis on spatiotemporal data.

temporal_block_maxima(da: xr.DataArray, freq: str, time_dim: str = 'time', min_periods: int | None = None, keep_time: bool = False) -> xr.DataArray

Extract temporal block maximum values from an xarray DataArray.

Uses time-based resampling to compute maximum values over temporal blocks.

Parameters:

Name	Type	Description	Default
da	DataArray	Input data array	required
freq	str	Frequency string for temporal grouping (e.g., 'YS' for annual, 'MS' for monthly, '7D' for weekly, '3h' for 3-hourly)	required
time_dim	str	Name of the time dimension	'time'
min_periods	int	Minimum number of observations required for a valid maximum. If None, all blocks are used regardless of missing data.	None
keep_time	bool	If True, attach a f"{time_dim}_of_max" coordinate holding the original timestamp at which each block maximum occurred (rather than only the resampled block label). Empty or min_periods-masked blocks get NaT. Useful for plotting the maxima at the moment they actually happened instead of at the block boundary.	False

Returns:

Type	Description
DataArray	Temporal block maximum values. When keep_time=True the result also carries a f"{time_dim}_of_max" coordinate with the argmax timestamps.

Examples:

>>> # Annual maxima
>>> annual_max = temporal_block_maxima(da, 'YS')

>>> # Monthly maxima with at least 20 observations per month
>>> monthly_max = temporal_block_maxima(da, 'MS', min_periods=20)

>>> # Seasonal maxima (3-month blocks)
>>> seasonal_max = temporal_block_maxima(da, '3MS')

>>> # Annual maxima labelled with the day each maximum occurred
>>> am = temporal_block_maxima(da, 'YS', keep_time=True)
>>> am['time_of_max']  # actual timestamps, NaT for empty blocks

Source code in src/xtremax/extraction/block_maxima.py

def temporal_block_maxima(
    da: xr.DataArray,
    freq: str,
    time_dim: str = "time",
    min_periods: int | None = None,
    keep_time: bool = False,
) -> xr.DataArray:
    """
    Extract temporal block maximum values from an xarray DataArray.

    Uses time-based resampling to compute maximum values over temporal blocks.

    Parameters
    ----------
    da : xr.DataArray
        Input data array
    freq : str
        Frequency string for temporal grouping (e.g., 'YS' for annual,
        'MS' for monthly, '7D' for weekly, '3h' for 3-hourly)
    time_dim : str, default 'time'
        Name of the time dimension
    min_periods : int, optional
        Minimum number of observations required for a valid maximum.
        If None, all blocks are used regardless of missing data.
    keep_time : bool, default False
        If True, attach a ``f"{time_dim}_of_max"`` coordinate holding the
        **original timestamp** at which each block maximum occurred (rather than
        only the resampled block label). Empty or ``min_periods``-masked blocks
        get ``NaT``. Useful for plotting the maxima at the moment they actually
        happened instead of at the block boundary.

    Returns
    -------
    xr.DataArray
        Temporal block maximum values. When ``keep_time=True`` the result also
        carries a ``f"{time_dim}_of_max"`` coordinate with the argmax timestamps.

    Examples
    --------
    >>> # Annual maxima
    >>> annual_max = temporal_block_maxima(da, 'YS')

    >>> # Monthly maxima with at least 20 observations per month
    >>> monthly_max = temporal_block_maxima(da, 'MS', min_periods=20)

    >>> # Seasonal maxima (3-month blocks)
    >>> seasonal_max = temporal_block_maxima(da, '3MS')

    >>> # Annual maxima labelled with the day each maximum occurred
    >>> am = temporal_block_maxima(da, 'YS', keep_time=True)
    >>> am['time_of_max']  # actual timestamps, NaT for empty blocks
    """
    resampler = da.resample({time_dim: freq})

    maxima = resampler.max(dim=time_dim, keep_attrs=True)
    if min_periods is not None:
        # Mask blocks with insufficient data.
        counts = resampler.count(dim=time_dim)
        maxima = maxima.where(counts >= min_periods)

    if keep_time:
        na = _na_for(np.asarray(da[time_dim].values).dtype)

        def _arg_time(block: xr.DataArray) -> xr.DataArray:
            stamps = xr.DataArray(np.asarray(block[time_dim].values), dims=[time_dim])
            return _argmax_coord(block, stamps, [time_dim], na)

        arg_time = da.resample({time_dim: freq}).map(_arg_time)
        # Align NaT with the maxima's own NaNs (e.g. min_periods masking).
        arg_time = arg_time.where(maxima.notnull())
        maxima = maxima.assign_coords({f"{time_dim}_of_max": arg_time})

    return maxima

spatial_block_maxima(da: xr.DataArray, block_size: int | dict[Hashable, int], dims: Hashable | Sequence[Hashable] = None, min_periods: int | None = None, boundary: str = 'trim', keep_coords: bool = False) -> xr.DataArray

Extract spatial block maximum values from an xarray DataArray.

Uses coarsening to compute maximum values over spatial blocks of fixed size.

Parameters:

Name	Type	Description	Default
da	DataArray	Input data array	required
block_size	int or dict	Size of spatial blocks. If int, uses same size for all specified dims. If dict, maps dimension names to block sizes (e.g., {'x': 10, 'y': 10})	required
dims	Hashable or Sequence[Hashable]	Dimension(s) to coarsen. If None and block_size is int, raises error. If None and block_size is dict, uses dict keys.	None
min_periods	int	Minimum number of valid observations required within each block. If None, computes maximum even with missing data.	None
boundary	str	How to handle boundaries. Options: 'trim', 'pad', 'exact'	'trim'
keep_coords	bool	If True, attach a f"{dim}_of_max" coordinate for each coarsened dim holding the original coordinate at which each block maximum occurred (i.e. where in the block the extreme is), rather than only the coarsened block label. Empty / masked blocks get NaN.	False

Returns:

Type	Description
DataArray	Spatial block maximum values with reduced resolution. When keep_coords=True the result also carries a f"{dim}_of_max" coordinate per coarsened dim with the argmax locations.

Examples:

>>> # 10x10 spatial blocks
>>> spatial_max = spatial_block_maxima(da, {'x': 10, 'y': 10})

>>> # 5x5 blocks along lat/lon dimensions
>>> spatial_max = spatial_block_maxima(da, 5, dims=['lat', 'lon'])

>>> # Require at least 50 valid points per block
>>> spatial_max = spatial_block_maxima(
...     da, {'x': 10, 'y': 10}, min_periods=50
... )

>>> # 5x5 blocks, keeping the (lat, lon) of each block's maximum
>>> sm = spatial_block_maxima(da, 5, dims=['lat', 'lon'], keep_coords=True)
>>> sm['lat_of_max'], sm['lon_of_max']

Source code in src/xtremax/extraction/block_maxima.py

def spatial_block_maxima(
    da: xr.DataArray,
    block_size: int | dict[Hashable, int],
    dims: Hashable | Sequence[Hashable] = None,
    min_periods: int | None = None,
    boundary: str = "trim",
    keep_coords: bool = False,
) -> xr.DataArray:
    """
    Extract spatial block maximum values from an xarray DataArray.

    Uses coarsening to compute maximum values over spatial blocks of fixed size.

    Parameters
    ----------
    da : xr.DataArray
        Input data array
    block_size : int or dict
        Size of spatial blocks. If int, uses same size for all specified dims.
        If dict, maps dimension names to block sizes (e.g., {'x': 10, 'y': 10})
    dims : Hashable or Sequence[Hashable], optional
        Dimension(s) to coarsen. If None and block_size is int,
        raises error. If None and block_size is dict, uses dict keys.
    min_periods : int, optional
        Minimum number of valid observations required within each block.
        If None, computes maximum even with missing data.
    boundary : str, default 'trim'
        How to handle boundaries. Options: 'trim', 'pad', 'exact'
    keep_coords : bool, default False
        If True, attach a ``f"{dim}_of_max"`` coordinate for each coarsened
        ``dim`` holding the **original coordinate** at which each block maximum
        occurred (i.e. *where* in the block the extreme is), rather than only the
        coarsened block label. Empty / masked blocks get NaN.

    Returns
    -------
    xr.DataArray
        Spatial block maximum values with reduced resolution. When
        ``keep_coords=True`` the result also carries a ``f"{dim}_of_max"``
        coordinate per coarsened ``dim`` with the argmax locations.

    Examples
    --------
    >>> # 10x10 spatial blocks
    >>> spatial_max = spatial_block_maxima(da, {'x': 10, 'y': 10})

    >>> # 5x5 blocks along lat/lon dimensions
    >>> spatial_max = spatial_block_maxima(da, 5, dims=['lat', 'lon'])

    >>> # Require at least 50 valid points per block
    >>> spatial_max = spatial_block_maxima(
    ...     da, {'x': 10, 'y': 10}, min_periods=50
    ... )

    >>> # 5x5 blocks, keeping the (lat, lon) of each block's maximum
    >>> sm = spatial_block_maxima(da, 5, dims=['lat', 'lon'], keep_coords=True)
    >>> sm['lat_of_max'], sm['lon_of_max']
    """
    # Handle block_size specification
    if isinstance(block_size, dict):
        coarsen_kwargs = block_size
    else:
        if dims is None:
            raise ValueError("Must specify 'dims' when block_size is an integer")
        if isinstance(dims, str):
            dims = [dims]
        coarsen_kwargs = {dim: block_size for dim in dims}

    # Create coarsen object
    coarsened = da.coarsen(dim=coarsen_kwargs, boundary=boundary)

    maxima = coarsened.max(keep_attrs=True)
    if min_periods is not None:
        counts = coarsened.count()
        maxima = maxima.where(counts >= min_periods)

    if keep_coords:
        cdims = list(coarsen_kwargs)
        taken = set(da.dims) | set(da.coords)
        win: dict[Hashable, str] = {}
        for d in cdims:
            win[d] = _unique_dim(f"_{d}_win", taken)
            taken.add(win[d])
        # Split each coarsened dim into (block, within-block) and find the
        # argmax over the flattened within-block window.
        con = coarsened.construct({d: (d, win[d]) for d in cdims})
        window_dims = [win[d] for d in cdims]
        # Align the argmax locations with the maxima's own NaNs, so empty and
        # ``min_periods``-masked blocks report no location.
        valid = maxima.notnull()
        # Build each coordinate's within-block grid, broadcast over every window
        # dim so the argmax runs over all of them at once — robust to any
        # number of coarsened dims.
        coord_windows = [
            xr.DataArray(da[d].values, dims=[d])
            .coarsen({d: coarsen_kwargs[d]}, boundary=boundary)
            .construct({d: (d, win[d])})
            for d in cdims
        ]
        for d, cwin in zip(cdims, xr.broadcast(*coord_windows), strict=True):
            d_of_max = _argmax_coord(con, cwin, window_dims, _na_for(cwin.dtype)).where(
                valid
            )
            maxima = maxima.assign_coords({f"{d}_of_max": d_of_max})

    return maxima

sliding_block_maxima(da: xr.DataArray, window_size: int, dim: str = 'time', stride: int = 1, min_periods: int | None = None, center: bool = False, keep_time: bool = False) -> xr.DataArray

Extract block maxima using a sliding window approach.

Sliding windows provide overlapping blocks, which helps: - Avoid sensitivity to arbitrary block boundaries - Increase the number of block maxima for better statistical estimation - Capture extreme events that might span fixed block boundaries

Parameters:

Name	Type	Description	Default
da	DataArray	Input data array	required
window_size	int	Size of the sliding window (number of observations)	required
dim	str	Dimension along which to slide the window	'time'
stride	int	Step size between windows. stride=1 gives fully overlapping windows, stride=window_size gives non-overlapping blocks (equivalent to coarsen)	1
min_periods	int	Minimum number of valid observations required for a valid maximum. If None, defaults to window_size.	None
center	bool	If True, set the window labels at the center of the window	False
keep_time	bool	If True, attach a f"{dim}_of_max" coordinate holding the original coordinate at which each window's maximum occurred (which can differ from the window's own label). Windows whose maximum is masked — no valid data, or fewer than min_periods observations — get NaN/NaT.	False

Returns:

Type	Description
DataArray	Sliding block maximum values. When keep_time=True the result also carries a f"{dim}_of_max" coordinate with the argmax coordinates.

Examples:

>>> # 365-day sliding window with daily stride
>>> sliding_max = sliding_block_maxima(da, window_size=365, stride=1)

>>> # Non-overlapping 30-day blocks (equivalent to fixed blocks)
>>> block_max = sliding_block_maxima(da, window_size=30, stride=30)

>>> # 90-day window with 7-day stride (weekly updates)
>>> weekly_max = sliding_block_maxima(da, window_size=90, stride=7)

>>> # Track the day each window maximum actually occurred
>>> sm = sliding_block_maxima(da, window_size=30, keep_time=True)
>>> sm['time_of_max']

Source code in src/xtremax/extraction/block_maxima.py

def sliding_block_maxima(
    da: xr.DataArray,
    window_size: int,
    dim: str = "time",
    stride: int = 1,
    min_periods: int | None = None,
    center: bool = False,
    keep_time: bool = False,
) -> xr.DataArray:
    """
    Extract block maxima using a sliding window approach.

    Sliding windows provide overlapping blocks, which helps:
    - Avoid sensitivity to arbitrary block boundaries
    - Increase the number of block maxima for better statistical estimation
    - Capture extreme events that might span fixed block boundaries

    Parameters
    ----------
    da : xr.DataArray
        Input data array
    window_size : int
        Size of the sliding window (number of observations)
    dim : str, default 'time'
        Dimension along which to slide the window
    stride : int, default 1
        Step size between windows. stride=1 gives fully overlapping windows,
        stride=window_size gives non-overlapping blocks (equivalent to coarsen)
    min_periods : int, optional
        Minimum number of valid observations required for a valid maximum.
        If None, defaults to window_size.
    center : bool, default False
        If True, set the window labels at the center of the window
    keep_time : bool, default False
        If True, attach a ``f"{dim}_of_max"`` coordinate holding the **original
        coordinate** at which each window's maximum occurred (which can differ
        from the window's own label). Windows whose maximum is masked — no valid
        data, or fewer than ``min_periods`` observations — get NaN/NaT.

    Returns
    -------
    xr.DataArray
        Sliding block maximum values. When ``keep_time=True`` the result also
        carries a ``f"{dim}_of_max"`` coordinate with the argmax coordinates.

    Examples
    --------
    >>> # 365-day sliding window with daily stride
    >>> sliding_max = sliding_block_maxima(da, window_size=365, stride=1)

    >>> # Non-overlapping 30-day blocks (equivalent to fixed blocks)
    >>> block_max = sliding_block_maxima(da, window_size=30, stride=30)

    >>> # 90-day window with 7-day stride (weekly updates)
    >>> weekly_max = sliding_block_maxima(da, window_size=90, stride=7)

    >>> # Track the day each window maximum actually occurred
    >>> sm = sliding_block_maxima(da, window_size=30, keep_time=True)
    >>> sm['time_of_max']
    """
    if min_periods is None:
        min_periods = window_size

    # Use rolling window
    rolling = da.rolling({dim: window_size}, center=center, min_periods=min_periods)
    maxima = rolling.max()

    if keep_time:
        wdim = _unique_dim("_w", set(da.dims) | set(da.coords))
        windows = rolling.construct(wdim)
        coord_windows = (
            xr.DataArray(da[dim].values, dims=[dim], coords={dim: da[dim]})
            .rolling({dim: window_size}, center=center)
            .construct(wdim)
        )
        # Align with the maxima's own NaNs so windows with no data — or too few
        # valid observations for ``min_periods`` — report no coordinate.
        arg_coord = _argmax_coord(
            windows, coord_windows, [wdim], _na_for(coord_windows.dtype)
        ).where(maxima.notnull())
        # Assigned before striding so the coordinate is subsampled with it.
        maxima = maxima.assign_coords({f"{dim}_of_max": arg_coord})

    # Apply stride by subsampling, starting from the first *valid*
    # window — the first position whose window holds at least
    # min_periods samples. With the default min_periods == window_size
    # that is the first complete window, so stride == window_size
    # reproduces non-overlapping blocks (coarsen equivalence); an
    # explicit smaller min_periods keeps the valid partial edge windows
    # instead of discarding them. Rolling labels sit at the window end
    # (center=False) or centre (center=True, window spanning
    # [i - w//2, i + (w-1)//2]).
    if stride > 1:
        if center:
            offset = max(0, min_periods - 1 - (window_size - 1) // 2)
        else:
            offset = min_periods - 1
        maxima = maxima.isel({dim: slice(offset, None, stride)})

    return maxima

r_largest_block_maxima(da: xr.DataArray, block_size: int | str, r: int = 3, dim: str = 'time', min_periods: int | None = None, keep_time: bool = False) -> xr.DataArray

Extract r-largest order statistics from each block.

Instead of just the maximum, extract the r largest values from each block. This provides more information for extreme value analysis and can improve statistical efficiency.

Parameters:

Name	Type	Description	Default
da	DataArray	Input data array	required
block_size	int or str	Size of blocks. If int, number of elements. If str, frequency string.	required
r	int	Number of largest values to extract from each block	3
dim	str	Dimension along which to compute order statistics	'time'
min_periods	int	Minimum number of observations required in a block	None
keep_time	bool	If True, attach a f"{dim}_of_max" coordinate (with the same order axis) holding the original coordinate at which each of the r largest values occurred. Missing slots get NaN/NaT.	False

Returns:

Type	Description
DataArray	Array with r-largest values from each block, with new dimension 'order'. For integer block_size, a trailing partial block (fewer than block_size samples) is trimmed and does not contribute values. When keep_time=True it also carries a f"{dim}_of_max" coordinate over (..., order) with the timestamps of those r values.

Examples:

>>> # 3 largest values from each year
>>> annual_r_largest = r_largest_block_maxima(da, 'YS', r=3)

>>> # Top 5 values from 100-element blocks
>>> block_r_largest = r_largest_block_maxima(da, 100, r=5)

>>> # ... and the day each of the 3 largest occurred
>>> rl = r_largest_block_maxima(da, 'YS', r=3, keep_time=True)
>>> rl['time_of_max']

Notes

The r-largest order statistics model is useful when you want to use more information than just the block maximum for fitting extreme value distributions.

Source code in src/xtremax/extraction/block_maxima.py

def r_largest_block_maxima(
    da: xr.DataArray,
    block_size: int | str,
    r: int = 3,
    dim: str = "time",
    min_periods: int | None = None,
    keep_time: bool = False,
) -> xr.DataArray:
    """
    Extract r-largest order statistics from each block.

    Instead of just the maximum, extract the r largest values from each block.
    This provides more information for extreme value analysis and can improve
    statistical efficiency.

    Parameters
    ----------
    da : xr.DataArray
        Input data array
    block_size : int or str
        Size of blocks. If int, number of elements. If str, frequency string.
    r : int, default 3
        Number of largest values to extract from each block
    dim : str, default 'time'
        Dimension along which to compute order statistics
    min_periods : int, optional
        Minimum number of observations required in a block
    keep_time : bool, default False
        If True, attach a ``f"{dim}_of_max"`` coordinate (with the same
        ``order`` axis) holding the **original coordinate** at which each of the
        r largest values occurred. Missing slots get NaN/NaT.

    Returns
    -------
    xr.DataArray
        Array with r-largest values from each block, with new dimension 'order'.
        For integer ``block_size``, a trailing partial block (fewer than
        ``block_size`` samples) is trimmed and does not contribute values.
        When ``keep_time=True`` it also carries a ``f"{dim}_of_max"`` coordinate
        over ``(..., order)`` with the timestamps of those r values.

    Examples
    --------
    >>> # 3 largest values from each year
    >>> annual_r_largest = r_largest_block_maxima(da, 'YS', r=3)

    >>> # Top 5 values from 100-element blocks
    >>> block_r_largest = r_largest_block_maxima(da, 100, r=5)

    >>> # ... and the day each of the 3 largest occurred
    >>> rl = r_largest_block_maxima(da, 'YS', r=3, keep_time=True)
    >>> rl['time_of_max']

    Notes
    -----
    The r-largest order statistics model is useful when you want to use
    more information than just the block maximum for fitting extreme
    value distributions.
    """

    def _top_r_1d(values: np.ndarray) -> np.ndarray:
        """Top-r descending from a 1-D array, NaN-padded to length r.

        Operating per 1-D slice is what keeps non-``dim`` axes independent
        — flattening a multi-dim block pools every non-``dim`` axis into
        one sample and corrupts per-site r-largest extraction.
        """
        clean = values[~np.isnan(values)]
        out = np.full(r, np.nan, dtype=float)
        if min_periods is not None and clean.size < min_periods:
            return out
        if clean.size == 0:
            return out
        sorted_asc = np.sort(clean)
        k = min(clean.size, r)
        out[:k] = sorted_asc[-k:][::-1]
        return out

    def _top_r_pos_1d(values: np.ndarray) -> np.ndarray:
        """Original indices of the r largest (descending), -1-padded to r."""
        valid = np.flatnonzero(~np.isnan(values))
        out = np.full(r, -1, dtype=np.int64)
        if min_periods is not None and valid.size < min_periods:
            return out
        if valid.size == 0:
            return out
        order = valid[np.argsort(values[valid])][::-1]
        k = min(order.size, r)
        out[:k] = order[:k]
        return out

    def _apply(values: xr.DataArray, func, out_dtype) -> xr.DataArray:
        return xr.apply_ufunc(
            func,
            values,
            input_core_dims=[[dim]],
            output_core_dims=[["order"]],
            vectorize=True,
            output_dtypes=[out_dtype],
            dask="parallelized",
            dask_gufunc_kwargs={
                "output_sizes": {"order": r},
                # Core-dim inputs chunked along `dim` must be
                # consolidated before the gufunc runs.
                "allow_rechunk": True,
            },
        )

    def _coord_from_pos(pos: xr.DataArray, coords: np.ndarray) -> xr.DataArray:
        # .where fills unselected slots with NaN (NaT for datetime coords).
        picked = coords[np.clip(pos.values, 0, None)]
        out = xr.DataArray(picked, dims=pos.dims, coords=pos.coords)
        return out.where(pos >= 0)

    if isinstance(block_size, str):
        # Time-based resampling; delegate to ``_apply`` over `dim` so
        # non-``dim`` axes are processed independently per group.
        groups = da.resample({dim: block_size})
        result = groups.map(lambda g: _apply(g, _top_r_1d, float))
        result = result.assign_coords(order=np.arange(1, r + 1))
        result.attrs = da.attrs
        if keep_time:
            # The block index of each top-r value, mapped to original coords
            # per block (positions are within-block, so resolve per group).
            def _group_times(g: xr.DataArray) -> xr.DataArray:
                pos = _apply(g, _top_r_pos_1d, np.int64)
                return _coord_from_pos(pos, np.asarray(g[dim].values))

            times = groups.map(_group_times)
            result = result.assign_coords({f"{dim}_of_max": times})
        return result
    else:
        # Fixed-size blocks. Label each position along `dim` with its
        # block id and groupby; for each block, run the 1-D selector per
        # non-``dim`` slice via ``_apply``. This preserves the multi-dim
        # structure instead of flattening every non-``dim`` axis into
        # the block-wise top-r pool.
        n_blocks = da.sizes[dim] // block_size
        if n_blocks == 0:
            raise ValueError(
                f"Series has {da.sizes[dim]} samples along {dim!r}, shorter "
                f"than one block of size {block_size}; no block maxima can "
                "be extracted."
            )
        trimmed_length = n_blocks * block_size
        trimmed = da.isel({dim: slice(0, trimmed_length)})
        block_ids = np.arange(trimmed_length) // block_size
        trimmed = trimmed.assign_coords(_block=(dim, block_ids))

        result = trimmed.groupby("_block").map(lambda g: _apply(g, _top_r_1d, float))
        block_coords = da[dim].values[:trimmed_length:block_size]
        result = result.rename({"_block": "block"}).assign_coords(
            block=block_coords, order=np.arange(1, r + 1)
        )
        result.attrs = da.attrs
        if keep_time:

            def _block_times(g: xr.DataArray) -> xr.DataArray:
                pos = _apply(g, _top_r_pos_1d, np.int64)
                return _coord_from_pos(pos, np.asarray(g[dim].values))

            times = trimmed.groupby("_block").map(_block_times)
            times = times.rename({"_block": "block"}).assign_coords(
                block=block_coords, order=np.arange(1, r + 1)
            )
            result = result.assign_coords({f"{dim}_of_max": times})
        return result

declustered_block_maxima(da: xr.DataArray, threshold: float, min_separation: int, dim: str = 'time', method: str = 'runs', keep_time: bool = False) -> xr.DataArray

Extract declustered block maxima by identifying independent extreme events.

Declustering removes temporal dependence in extremes by ensuring sufficient separation between events. This is crucial for extreme value theory which assumes independence.

Parameters:

Name	Type	Description	Default
da	DataArray	Input data array	required
threshold	float	Threshold value for identifying extreme events	required
min_separation	int	Minimum separation (in time steps) between independent events	required
dim	str	Dimension along which to decluster	'time'
method	str	Declustering method: - 'runs': Select maximum from each exceedance run, merging runs whose below-threshold gap is shorter than min_separation - 'separation': Enforce minimum time separation between peaks	'runs'
keep_time	bool	If True, attach a f"{dim}_of_max" coordinate holding the original timestamp at which each cluster's maximum occurred. The 'runs' method otherwise labels the representative at the end of its run, so this recovers the true peak time; for 'separation' the peaks already sit at their original positions.	False

Returns:

Type	Description
DataArray	Declustered block maxima with reduced temporal dependence.

Examples:

>>> # Extract independent storm maxima (3-day separation)
>>> storm_max = declustered_block_maxima(
...     da, threshold=100, min_separation=3, method='runs'
... )

>>> # Heat wave maxima with 7-day separation
>>> heatwave_max = declustered_block_maxima(
...     da, threshold=35, min_separation=7, method='separation'
... )

Notes

The runs method identifies continuous periods above threshold and takes the maximum from each run. The separation method ensures peaks are separated by at least min_separation time steps.

Source code in src/xtremax/extraction/block_maxima.py

def declustered_block_maxima(
    da: xr.DataArray,
    threshold: float,
    min_separation: int,
    dim: str = "time",
    method: str = "runs",
    keep_time: bool = False,
) -> xr.DataArray:
    """
    Extract declustered block maxima by identifying independent extreme events.

    Declustering removes temporal dependence in extremes by ensuring sufficient
    separation between events. This is crucial for extreme value theory which
    assumes independence.

    Parameters
    ----------
    da : xr.DataArray
        Input data array
    threshold : float
        Threshold value for identifying extreme events
    min_separation : int
        Minimum separation (in time steps) between independent events
    dim : str, default 'time'
        Dimension along which to decluster
    method : str, default 'runs'
        Declustering method:
        - 'runs': Select maximum from each exceedance run, merging runs
          whose below-threshold gap is shorter than ``min_separation``
        - 'separation': Enforce minimum time separation between peaks
    keep_time : bool, default False
        If True, attach a ``f"{dim}_of_max"`` coordinate holding the **original
        timestamp** at which each cluster's maximum occurred. The ``'runs'``
        method otherwise labels the representative at the *end* of its run, so
        this recovers the true peak time; for ``'separation'`` the peaks already
        sit at their original positions.

    Returns
    -------
    xr.DataArray
        Declustered block maxima with reduced temporal dependence.

    Examples
    --------
    >>> # Extract independent storm maxima (3-day separation)
    >>> storm_max = declustered_block_maxima(
    ...     da, threshold=100, min_separation=3, method='runs'
    ... )

    >>> # Heat wave maxima with 7-day separation
    >>> heatwave_max = declustered_block_maxima(
    ...     da, threshold=35, min_separation=7, method='separation'
    ... )

    Notes
    -----
    The runs method identifies continuous periods above threshold and takes
    the maximum from each run. The separation method ensures peaks are
    separated by at least min_separation time steps.
    """
    from xtremax.extraction.decluster import decluster_runs, decluster_separation

    if method == "runs":
        # Delegate so run IDs are computed per non-`dim` slice and cannot
        # collide across batch rows (e.g. different sites).
        # `min_separation` is the runs parameter r: below-threshold gaps
        # shorter than it are intra-cluster.
        reduced = decluster_runs(
            da,
            threshold=threshold,
            dim=dim,
            reduction="max",
            min_separation=min_separation,
        )
        if keep_time:
            reduced = _attach_runs_argtime(reduced, da, threshold, min_separation, dim)
        return reduced.dropna(dim, how="all")

    elif method == "separation":
        # Delegate to the separation-based declustering helper, which
        # actually applies `min_separation` in units of original time steps.
        peaks = decluster_separation(
            da, threshold=threshold, min_separation=min_separation, dim=dim
        )
        if keep_time:
            # Separation keeps peaks at their own positions, so the maximum's
            # coordinate is the position itself.
            arg = xr.full_like(peaks, np.nan, dtype=da[dim].dtype)
            arg = arg.where(peaks.isnull(), peaks[dim])
            peaks = peaks.assign_coords({f"{dim}_of_max": arg})
        return peaks

    else:
        raise ValueError(f"Unknown method: {method}. Use 'runs' or 'separation'")

Threshold selection

Peaks-over-threshold keeps every exceedance above \(u\), so the threshold choice is the bias–variance tradeoff: too low and the GPD asymptotics fail, too high and you run out of data. A constant threshold is only defensible for stationary data — under trends or seasonality, use the time-varying variants so the exceedance rate stays roughly constant.

Extraction utilities for extreme values on xarray data.

Block maxima, threshold selection, and declustering — the three standard preprocessing steps for extreme value analysis on spatiotemporal data.

constant_threshold(da: xr.DataArray, threshold: float) -> float

Return a constant threshold value.

This is the simplest threshold selection method, using a fixed value across all times and locations.

Parameters:

Name	Type	Description	Default
da	DataArray	Input data array (included for API consistency)	required
threshold	float	Constant threshold value	required

Returns:

Type	Description
float	The threshold value

Examples:

>>> # Use a constant threshold of 100
>>> thresh = constant_threshold(da, 100)

Source code in src/xtremax/extraction/threshold.py

def constant_threshold(da: xr.DataArray, threshold: float) -> float:
    """
    Return a constant threshold value.

    This is the simplest threshold selection method, using a fixed value
    across all times and locations.

    Parameters
    ----------
    da : xr.DataArray
        Input data array (included for API consistency)
    threshold : float
        Constant threshold value

    Returns
    -------
    float
        The threshold value

    Examples
    --------
    >>> # Use a constant threshold of 100
    >>> thresh = constant_threshold(da, 100)
    """
    return threshold

quantile_threshold(da: xr.DataArray, quantile: float, dim: Hashable | Sequence[Hashable] | None = None) -> float | xr.DataArray

Calculate threshold based on a percentile of the data.

Common choices are 0.95 (95th percentile) or 0.99 (99th percentile). Higher quantiles result in fewer but more extreme exceedances.

Parameters:

Name	Type	Description	Default
da	DataArray	Input data array	required
quantile	float	Quantile value between 0 and 1 (e.g., 0.95 for 95th percentile)	required
dim	Hashable or Sequence[Hashable]	Dimension(s) over which to compute quantile. If None, computes over all dimensions returning a scalar.	None

Returns:

Type	Description
float or DataArray	Threshold value(s)

Examples:

>>> # 95th percentile threshold
>>> thresh = quantile_threshold(da, 0.95)

>>> # 99th percentile, computed separately for each location
>>> thresh = quantile_threshold(da, 0.99, dim='time')

Source code in src/xtremax/extraction/threshold.py

def quantile_threshold(
    da: xr.DataArray,
    quantile: float,
    dim: Hashable | Sequence[Hashable] | None = None,
) -> float | xr.DataArray:
    """
    Calculate threshold based on a percentile of the data.

    Common choices are 0.95 (95th percentile) or 0.99 (99th percentile).
    Higher quantiles result in fewer but more extreme exceedances.

    Parameters
    ----------
    da : xr.DataArray
        Input data array
    quantile : float
        Quantile value between 0 and 1 (e.g., 0.95 for 95th percentile)
    dim : Hashable or Sequence[Hashable], optional
        Dimension(s) over which to compute quantile. If None, computes
        over all dimensions returning a scalar.

    Returns
    -------
    float or xr.DataArray
        Threshold value(s)

    Examples
    --------
    >>> # 95th percentile threshold
    >>> thresh = quantile_threshold(da, 0.95)

    >>> # 99th percentile, computed separately for each location
    >>> thresh = quantile_threshold(da, 0.99, dim='time')
    """
    # Drop the scalar `quantile` coordinate xarray attaches — it leaks
    # into downstream arithmetic (any later combination with another
    # quantile-carrying array would conflict on the coord).
    result = da.quantile(quantile, dim=dim).drop_vars("quantile")
    if result.ndim == 0:
        # Fully reduced: honour the documented `float` return type.
        return float(result)
    return result

rolling_threshold(da: xr.DataArray, quantile: float, window_size: int = 15, time_dim: str = 'time', center: bool = True, min_periods: int | None = None) -> xr.DataArray

Calculate time-varying threshold using a rolling window.

Uses a rolling window to compute a threshold that varies smoothly over time. Useful for removing short-term variations while preserving longer-term trends.

Parameters:

Name	Type	Description	Default
da	DataArray	Input data array with time dimension	required
quantile	float	Quantile value between 0 and 1	required
window_size	int	Size of rolling window (in time steps) for computing the threshold	15
time_dim	str	Name of the time dimension	'time'
center	bool	If True, set window labels at center of window	True
min_periods	int	Minimum number of observations required in window. If None, uses window_size. Windows with fewer valid observations (e.g. at the series edges) return NaN.	None

Returns:

Type	Description
DataArray	Time-varying threshold values with smooth temporal variation. Positions whose window holds fewer than min_periods valid observations are NaN.

Examples:

>>> # Rolling 15-day 95th percentile threshold
>>> thresh = rolling_threshold(da, 0.95, window_size=15)

>>> # Rolling 30-day 99th percentile threshold
>>> thresh = rolling_threshold(da, 0.99, window_size=30)

>>> # Non-centered window (past values only)
>>> thresh = rolling_threshold(da, 0.95, window_size=30, center=False)

Notes

Rolling window approach provides smoother temporal variation than groupby methods. Good for removing high-frequency variations while keeping the seasonal cycle.

Source code in src/xtremax/extraction/threshold.py

def rolling_threshold(
    da: xr.DataArray,
    quantile: float,
    window_size: int = 15,
    time_dim: str = "time",
    center: bool = True,
    min_periods: int | None = None,
) -> xr.DataArray:
    """
    Calculate time-varying threshold using a rolling window.

    Uses a rolling window to compute a threshold that varies smoothly
    over time. Useful for removing short-term variations while preserving
    longer-term trends.

    Parameters
    ----------
    da : xr.DataArray
        Input data array with time dimension
    quantile : float
        Quantile value between 0 and 1
    window_size : int, default 15
        Size of rolling window (in time steps) for computing the threshold
    time_dim : str, default 'time'
        Name of the time dimension
    center : bool, default True
        If True, set window labels at center of window
    min_periods : int, optional
        Minimum number of observations required in window. If None, uses
        window_size. Windows with fewer valid observations (e.g. at the
        series edges) return NaN.

    Returns
    -------
    xr.DataArray
        Time-varying threshold values with smooth temporal variation.
        Positions whose window holds fewer than ``min_periods`` valid
        observations are NaN.

    Examples
    --------
    >>> # Rolling 15-day 95th percentile threshold
    >>> thresh = rolling_threshold(da, 0.95, window_size=15)

    >>> # Rolling 30-day 99th percentile threshold
    >>> thresh = rolling_threshold(da, 0.99, window_size=30)

    >>> # Non-centered window (past values only)
    >>> thresh = rolling_threshold(da, 0.95, window_size=30, center=False)

    Notes
    -----
    Rolling window approach provides smoother temporal variation than
    groupby methods. Good for removing high-frequency variations while
    keeping the seasonal cycle.
    """
    if min_periods is None:
        min_periods = window_size

    rolling = da.rolling(
        {time_dim: window_size}, center=center, min_periods=min_periods
    )
    # xarray removed DataArrayRolling.quantile; materialize the window axis
    # with construct and take the quantile along it. `construct` does not
    # apply `min_periods` (and `.quantile` skips NaN), so partial edge
    # windows must be masked explicitly — otherwise they'd return
    # quantiles of only 1..window_size-1 samples.
    windowed = rolling.construct("_window")
    result = windowed.quantile(quantile, dim="_window").drop_vars("quantile")
    valid = windowed.notnull().sum("_window") >= min_periods
    return result.where(valid)

seasonal_threshold(da: xr.DataArray, quantile: float = 0.95, time_dim: str = 'time') -> xr.DataArray

Calculate seasonal threshold (backward compatibility wrapper).

This is a convenience wrapper around temporal_threshold with groupby='time.season'. For more control, use temporal_threshold directly.

Parameters:

Name	Type	Description	Default
da	DataArray	Input data array	required
quantile	float	Quantile value between 0 and 1	0.95
time_dim	str	Name of the time dimension	'time'

Returns:

Type	Description
DataArray	Seasonal threshold values

Examples:

>>> # 95th percentile by season
>>> thresh = seasonal_threshold(da, 0.95)

See Also

temporal_threshold : More flexible temporal grouping rolling_threshold : Smooth rolling window threshold

Source code in src/xtremax/extraction/threshold.py

def seasonal_threshold(
    da: xr.DataArray, quantile: float = 0.95, time_dim: str = "time"
) -> xr.DataArray:
    """
    Calculate seasonal threshold (backward compatibility wrapper).

    This is a convenience wrapper around temporal_threshold with
    groupby='time.season'. For more control, use temporal_threshold directly.

    Parameters
    ----------
    da : xr.DataArray
        Input data array
    quantile : float, default 0.95
        Quantile value between 0 and 1
    time_dim : str, default 'time'
        Name of the time dimension

    Returns
    -------
    xr.DataArray
        Seasonal threshold values

    Examples
    --------
    >>> # 95th percentile by season
    >>> thresh = seasonal_threshold(da, 0.95)

    See Also
    --------
    temporal_threshold : More flexible temporal grouping
    rolling_threshold : Smooth rolling window threshold
    """
    # Derive the groupby key from `time_dim` instead of hard-coding
    # "time.season"; otherwise datasets whose temporal axis is named
    # something else (e.g. "date") fail despite the API claiming support
    # via the `time_dim` argument.
    return temporal_threshold(
        da, quantile, groupby=f"{time_dim}.season", time_dim=time_dim
    )

temporal_threshold(da: xr.DataArray, quantile: float, groupby: str = 'time.season', time_dim: str = 'time', method: str = 'quantile') -> xr.DataArray

Calculate time-varying threshold based on temporal grouping.

Computes threshold separately for each temporal group (season, month, year, week, dayofyear, etc.) using groupby operations. More flexible than rolling window approach.

Parameters:

Name	Type	Description	Default
da	DataArray	Input data array with time dimension	required
quantile	float	Quantile value between 0 and 1 (e.g., 0.95 for 95th percentile)	required
groupby	str	Grouping specification. Common options: - 'time.season' : seasonal (DJF, MAM, JJA, SON) - 'time.month' : monthly (1-12) - 'time.year' : annual - 'time.dayofyear' : day of year (1-365/366) - 'time.hour' : hour of day (For week-of-year grouping, precompute an isocalendar week coordinate; the 'time.week' accessor is deprecated.)	'time.season'
time_dim	str	Name of the time dimension	'time'
method	str	Aggregation method: 'quantile', 'mean', 'median', 'max', 'min'	'quantile'

Returns:

Type	Description
DataArray	Time-varying threshold values, aligned with original time coordinate

Examples:

>>> # Seasonal 95th percentile thresholds
>>> thresh = temporal_threshold(da, 0.95, groupby='time.season')

>>> # Monthly 99th percentile thresholds
>>> thresh = temporal_threshold(da, 0.99, groupby='time.month')

>>> # Day-of-year thresholds (smooth seasonal cycle)
>>> thresh = temporal_threshold(da, 0.95, groupby='time.dayofyear')

>>> # Annual maximum
>>> thresh = temporal_threshold(da, 0.95, groupby='time.year', method='max')

>>> # Day-of-year mean thresholds
>>> thresh = temporal_threshold(da, 0.9, groupby='time.dayofyear', method='mean')

Notes

The function groups data by the specified temporal accessor (e.g., month), computes the threshold for each group, then broadcasts back to the original time coordinate. This is useful for capturing regular temporal patterns in the tail distribution.

For day-of-year grouping, this gives a smooth seasonal cycle of thresholds. For month/season, this gives discrete jumps between periods.

Source code in src/xtremax/extraction/threshold.py

def temporal_threshold(
    da: xr.DataArray,
    quantile: float,
    groupby: str = "time.season",
    time_dim: str = "time",
    method: str = "quantile",
) -> xr.DataArray:
    """
    Calculate time-varying threshold based on temporal grouping.

    Computes threshold separately for each temporal group (season, month,
    year, week, dayofyear, etc.) using groupby operations. More flexible
    than rolling window approach.

    Parameters
    ----------
    da : xr.DataArray
        Input data array with time dimension
    quantile : float
        Quantile value between 0 and 1 (e.g., 0.95 for 95th percentile)
    groupby : str, default 'time.season'
        Grouping specification. Common options:
        - 'time.season' : seasonal (DJF, MAM, JJA, SON)
        - 'time.month' : monthly (1-12)
        - 'time.year' : annual
        - 'time.dayofyear' : day of year (1-365/366)
        - 'time.hour' : hour of day
        (For week-of-year grouping, precompute an isocalendar week
        coordinate; the ``'time.week'`` accessor is deprecated.)
    time_dim : str, default 'time'
        Name of the time dimension
    method : str, default 'quantile'
        Aggregation method: 'quantile', 'mean', 'median', 'max', 'min'

    Returns
    -------
    xr.DataArray
        Time-varying threshold values, aligned with original time coordinate

    Examples
    --------
    >>> # Seasonal 95th percentile thresholds
    >>> thresh = temporal_threshold(da, 0.95, groupby='time.season')

    >>> # Monthly 99th percentile thresholds
    >>> thresh = temporal_threshold(da, 0.99, groupby='time.month')

    >>> # Day-of-year thresholds (smooth seasonal cycle)
    >>> thresh = temporal_threshold(da, 0.95, groupby='time.dayofyear')

    >>> # Annual maximum
    >>> thresh = temporal_threshold(da, 0.95, groupby='time.year', method='max')

    >>> # Day-of-year mean thresholds
    >>> thresh = temporal_threshold(da, 0.9, groupby='time.dayofyear', method='mean')

    Notes
    -----
    The function groups data by the specified temporal accessor (e.g., month),
    computes the threshold for each group, then broadcasts back to the original
    time coordinate. This is useful for capturing regular temporal patterns in
    the tail distribution.

    For day-of-year grouping, this gives a smooth seasonal cycle of thresholds.
    For month/season, this gives discrete jumps between periods.
    """
    # Group by temporal accessor
    grouped = da.groupby(groupby)

    # Compute threshold for each group
    if method == "quantile":
        threshold_by_group = grouped.quantile(quantile, dim=time_dim).drop_vars(
            "quantile"
        )
    elif method == "mean":
        threshold_by_group = grouped.mean(dim=time_dim)
    elif method == "median":
        threshold_by_group = grouped.median(dim=time_dim)
    elif method == "max":
        threshold_by_group = grouped.max(dim=time_dim)
    elif method == "min":
        threshold_by_group = grouped.min(dim=time_dim)
    else:
        raise ValueError(
            f"Unknown method: {method}. "
            f"Choose from: 'quantile', 'mean', 'median', 'max', 'min'"
        )

    # Broadcast threshold back to original time coordinate
    # This ensures threshold has same shape as input
    threshold_expanded = threshold_by_group.sel({groupby.split(".")[-1]: da[groupby]})

    return threshold_expanded

Declustering & the extremal index

Raw exceedances arrive in clusters — one storm produces several consecutive threshold crossings — which violates the independence the GPD likelihood assumes. Declustering reduces each cluster to one peak; the extremal index \(\theta \in (0, 1]\) quantifies how much clustering there was (\(\theta = 1\) means none).

Extraction utilities for extreme values on xarray data.

Block maxima, threshold selection, and declustering — the three standard preprocessing steps for extreme value analysis on spatiotemporal data.

decluster_runs(da: xr.DataArray, threshold: float | xr.DataArray, dim: str = 'time', reduction: str = 'max', min_separation: int = 1) -> xr.DataArray

Decluster exceedances using the runs method.

Identifies runs of exceedances and applies a reduction function to each run. This ensures temporal independence by treating each cluster as a single event.

Parameters:

Name	Type	Description	Default
da	DataArray	Input data array	required
threshold	float or DataArray	Threshold value(s)	required
dim	str	Dimension along which to identify runs	'time'
reduction	str	Reduction function to apply to each run: 'max', 'mean', 'sum', 'min'	'max'
min_separation	int	Runs-declustering separation parameter r: two exceedances belong to the same cluster when fewer than min_separation consecutive below-threshold steps separate them. The default 1 reproduces classical contiguous runs (any single dip below the threshold ends the run).	1

Returns:

Type	Description
DataArray	Array with the same shape as da. Run representatives are placed at the last exceedance position of each run along dim and all other positions are NaN. For multi-dimensional inputs, each non-dim slice is declustered independently (so run IDs never collide across slices).

Examples:

>>> # Maximum value from each exceedance cluster
>>> declustered = decluster_runs(da, threshold=100, reduction='max')

>>> # Merge exceedances separated by gaps of fewer than 3 steps
>>> declustered = decluster_runs(da, threshold=100, min_separation=3)

Notes

The runs method identifies periods where values exceed the threshold, merging exceedances whose below-threshold gap is shorter than min_separation into one cluster, and extracts one representative value per cluster. The reduction is applied to the cluster's exceedance values only — intra-cluster below-threshold samples do not contribute.

Source code in src/xtremax/extraction/decluster.py

def decluster_runs(
    da: xr.DataArray,
    threshold: float | xr.DataArray,
    dim: str = "time",
    reduction: str = "max",
    min_separation: int = 1,
) -> xr.DataArray:
    """
    Decluster exceedances using the runs method.

    Identifies runs of exceedances and applies a reduction function to
    each run. This ensures temporal independence by treating each
    cluster as a single event.

    Parameters
    ----------
    da : xr.DataArray
        Input data array
    threshold : float or xr.DataArray
        Threshold value(s)
    dim : str, default 'time'
        Dimension along which to identify runs
    reduction : str, default 'max'
        Reduction function to apply to each run: 'max', 'mean', 'sum', 'min'
    min_separation : int, default 1
        Runs-declustering separation parameter *r*: two exceedances
        belong to the same cluster when fewer than ``min_separation``
        consecutive below-threshold steps separate them. The default 1
        reproduces classical contiguous runs (any single dip below the
        threshold ends the run).

    Returns
    -------
    xr.DataArray
        Array with the same shape as ``da``. Run representatives are placed
        at the last exceedance position of each run along ``dim`` and all
        other positions are NaN. For multi-dimensional inputs, each
        non-``dim`` slice is declustered independently (so run IDs never
        collide across slices).

    Examples
    --------
    >>> # Maximum value from each exceedance cluster
    >>> declustered = decluster_runs(da, threshold=100, reduction='max')

    >>> # Merge exceedances separated by gaps of fewer than 3 steps
    >>> declustered = decluster_runs(da, threshold=100, min_separation=3)

    Notes
    -----
    The runs method identifies periods where values exceed the
    threshold, merging exceedances whose below-threshold gap is shorter
    than ``min_separation`` into one cluster, and extracts one
    representative value per cluster. The reduction is applied to the
    cluster's *exceedance* values only — intra-cluster below-threshold
    samples do not contribute.
    """
    if reduction not in _REDUCERS:
        raise ValueError(f"Unknown reduction: {reduction}")
    if min_separation < 1:
        raise ValueError(f"min_separation must be >= 1, got {min_separation}")
    reducer = _REDUCERS[reduction]

    exceedances = da > threshold

    def _decluster_runs_1d(values: np.ndarray, exc: np.ndarray) -> np.ndarray:
        """Reduce each run of exceedances to a single value.

        The representative is stored at the last exceedance position of
        its run; all other entries are NaN. Operating per 1-D slice keeps
        run IDs from different batch rows (e.g. different ``site``s) from
        colliding.
        """
        out = np.full(values.shape[0], np.nan, dtype=float)
        positions = np.flatnonzero(exc)
        if positions.size == 0:
            return out
        # A gap of >= min_separation below-threshold steps between
        # consecutive exceedances starts a new cluster.
        breaks = np.flatnonzero(np.diff(positions) - 1 >= min_separation) + 1
        for cluster in np.split(positions, breaks):
            out[cluster[-1]] = reducer(values[cluster])
        return out

    return xr.apply_ufunc(
        _decluster_runs_1d,
        da,
        exceedances,
        input_core_dims=[[dim], [dim]],
        output_core_dims=[[dim]],
        vectorize=True,
        output_dtypes=[float],
    )

decluster_separation(da: xr.DataArray, threshold: float | xr.DataArray, min_separation: int, dim: str = 'time') -> xr.DataArray

Decluster by enforcing minimum separation between peaks.

Identifies peaks above threshold and ensures they are separated by at least min_separation time steps. When peaks are too close, keeps the larger one.

Parameters:

Name	Type	Description	Default
da	DataArray	Input data array	required
threshold	float or DataArray	Threshold value(s)	required
min_separation	int	Minimum separation (in time steps) between independent peaks	required
dim	str	Dimension along which to enforce separation	'time'

Returns:

Type	Description
DataArray	Declustered peak values

Examples:

>>> # Peaks separated by at least 3 days
>>> declustered = decluster_separation(da, threshold=100, min_separation=3)

>>> # Storm events with 7-day minimum separation
>>> storms = decluster_separation(da, threshold=95, min_separation=7)

Notes

This method first identifies all local maxima above threshold, then iteratively removes peaks that are too close to larger peaks.

Plateaus (runs of equal consecutive values) count as a single peak when the values adjacent to the entire run are strictly lower — the peak is stamped at the run's first sample. Quantized observational data (rounded gauges, integer records) hit this case routinely; requiring strict inequality between adjacent samples used to drop such events entirely, while a naive adjacent-sample tie-break would also mark rising shelves (e.g. [0, 5, 5, 6]) as peaks.

Source code in src/xtremax/extraction/decluster.py

def decluster_separation(
    da: xr.DataArray,
    threshold: float | xr.DataArray,
    min_separation: int,
    dim: str = "time",
) -> xr.DataArray:
    """
    Decluster by enforcing minimum separation between peaks.

    Identifies peaks above threshold and ensures they are separated by
    at least min_separation time steps. When peaks are too close, keeps
    the larger one.

    Parameters
    ----------
    da : xr.DataArray
        Input data array
    threshold : float or xr.DataArray
        Threshold value(s)
    min_separation : int
        Minimum separation (in time steps) between independent peaks
    dim : str, default 'time'
        Dimension along which to enforce separation

    Returns
    -------
    xr.DataArray
        Declustered peak values

    Examples
    --------
    >>> # Peaks separated by at least 3 days
    >>> declustered = decluster_separation(da, threshold=100, min_separation=3)

    >>> # Storm events with 7-day minimum separation
    >>> storms = decluster_separation(da, threshold=95, min_separation=7)

    Notes
    -----
    This method first identifies all local maxima above threshold,
    then iteratively removes peaks that are too close to larger peaks.

    Plateaus (runs of equal consecutive values) count as a single peak
    when the values adjacent to the entire run are strictly lower — the
    peak is stamped at the run's first sample. Quantized observational
    data (rounded gauges, integer records) hit this case routinely;
    requiring strict inequality between adjacent samples used to drop
    such events entirely, while a naive adjacent-sample tie-break would
    also mark rising shelves (e.g. ``[0, 5, 5, 6]``) as peaks.
    """
    # Identify local maxima above threshold.
    #
    # `da.shift` fills the edge with NaN, so comparisons at the first /
    # last positions along `dim` silently evaluate to False. Previously
    # this meant boundary elements could NEVER be marked as peaks — a
    # silent event drop for short windows or chunked workflows. Treat a
    # missing neighbour as "no constraint from that side": a boundary
    # position is a peak if it is strictly greater than its single
    # existing neighbour (and above threshold).
    # Plateau-aware peak detection (equal consecutive maxima): group the
    # series into maximal runs of equal values; a run is a peak iff the
    # values adjacent to the *entire* run are strictly lower (a missing
    # neighbour at the series boundary is no constraint). The peak is
    # stamped at the run's first sample. Comparing only adjacent samples
    # either dropped flat-topped exceedances entirely (strict on both
    # sides) or marked rising shelves like [0, 5, 5, 6] as peaks
    # (>= on one side).
    exceedances = da > threshold

    def _plateau_peaks_1d(values: np.ndarray, exc: np.ndarray) -> np.ndarray:
        n = values.shape[0]
        peaks = np.zeros(n, dtype=bool)
        if n == 0:
            return peaks
        # Run starts: first index, then every value change. NaN != NaN,
        # so NaN samples form their own (never-exceeding) runs.
        change = np.concatenate(([True], values[1:] != values[:-1]))
        starts = np.flatnonzero(change)
        ends = np.append(starts[1:] - 1, n - 1)
        left_ok = (starts == 0) | (values[np.maximum(starts - 1, 0)] < values[starts])
        right_ok = (ends == n - 1) | (
            values[np.minimum(ends + 1, n - 1)] < values[starts]
        )
        peak_starts = starts[left_ok & right_ok & exc[starts]]
        peaks[peak_starts] = True
        return peaks

    is_peak = xr.apply_ufunc(
        _plateau_peaks_1d,
        da,
        exceedances,
        input_core_dims=[[dim], [dim]],
        output_core_dims=[[dim]],
        vectorize=True,
        output_dtypes=[bool],
    )

    def _select_separated_peaks_1d(
        values: np.ndarray, peak_mask: np.ndarray
    ) -> np.ndarray:
        """Greedily pick peaks so retained maxima are ≥ min_separation apart.

        Works purely in positional indices along the core axis, which is
        what `min_separation` is defined against (time steps, not
        elements of the filtered peak list).
        """
        peak_positions = np.flatnonzero(peak_mask)
        selected = np.zeros_like(peak_mask, dtype=bool)
        if peak_positions.size == 0:
            return selected
        peak_values = values[peak_positions]
        order = np.argsort(peak_values)[::-1]
        chosen: list[int] = []
        for idx in peak_positions[order]:
            pos = int(idx)
            if all(abs(pos - s) >= min_separation for s in chosen):
                chosen.append(pos)
                selected[pos] = True
        return selected

    # `apply_ufunc` with core dim == `dim` vectorises the 1-D selector
    # over any remaining batch dims (e.g. `site`), so each batch row
    # gets its own separation-filtered peak set and peak indices stay
    # valid positions along `dim`.
    selected_mask = xr.apply_ufunc(
        _select_separated_peaks_1d,
        da,
        is_peak,
        input_core_dims=[[dim], [dim]],
        output_core_dims=[[dim]],
        vectorize=True,
        output_dtypes=[bool],
    )

    return da.where(selected_mask, drop=True)

estimate_extremal_index(da: xr.DataArray, threshold: float | xr.DataArray, dim: str = 'time', method: str = 'runs', min_separation: int = 1) -> float | xr.DataArray

Estimate the extremal index (theta), a measure of clustering.

The extremal index θ ∈ [0, 1] measures the degree of clustering in extremes: - θ = 1: no clustering (independent extremes) - θ < 1: clustering present (dependent extremes) - θ = 0: complete clustering (perfect dependence)

Parameters:

Name	Type	Description	Default
da	DataArray	Input data array	required
threshold	float or DataArray	Threshold value for identifying extremes	required
dim	str	Dimension along which to estimate clustering	'time'
method	str	Estimation method: 'runs' or 'intervals'	'runs'
min_separation	int	Runs parameter r: exceedances separated by fewer than min_separation consecutive below-threshold steps count as the same cluster. The default 1 reproduces the classical contiguous-runs estimator.	1

Returns:

Type	Description
float or DataArray	Extremal index estimate (between 0 and 1). For 1-D input a scalar float. For inputs with additional dimensions, a DataArray of per-slice estimates (one θ per batch row, e.g. per site).

Examples:

>>> # Estimate clustering strength
>>> theta = estimate_extremal_index(da, threshold=100)
>>> print(f"Extremal index: {theta:.3f}")
>>> # theta ≈ 0.7 suggests moderate clustering

>>> # Use high threshold (99th percentile)
>>> from xtremax.extraction.threshold import quantile_threshold
>>> thresh = quantile_threshold(da, 0.99)
>>> theta = estimate_extremal_index(da, thresh)

Notes

The runs estimator uses: θ = (number of clusters) / (number of exceedances) A smaller θ indicates more clustering. Counting is performed independently per non-dim slice so cluster labels from different slices cannot collide.

Source code in src/xtremax/extraction/decluster.py

def estimate_extremal_index(
    da: xr.DataArray,
    threshold: float | xr.DataArray,
    dim: str = "time",
    method: str = "runs",
    min_separation: int = 1,
) -> float | xr.DataArray:
    """
    Estimate the extremal index (theta), a measure of clustering.

    The extremal index θ ∈ [0, 1] measures the degree of clustering in extremes:
    - θ = 1: no clustering (independent extremes)
    - θ < 1: clustering present (dependent extremes)
    - θ = 0: complete clustering (perfect dependence)

    Parameters
    ----------
    da : xr.DataArray
        Input data array
    threshold : float or xr.DataArray
        Threshold value for identifying extremes
    dim : str, default 'time'
        Dimension along which to estimate clustering
    method : str, default 'runs'
        Estimation method: 'runs' or 'intervals'
    min_separation : int, default 1
        Runs parameter *r*: exceedances separated by fewer than
        ``min_separation`` consecutive below-threshold steps count as
        the same cluster. The default 1 reproduces the classical
        contiguous-runs estimator.

    Returns
    -------
    float or xr.DataArray
        Extremal index estimate (between 0 and 1). For 1-D input a scalar
        ``float``. For inputs with additional dimensions, a DataArray of
        per-slice estimates (one θ per batch row, e.g. per site).

    Examples
    --------
    >>> # Estimate clustering strength
    >>> theta = estimate_extremal_index(da, threshold=100)
    >>> print(f"Extremal index: {theta:.3f}")
    >>> # theta ≈ 0.7 suggests moderate clustering

    >>> # Use high threshold (99th percentile)
    >>> from xtremax.extraction.threshold import quantile_threshold
    >>> thresh = quantile_threshold(da, 0.99)
    >>> theta = estimate_extremal_index(da, thresh)

    Notes
    -----
    The runs estimator uses: θ = (number of clusters) / (number of exceedances)
    A smaller θ indicates more clustering. Counting is performed
    independently per non-``dim`` slice so cluster labels from different
    slices cannot collide.
    """
    if method == "intervals":
        raise NotImplementedError("Intervals method not yet implemented")
    if method != "runs":
        raise ValueError(f"Unknown method: {method}")
    if min_separation < 1:
        raise ValueError(f"min_separation must be >= 1, got {min_separation}")

    exceedances = da > threshold

    def _theta_runs_1d(exc: np.ndarray) -> float:
        positions = np.flatnonzero(exc)
        if positions.size == 0:
            return np.nan
        # Count clusters: a new cluster starts whenever the gap of
        # consecutive below-threshold steps reaches min_separation.
        n_runs = 1 + int(np.sum(np.diff(positions) - 1 >= min_separation))
        return n_runs / positions.size

    theta = xr.apply_ufunc(
        _theta_runs_1d,
        exceedances,
        input_core_dims=[[dim]],
        output_core_dims=[[]],
        vectorize=True,
        output_dtypes=[float],
    )
    # Preserve the historical 1-D scalar return type.
    if theta.ndim == 0:
        return float(theta.values)
    return theta

Quantile-regression thresholds

A threshold that follows a covariate-driven trend, fitted by quantile regression rather than taken as a fixed quantile — the right tool when sea level rise or warming shifts the whole distribution under you.

These require scikit-learn and are only exported when the xtremax[threshold] extra is installed:

pip install "xtremax[threshold]"

Quantile regression utilities for extremes thresholding.

Provides two layers:

1. XarrayQuantileRegressor — a thin sklearn-compatible wrapper around sklearn.linear_model.QuantileRegressor that accepts and returns xr.DataArray objects, preserving coordinates through fit → predict.

2. quantile_regression_threshold — domain-specific function that builds the feature matrix (numeric time, optional covariates) and delegates to the regressor, returning a time-varying threshold suitable for peak-over-threshold analysis.

Example

import xarray as xr, numpy as np, pandas as pd da = xr.DataArray( ... np.random.default_rng(0).standard_normal(200), ... dims="time", ... coords={"time": pd.date_range("2000", periods=200, freq="D")}, ... ) threshold = quantile_regression_threshold(da, quantile=0.95, time_dim="time")

XarrayQuantileRegressor(quantile: float = 0.5, alpha: float = 0.0, solver: str = 'highs', fit_intercept: bool = True)

Bases: BaseEstimator, RegressorMixin

Quantile regressor that operates on :class:xr.DataArray inputs.

Thin wrapper around :class:sklearn.linear_model.QuantileRegressor. Both fit and predict expect DataArray objects; coordinate metadata on the leading dimension is preserved through the round-trip.

Parameters:

Name	Type	Description	Default
quantile	float	Target quantile ∈ (0, 1).	0.5
alpha	float	L1 regularisation strength. 0 gives the classical unregularised linear quantile regression (matching the default behaviour of statsmodels.QuantReg).	0.0
solver	str	Solver forwarded to QuantileRegressor. "highs" (the default) is robust for most problem sizes.	'highs'
fit_intercept	bool	Whether the sklearn estimator should fit its own intercept. Set to False if you supply a constant column yourself.	True

Source code in src/xtremax/extraction/quantile_regression.py

def __init__(
    self,
    quantile: float = 0.5,
    alpha: float = 0.0,
    solver: str = "highs",
    fit_intercept: bool = True,
):
    self.quantile = quantile
    self.alpha = alpha
    self.solver = solver
    self.fit_intercept = fit_intercept

coef_: np.ndarray property

Regression coefficients (delegates to the inner model).

intercept_: float property

Intercept (delegates to the inner model).

fit(X: xr.DataArray, y: xr.DataArray) -> XarrayQuantileRegressor

Fit the quantile regression model.

Parameters:

Name	Type	Description	Default
X	(DataArray, shape(n_samples) or (n_samples, n_features))	Feature matrix.	required
y	(DataArray, shape(n_samples))	Target values.	required

Returns:

Type	Description
self

Source code in src/xtremax/extraction/quantile_regression.py

def fit(self, X: xr.DataArray, y: xr.DataArray) -> XarrayQuantileRegressor:
    """Fit the quantile regression model.

    Parameters
    ----------
    X : xr.DataArray, shape (n_samples,) or (n_samples, n_features)
        Feature matrix.
    y : xr.DataArray, shape (n_samples,)
        Target values.

    Returns
    -------
    self
    """
    X_np = X.values
    if X_np.ndim == 1:
        X_np = X_np.reshape(-1, 1)

    self.model_ = _SklearnQuantileRegressor(
        quantile=self.quantile,
        alpha=self.alpha,
        solver=self.solver,
        fit_intercept=self.fit_intercept,
    )
    self.model_.fit(X_np, y.values.ravel())
    self.n_features_in_ = X_np.shape[1]
    return self

predict(X: xr.DataArray) -> xr.DataArray

Predict quantile values.

Parameters:

Name	Type	Description	Default
X	DataArray	Feature matrix (same schema as training X).	required

Returns:

Type	Description
DataArray	Predicted quantile values with the leading dimension's coordinates preserved from X.

Source code in src/xtremax/extraction/quantile_regression.py

def predict(self, X: xr.DataArray) -> xr.DataArray:
    """Predict quantile values.

    Parameters
    ----------
    X : xr.DataArray
        Feature matrix (same schema as training *X*).

    Returns
    -------
    xr.DataArray
        Predicted quantile values with the leading dimension's
        coordinates preserved from *X*.
    """
    check_is_fitted(self, "model_")

    first_dim = X.dims[0]
    X_np = X.values
    if X_np.ndim == 1:
        X_np = X_np.reshape(-1, 1)

    y_hat = self.model_.predict(X_np)
    return xr.DataArray(
        y_hat,
        dims=[first_dim],
        coords={first_dim: X.coords[first_dim]},
    )

quantile_regression_threshold(da: xr.DataArray, quantile: float, time_dim: str, covariates: xr.DataArray | None = None, *, alpha: float = 0.0, solver: str = 'highs') -> xr.DataArray

Estimate a time-varying threshold via linear quantile regression.

Fits the model

.. math::

Q_{\tau}(y \mid t, \mathbf{z}) = \beta_0
    + \beta_1 t + \boldsymbol{\beta}_z^\top \mathbf{z}

where t is numeric time and z are optional covariates, then returns the fitted quantile surface evaluated at the observed times. The result is suitable as a non-stationary threshold for peak-over-threshold (POT) extreme value analysis.

Parameters:

Name	Type	Description	Default
da	DataArray	Response variable (1-D along time_dim).	required
quantile	float	Target quantile τ ∈ (0, 1), e.g. 0.95 for the 95th percentile.	required
time_dim	str	Name of the time dimension in da.	required
covariates	DataArray	Additional covariates aligned along time_dim.	None
alpha	float	L1 regularisation passed to the regressor (0 = unregularised).	0.0
solver	str	LP solver, forwarded to QuantileRegressor.	'highs'

Returns:

Type	Description
DataArray	Fitted threshold values with the same coords as da.

Examples:

>>> # 95th-percentile linear-trend threshold
>>> threshold = quantile_regression_threshold(da, 0.95, "time")

>>> # With a covariate (e.g., sea-surface temperature)
>>> threshold = quantile_regression_threshold(
...     da, 0.99, "time", covariates=sst
... )

Source code in src/xtremax/extraction/quantile_regression.py

def quantile_regression_threshold(
    da: xr.DataArray,
    quantile: float,
    time_dim: str,
    covariates: xr.DataArray | None = None,
    *,
    alpha: float = 0.0,
    solver: str = "highs",
) -> xr.DataArray:
    """Estimate a time-varying threshold via linear quantile regression.

    Fits the model

    .. math::

        Q_{\\tau}(y \\mid t, \\mathbf{z}) = \\beta_0
            + \\beta_1 t + \\boldsymbol{\\beta}_z^\\top \\mathbf{z}

    where *t* is numeric time and **z** are optional covariates, then
    returns the fitted quantile surface evaluated at the observed times.
    The result is suitable as a non-stationary threshold for
    peak-over-threshold (POT) extreme value analysis.

    Parameters
    ----------
    da : xr.DataArray
        Response variable (1-D along *time_dim*).
    quantile : float
        Target quantile τ ∈ (0, 1), e.g. 0.95 for the 95th percentile.
    time_dim : str
        Name of the time dimension in *da*.
    covariates : xr.DataArray, optional
        Additional covariates aligned along *time_dim*.
    alpha : float
        L1 regularisation passed to the regressor (0 = unregularised).
    solver : str
        LP solver, forwarded to ``QuantileRegressor``.

    Returns
    -------
    xr.DataArray
        Fitted threshold values with the same coords as *da*.

    Examples
    --------
    >>> # 95th-percentile linear-trend threshold
    >>> threshold = quantile_regression_threshold(da, 0.95, "time")

    >>> # With a covariate (e.g., sea-surface temperature)
    >>> threshold = quantile_regression_threshold(
    ...     da, 0.99, "time", covariates=sst
    ... )
    """
    X = _build_feature_matrix(da, time_dim, covariates)

    reg = XarrayQuantileRegressor(
        quantile=quantile,
        alpha=alpha,
        solver=solver,
        fit_intercept=True,
    )
    reg.fit(X, da)
    return reg.predict(X)