Metrics¶

xrops.metrics groups evaluation metrics by scientific diagnostic family. Each family pairs Layer-0 xarray functions with Layer-1 Operator wrappers that compose into Sequential pipelines and Graph networks.

For the design rationale, see the validation design doc and the validation API map.

Diagnostic families¶

Submodule	Family	Status
`metrics.pixel`	Pointwise scalar errors	shipped
`metrics.spectral`	Power-spectrum scores, resolved scale	shipped
`metrics.multiscale`	Region- and band-conditioned skill	shipped
`metrics.physical`	Physical-balance / conservation residuals	shipped
`metrics.masked`	Coverage-aware masked wrappers	shipped
`metrics.distributional`	Distributional distances (CRPS, Wasserstein, …)	shipped
`metrics.structural`	Structural / geometric similarity (SSIM, …)	partial
`metrics.probabilistic`	Ensemble calibration	partial
`metrics.object`	Event / object verification (POD, FAR, IoU)	partial
`metrics.forecast`	Lead-time skill	stub
`metrics.lagrangian`	Trajectory / transport metrics	stub

Ergonomic re-exports

The Layer-1 operators are re-exported flat from the package root, and a dedicated xrops.metrics.operators module mirrors them:

from xrops.metrics import RMSE, MSE, MAE, NRMSE, Bias, Correlation, R2Score, PSDScore
# equivalently:
from xrops.metrics.operators import RMSE, PSDScore

Pixel — pointwise errors¶

Per-pixel scalar error and skill: MSE, RMSE, MAE, bias, normalized RMSE (and its skill-score form), correlation, and R².

MSE ¶

MSE(variable: str, dims: Dim)

Bases: _PixelMetricOp

Mean-squared-error operator.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/pixel.py

def __init__(self, variable: str, dims: Dim):
    self.variable = variable
    self.dims = dims if isinstance(dims, str) else list(dims)

RMSE ¶

RMSE(variable: str, dims: Dim)

Bases: _PixelMetricOp

Root-mean-squared-error operator.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/pixel.py

def __init__(self, variable: str, dims: Dim):
    self.variable = variable
    self.dims = dims if isinstance(dims, str) else list(dims)

MAE ¶

MAE(variable: str, dims: Dim)

Bases: _PixelMetricOp

Mean-absolute-error operator.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/pixel.py

def __init__(self, variable: str, dims: Dim):
    self.variable = variable
    self.dims = dims if isinstance(dims, str) else list(dims)

Bias ¶

Bias(variable: str, dims: Dim)

Bases: _PixelMetricOp

Mean-bias operator (<pred - ref>).

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/pixel.py

def __init__(self, variable: str, dims: Dim):
    self.variable = variable
    self.dims = dims if isinstance(dims, str) else list(dims)

NRMSE ¶

NRMSE(variable: str, dims: Dim)

Bases: _PixelMetricOp

Normalized-RMSE operator (1 - RMSE / sqrt(<ref^2>)).

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/pixel.py

def __init__(self, variable: str, dims: Dim):
    self.variable = variable
    self.dims = dims if isinstance(dims, str) else list(dims)

NRMSEScore ¶

NRMSEScore(variable: str, dims: Dim)

Bases: _PixelMetricOp

Mercator/OceanBench NRMSE skill operator (1 - RMSE / std(ref)).

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/pixel.py

def __init__(self, variable: str, dims: Dim):
    self.variable = variable
    self.dims = dims if isinstance(dims, str) else list(dims)

Correlation ¶

Correlation(variable: str, dims: Dim)

Bases: _PixelMetricOp

Pearson-correlation operator.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/pixel.py

def __init__(self, variable: str, dims: Dim):
    self.variable = variable
    self.dims = dims if isinstance(dims, str) else list(dims)

R2Score ¶

R2Score(variable: str, dims: Dim)

Bases: _PixelMetricOp

Coefficient-of-determination operator (1 - SS_res / SS_tot).

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/pixel.py

def __init__(self, variable: str, dims: Dim):
    self.variable = variable
    self.dims = dims if isinstance(dims, str) else list(dims)

Functions¶

pixel ¶

Pointwise (pixel-level) evaluation metrics — public re-export.

Layer-0 functions: :func:mse, :func:rmse, :func:nrmse, :func:nrmse_score, :func:mae, :func:bias, :func:correlation, :func:r2_score.

Layer-1 operators: :class:MSE, :class:RMSE, :class:NRMSE, :class:NRMSEScore, :class:MAE, :class:Bias, :class:Correlation, :class:R2Score.

Implementation lives in :mod:xrops.metrics._src.pixel.

MAE ¶

MAE(variable: str, dims: Dim)

Bases: _PixelMetricOp

Mean-absolute-error operator.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/pixel.py

def __init__(self, variable: str, dims: Dim):
    self.variable = variable
    self.dims = dims if isinstance(dims, str) else list(dims)

MSE ¶

MSE(variable: str, dims: Dim)

Bases: _PixelMetricOp

Mean-squared-error operator.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/pixel.py

def __init__(self, variable: str, dims: Dim):
    self.variable = variable
    self.dims = dims if isinstance(dims, str) else list(dims)

NRMSE ¶

NRMSE(variable: str, dims: Dim)

Bases: _PixelMetricOp

Normalized-RMSE operator (1 - RMSE / sqrt(<ref^2>)).

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/pixel.py

def __init__(self, variable: str, dims: Dim):
    self.variable = variable
    self.dims = dims if isinstance(dims, str) else list(dims)

RMSE ¶

RMSE(variable: str, dims: Dim)

Bases: _PixelMetricOp

Root-mean-squared-error operator.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/pixel.py

def __init__(self, variable: str, dims: Dim):
    self.variable = variable
    self.dims = dims if isinstance(dims, str) else list(dims)

Bias ¶

Bias(variable: str, dims: Dim)

Bases: _PixelMetricOp

Mean-bias operator (<pred - ref>).

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/pixel.py

def __init__(self, variable: str, dims: Dim):
    self.variable = variable
    self.dims = dims if isinstance(dims, str) else list(dims)

Correlation ¶

Correlation(variable: str, dims: Dim)

Bases: _PixelMetricOp

Pearson-correlation operator.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/pixel.py

def __init__(self, variable: str, dims: Dim):
    self.variable = variable
    self.dims = dims if isinstance(dims, str) else list(dims)

NRMSEScore ¶

NRMSEScore(variable: str, dims: Dim)

Bases: _PixelMetricOp

Mercator/OceanBench NRMSE skill operator (1 - RMSE / std(ref)).

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/pixel.py

def __init__(self, variable: str, dims: Dim):
    self.variable = variable
    self.dims = dims if isinstance(dims, str) else list(dims)

R2Score ¶

R2Score(variable: str, dims: Dim)

Bases: _PixelMetricOp

Coefficient-of-determination operator (1 - SS_res / SS_tot).

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/pixel.py

def __init__(self, variable: str, dims: Dim):
    self.variable = variable
    self.dims = dims if isinstance(dims, str) else list(dims)

bias ¶

bias(pred: DataArray, ref: DataArray, *, dim: Dim) -> DataArray

Mean bias <pred - ref> reduced over dim.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/pixel.py

def bias(pred: xr.DataArray, ref: xr.DataArray, *, dim: Dim) -> xr.DataArray:
    """Mean bias ``<pred - ref>`` reduced over ``dim``."""
    return _apply_pixel_kernel(_kernels.bias, pred, ref, dim)

correlation ¶

correlation(pred: DataArray, ref: DataArray, *, dim: Dim) -> DataArray

Pearson correlation between prediction and reference over dim.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/pixel.py

def correlation(pred: xr.DataArray, ref: xr.DataArray, *, dim: Dim) -> xr.DataArray:
    """Pearson correlation between prediction and reference over ``dim``."""
    return _apply_pixel_kernel(_kernels.correlation, pred, ref, dim)

mae ¶

mae(pred: DataArray, ref: DataArray, *, dim: Dim) -> DataArray

Mean absolute error reduced over dim.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/pixel.py

def mae(pred: xr.DataArray, ref: xr.DataArray, *, dim: Dim) -> xr.DataArray:
    """Mean absolute error reduced over ``dim``."""
    return _apply_pixel_kernel(_kernels.mae, pred, ref, dim)

mse ¶

mse(pred: DataArray, ref: DataArray, *, dim: Dim) -> DataArray

Mean squared error reduced over dim.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/pixel.py

def mse(pred: xr.DataArray, ref: xr.DataArray, *, dim: Dim) -> xr.DataArray:
    """Mean squared error reduced over ``dim``."""
    return _apply_pixel_kernel(_kernels.mse, pred, ref, dim)

nrmse ¶

nrmse(pred: DataArray, ref: DataArray, *, dim: Dim) -> DataArray

Normalized RMSE: 1 - RMSE / sqrt(<ref^2>).

Returns a score in (-inf, 1] where 1 means a perfect match and 0 means the prediction is as wrong as a zero prediction.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/pixel.py

def nrmse(pred: xr.DataArray, ref: xr.DataArray, *, dim: Dim) -> xr.DataArray:
    """Normalized RMSE: ``1 - RMSE / sqrt(<ref^2>)``.

    Returns a score in ``(-inf, 1]`` where 1 means a perfect match and 0
    means the prediction is as wrong as a zero prediction.
    """
    return _apply_pixel_kernel(_kernels.nrmse, pred, ref, dim)

nrmse_score ¶

nrmse_score(pred: DataArray, ref: DataArray, *, dim: Dim) -> DataArray

Mercator/OceanBench-flavour NRMSE skill: 1 - RMSE / std(ref).

Differs from :func:nrmse in the denominator: this uses the anomaly magnitude std(ref) while :func:nrmse uses the raw signal magnitude sqrt(<ref^2>). For zero-mean references the two agree.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/pixel.py

def nrmse_score(pred: xr.DataArray, ref: xr.DataArray, *, dim: Dim) -> xr.DataArray:
    """Mercator/OceanBench-flavour NRMSE skill: ``1 - RMSE / std(ref)``.

    Differs from :func:`nrmse` in the denominator: this uses the
    anomaly magnitude ``std(ref)`` while :func:`nrmse` uses the raw
    signal magnitude ``sqrt(<ref^2>)``. For zero-mean references the
    two agree.
    """
    return _apply_pixel_kernel(_kernels.nrmse_score, pred, ref, dim)

r2_score ¶

r2_score(pred: DataArray, ref: DataArray, *, dim: Dim) -> DataArray

Coefficient of determination: 1 - SS_res / SS_tot.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/pixel.py

def r2_score(pred: xr.DataArray, ref: xr.DataArray, *, dim: Dim) -> xr.DataArray:
    """Coefficient of determination: ``1 - SS_res / SS_tot``."""
    return _apply_pixel_kernel(_kernels.r2_score, pred, ref, dim)

rmse ¶

rmse(pred: DataArray, ref: DataArray, *, dim: Dim) -> DataArray

Root mean squared error reduced over dim.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/pixel.py

def rmse(pred: xr.DataArray, ref: xr.DataArray, *, dim: Dim) -> xr.DataArray:
    """Root mean squared error reduced over ``dim``."""
    return _apply_pixel_kernel(_kernels.rmse, pred, ref, dim)

Spectral — power-spectrum skill¶

Power-spectral-density scores and the resolved-scale diagnostic (the wavelength at which a prediction's PSD score crosses a threshold).

PSDScore ¶

PSDScore(variable: str, psd_dims: Sequence[str], avg_dims: Sequence[str] | None = None, isotropic: bool = False, **kwargs: Any)

Bases: Operator

Two-input PSD score operator.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/spectral.py

def __init__(
    self,
    variable: str,
    psd_dims: Sequence[str],
    avg_dims: Sequence[str] | None = None,
    isotropic: bool = False,
    **kwargs: Any,
):
    self.variable = variable
    self.psd_dims = list(psd_dims)
    self.avg_dims = None if avg_dims is None else list(avg_dims)
    self.isotropic = isotropic
    self.kwargs = dict(kwargs)

SegmentedPSDScore ¶

SegmentedPSDScore(*, var_ref: str, var_pred: str, npt: int, dim: str = 'num_lines', overlap: float = 0.5, max_gap: timedelta64 | None = _DEFAULT_MAX_GAP, spacing_km: float | None = None, lon: str = 'lon', lat: str = 'lat', time: str = 'time', window: str | tuple[str, float] | ArrayLike = 'hann', scaling: str = 'density')

Bases: Operator

Single-input operator for along-track segmented PSD score.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/segmented_psd.py

def __init__(
    self,
    *,
    var_ref: str,
    var_pred: str,
    npt: int,
    dim: str = "num_lines",
    overlap: float = 0.5,
    max_gap: np.timedelta64 | None = _DEFAULT_MAX_GAP,
    spacing_km: float | None = None,
    lon: str = "lon",
    lat: str = "lat",
    time: str = "time",
    window: str | tuple[str, float] | ArrayLike = "hann",
    scaling: str = "density",
):
    self.var_ref = var_ref
    self.var_pred = var_pred
    self.npt = npt
    self.dim = dim
    self.overlap = overlap
    self.max_gap = max_gap
    self.spacing_km = spacing_km
    self.lon = lon
    self.lat = lat
    self.time = time
    self.window = window
    self.scaling = scaling

WaveletPSDScore ¶

WaveletPSDScore(variable: str, scales: Sequence[float] | DataArray, *, dim: tuple[str, str] = ('y', 'x'), x0: float = 50000.0, ntheta: int = 16, k0: float = 1.0, isotropic: bool = True)

Bases: Operator

Two-input localized wavelet PSD score operator.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/spectral.py

def __init__(
    self,
    variable: str,
    scales: Sequence[float] | xr.DataArray,
    *,
    dim: tuple[str, str] = ("y", "x"),
    x0: float = 50e3,
    ntheta: int = 16,
    k0: float = 1.0,
    isotropic: bool = True,
) -> None:
    self.variable = variable
    self.scales = scales
    self.dim = tuple(dim)
    self.x0 = float(x0)
    self.ntheta = int(ntheta)
    self.k0 = float(k0)
    self.isotropic = bool(isotropic)

Functions¶

spectral ¶

Spectral evaluation metrics — public re-export.

Layer-0 functions: :func:psd_error, :func:psd_score, :func:psd_score_spacetime, :func:along_track_psd_score, :func:psd_score_by_region, :func:resolved_scale, :func:resolved_scale_2d, :func:find_intercept_1D, :func:find_intercept_2D, :func:evaluate_by_frequency_band, :func:band_limited_rmse, :func:wavelet_psd_score, :func:wavelet_resolved_scale_map.

Layer-1 operators: :class:PSDScore, :class:FrequencyBandSkill, :class:BandLimitedRMSE, :class:SegmentedPSDScore.

Implementation lives in :mod:xrops.metrics._src.spectral for the core PSD primitives, with the gap-tolerant segmented along-track helpers (:func:along_track_psd_score, :func:psd_score_by_region, :class:SegmentedPSDScore) defined in :mod:xrops.metrics._src.segmented_psd and the 2-D space-time helpers (:func:psd_score_spacetime) in :mod:xrops.metrics._src.composite.

SegmentedPSDScore ¶

SegmentedPSDScore(*, var_ref: str, var_pred: str, npt: int, dim: str = 'num_lines', overlap: float = 0.5, max_gap: timedelta64 | None = _DEFAULT_MAX_GAP, spacing_km: float | None = None, lon: str = 'lon', lat: str = 'lat', time: str = 'time', window: str | tuple[str, float] | ArrayLike = 'hann', scaling: str = 'density')

Bases: Operator

Single-input operator for along-track segmented PSD score.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/segmented_psd.py

def __init__(
    self,
    *,
    var_ref: str,
    var_pred: str,
    npt: int,
    dim: str = "num_lines",
    overlap: float = 0.5,
    max_gap: np.timedelta64 | None = _DEFAULT_MAX_GAP,
    spacing_km: float | None = None,
    lon: str = "lon",
    lat: str = "lat",
    time: str = "time",
    window: str | tuple[str, float] | ArrayLike = "hann",
    scaling: str = "density",
):
    self.var_ref = var_ref
    self.var_pred = var_pred
    self.npt = npt
    self.dim = dim
    self.overlap = overlap
    self.max_gap = max_gap
    self.spacing_km = spacing_km
    self.lon = lon
    self.lat = lat
    self.time = time
    self.window = window
    self.scaling = scaling

BandLimitedRMSE ¶

BandLimitedRMSE(variable: str, dims: Sequence[str], bands: Mapping[str, tuple[float, float]], *, coord_spacing: Mapping[str, float] | None = None)

Bases: FrequencyBandSkill

Fixed-metric convenience subclass: FrequencyBandSkill with RMSE.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/spectral.py

def __init__(
    self,
    variable: str,
    dims: Sequence[str],
    bands: Mapping[str, tuple[float, float]],
    *,
    coord_spacing: Mapping[str, float] | None = None,
) -> None:
    from xrops.metrics._src.pixel import RMSE

    super().__init__(
        variable,
        dims,
        bands,
        metric=RMSE(variable, dims=tuple(dims)),
        coord_spacing=coord_spacing,
    )

FrequencyBandSkill ¶

FrequencyBandSkill(variable: str, dims: Sequence[str], bands: Mapping[str, tuple[float, float]], *, metric: Operator | None = None, coord_spacing: Mapping[str, float] | None = None)

Bases: Operator

Apply an inner metric to band-passed (pred, ref) pairs.

Parameters:

Name	Type	Description	Default
`variable`	`str`	Variable to score.	required
`dims`	`Sequence[str]`	Dimensions to FFT along (passed through to inner metric as well when `metric` is omitted).	required
`bands`	`Mapping[str, tuple[float, float]]`	`{name: (low, high)}` in physical coord units.	required
`metric`	`Operator \| None`	Inner :class:`Operator` (default `RMSE(variable, dims=dims)`).	`None`
`coord_spacing`	`Mapping[str, float] \| None`	Optional per-dim spacing override (in the same physical units as `bands`).	`None`

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/spectral.py

def __init__(
    self,
    variable: str,
    dims: Sequence[str],
    bands: Mapping[str, tuple[float, float]],
    *,
    metric: Operator | None = None,
    coord_spacing: Mapping[str, float] | None = None,
) -> None:
    if metric is None:
        from xrops.metrics._src.pixel import RMSE

        metric = RMSE(variable, dims=tuple(dims))
    elif not isinstance(metric, Operator):
        raise TypeError(
            "metric must be an Operator instance (e.g. RMSE(...)) so its "
            "configuration is introspectable."
        )
    _validate_bands(bands)
    self.variable = variable
    self.dims = tuple(dims)
    self.bands = {k: (float(v[0]), float(v[1])) for k, v in bands.items()}
    self.metric = metric
    self.coord_spacing = (
        None
        if coord_spacing is None
        else {k: float(v) for k, v in coord_spacing.items()}
    )

PSDScore ¶

PSDScore(variable: str, psd_dims: Sequence[str], avg_dims: Sequence[str] | None = None, isotropic: bool = False, **kwargs: Any)

Bases: Operator

Two-input PSD score operator.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/spectral.py

def __init__(
    self,
    variable: str,
    psd_dims: Sequence[str],
    avg_dims: Sequence[str] | None = None,
    isotropic: bool = False,
    **kwargs: Any,
):
    self.variable = variable
    self.psd_dims = list(psd_dims)
    self.avg_dims = None if avg_dims is None else list(avg_dims)
    self.isotropic = isotropic
    self.kwargs = dict(kwargs)

WaveletPSDScore ¶

WaveletPSDScore(variable: str, scales: Sequence[float] | DataArray, *, dim: tuple[str, str] = ('y', 'x'), x0: float = 50000.0, ntheta: int = 16, k0: float = 1.0, isotropic: bool = True)

Bases: Operator

Two-input localized wavelet PSD score operator.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/spectral.py

def __init__(
    self,
    variable: str,
    scales: Sequence[float] | xr.DataArray,
    *,
    dim: tuple[str, str] = ("y", "x"),
    x0: float = 50e3,
    ntheta: int = 16,
    k0: float = 1.0,
    isotropic: bool = True,
) -> None:
    self.variable = variable
    self.scales = scales
    self.dim = tuple(dim)
    self.x0 = float(x0)
    self.ntheta = int(ntheta)
    self.k0 = float(k0)
    self.isotropic = bool(isotropic)

psd_score_spacetime ¶

psd_score_spacetime(pred: DataArray, ref: DataArray, *, space_dim: str = 'lon', time_dim: str = 'time', avg_dims: Sequence[str] | None = ('lat',), level: float = 0.5, **xrft_kwargs: Any) -> tuple[Dataset, dict[str, float]]

Compute a 2-D space-time PSD score and resolved-scale summary.

Parameters:

Name	Type	Description	Default
`pred`	`DataArray`	Prediction DataArray.	required
`ref`	`DataArray`	Reference DataArray (same shape and coords as `pred`).	required
`space_dim`	`str`	Spatial dimension used in the PSD.	`'lon'`
`time_dim`	`str`	Temporal dimension used in the PSD.	`'time'`
`avg_dims`	`Sequence[str] \| None`	Optional dimensions averaged out after the PSD.	`('lat',)`
`level`	`float`	Threshold used by :func:`resolved_scale_2d`.	`0.5`
`**xrft_kwargs`	`Any`	Extra keyword arguments forwarded to :func:`psd_score`.	`{}`

Returns:

Type	Description
`Dataset`	Tuple `(score, summary)` where `score` is a Dataset with a
`dict[str, float]`	`"score"` variable on `(freq_<space_dim>, freq_<time_dim>)`
`tuple[Dataset, dict[str, float]]`	and `summary` contains the min/max resolved wavelengths.

Examples:

>>> score, summary = psd_score_spacetime(pred_da, ref_da)

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/composite.py

def psd_score_spacetime(
    pred: xr.DataArray,
    ref: xr.DataArray,
    *,
    space_dim: str = "lon",
    time_dim: str = "time",
    avg_dims: Sequence[str] | None = ("lat",),
    level: float = 0.5,
    **xrft_kwargs: Any,
) -> tuple[xr.Dataset, dict[str, float]]:
    """Compute a 2-D space-time PSD score and resolved-scale summary.

    Args:
        pred: Prediction DataArray.
        ref: Reference DataArray (same shape and coords as ``pred``).
        space_dim: Spatial dimension used in the PSD.
        time_dim: Temporal dimension used in the PSD.
        avg_dims: Optional dimensions averaged out after the PSD.
        level: Threshold used by :func:`resolved_scale_2d`.
        **xrft_kwargs: Extra keyword arguments forwarded to
            :func:`psd_score`.

    Returns:
        Tuple ``(score, summary)`` where ``score`` is a Dataset with a
        ``"score"`` variable on ``(freq_<space_dim>, freq_<time_dim>)``
        and ``summary`` contains the min/max resolved wavelengths.

    Examples:
        ```pycon
        >>> score, summary = psd_score_spacetime(pred_da, ref_da)
        ```
    """
    if xrft_kwargs.get("isotropic"):
        # Isotropic mode collapses the two spatial axes into a single
        # radial frequency, so the freq_<space_dim>/freq_<time_dim>
        # plane this function indexes into below would not exist.
        raise ValueError(
            "psd_score_spacetime requires the 2-D (space, time) PSD; "
            "isotropic=True is not supported here. Use psd_score "
            "directly for isotropic spectra."
        )
    freq_space_dim = f"freq_{space_dim}"
    freq_time_dim = f"freq_{time_dim}"
    score = psd_score(
        pred,
        ref,
        psd_dims=(space_dim, time_dim),
        avg_dims=avg_dims,
        **xrft_kwargs,
    )
    positive = (score[freq_space_dim] > 0.0) & (score[freq_time_dim] > 0.0)
    score = score.where(positive, drop=True).transpose(freq_space_dim, freq_time_dim)
    summary = resolved_scale_2d(
        score,
        level=level,
        space_dim=freq_space_dim,
        time_dim=freq_time_dim,
    )
    return score, summary

along_track_psd_score ¶

along_track_psd_score(pred: DataArray, ref: DataArray, *, lon: DataArray, lat: DataArray, time: DataArray | None = None, dim: str = 'num_lines', npt: int, overlap: float = 0.5, max_gap: timedelta64 | None = _DEFAULT_MAX_GAP, spacing_km: float | None = None, window: str | tuple[str, float] | ArrayLike = 'hann', scaling: str = 'density') -> Dataset

Compute per-window along-track PSD score from prediction/reference tracks.

Parameters:

Name	Type	Description	Default
`pred`	`DataArray`	Prediction DataArray (1-D along `dim`).	required
`ref`	`DataArray`	Reference DataArray (1-D along `dim`).	required
`lon`	`DataArray`	Longitude DataArray aligned with `pred`/`ref` on `dim`.	required
`lat`	`DataArray`	Latitude DataArray aligned with `pred`/`ref` on `dim`.	required
`time`	`DataArray \| None`	Time DataArray; required when `max_gap` is not `None`.	`None`
`dim`	`str`	Track dimension name (default `"num_lines"`).	`'num_lines'`
`npt`	`int`	Window length in samples.	required
`overlap`	`float`	Welch overlap fraction.	`0.5`
`max_gap`	`timedelta64 \| None`	Maximum allowed time gap before splitting a segment; pass `None` to disable gap detection.	`_DEFAULT_MAX_GAP`
`spacing_km`	`float \| None`	Override along-track spacing; otherwise inferred from `lon`/`lat` via haversine.	`None`
`window`	`str \| tuple[str, float] \| ArrayLike`	Welch window spec.	`'hann'`
`scaling`	`str`	Welch `scaling` (`"density"` or `"spectrum"`).	`'density'`

Returns:

Type	Description
`Dataset`	Per-segment Dataset with `psd_ref`, `psd_pred`, `psd_err`,
`Dataset`	`psd_score` and `coherence` variables.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/segmented_psd.py

def along_track_psd_score(
    pred: xr.DataArray,
    ref: xr.DataArray,
    *,
    lon: xr.DataArray,
    lat: xr.DataArray,
    time: xr.DataArray | None = None,
    dim: str = "num_lines",
    npt: int,
    overlap: float = 0.5,
    max_gap: np.timedelta64 | None = _DEFAULT_MAX_GAP,
    spacing_km: float | None = None,
    window: str | tuple[str, float] | ArrayLike = "hann",
    scaling: str = "density",
) -> xr.Dataset:
    """Compute per-window along-track PSD score from prediction/reference tracks.

    Args:
        pred: Prediction DataArray (1-D along ``dim``).
        ref: Reference DataArray (1-D along ``dim``).
        lon: Longitude DataArray aligned with ``pred``/``ref`` on ``dim``.
        lat: Latitude DataArray aligned with ``pred``/``ref`` on ``dim``.
        time: Time DataArray; required when ``max_gap`` is not ``None``.
        dim: Track dimension name (default ``"num_lines"``).
        npt: Window length in samples.
        overlap: Welch overlap fraction.
        max_gap: Maximum allowed time gap before splitting a segment; pass
            ``None`` to disable gap detection.
        spacing_km: Override along-track spacing; otherwise inferred from
            ``lon``/``lat`` via haversine.
        window: Welch window spec.
        scaling: Welch ``scaling`` (``"density"`` or ``"spectrum"``).

    Returns:
        Per-segment Dataset with ``psd_ref``, ``psd_pred``, ``psd_err``,
        ``psd_score`` and ``coherence`` variables.
    """
    if dim not in pred.dims or dim not in ref.dims:
        raise ValueError(f"{dim!r} is not a dimension of pred/ref.")

    ref_values = np.asarray(ref.transpose(dim).values, dtype=float)
    pred_values = np.asarray(pred.transpose(dim).values, dtype=float)
    if ref_values.ndim != 1 or pred_values.ndim != 1:
        raise ValueError("pred and ref must be 1-D along the track dimension.")
    if ref_values.shape != pred_values.shape:
        raise ValueError("pred and ref must have the same shape.")

    lon_values = _coord_values_from_da(lon, name="lon", dim=dim)
    lat_values = _coord_values_from_da(lat, name="lat", dim=dim)
    spacing = (
        _median_dx_km(lon_values, lat_values)
        if spacing_km is None
        else float(spacing_km)
    )
    fs = 1.0 / spacing
    if max_gap is not None:
        if time is None:
            raise ValueError("time DataArray is required when max_gap is not None.")
        gaps = _gap_indices_from_da(time, dim=dim, max_gap=max_gap)
    else:
        gaps = np.empty(0, dtype=int)
    bounds = _segment_bounds(
        ref_values.size, npt=npt, overlap=overlap, gap_indices=gaps
    )

    ref_segments = _segment_stack(ref_values, bounds)
    pred_segments = _segment_stack(pred_values, bounds)
    lon_segments = _segment_stack(lon_values, bounds)
    lat_segments = _segment_stack(lat_values, bounds)
    if ref_segments.size == 0:
        freqs, psd_ref, psd_pred, psd_err, coherence = _empty_spectra(
            npt=npt, fs=fs, window=window, scaling=scaling
        )
        segment_lons = np.empty(0, dtype=float)
        segment_lats = np.empty(0, dtype=float)
    else:
        finite = np.all(
            np.isfinite(ref_segments)
            & np.isfinite(pred_segments)
            & np.isfinite(lon_segments)
            & np.isfinite(lat_segments),
            axis=1,
        )
        bounds = [bounds[i] for i in np.flatnonzero(finite)]
        ref_segments = ref_segments[finite]
        pred_segments = pred_segments[finite]
        if ref_segments.shape[0] == 0:
            freqs, psd_ref, psd_pred, psd_err, coherence = _empty_spectra(
                npt=npt, fs=fs, window=window, scaling=scaling
            )
            segment_lons = np.empty(0, dtype=float)
            segment_lats = np.empty(0, dtype=float)
        else:
            err_segments = pred_segments - ref_segments
            freqs, psd_ref = signal.welch(
                ref_segments,
                fs=fs,
                window=window,
                nperseg=npt,
                noverlap=0,
                scaling=scaling,
                axis=-1,
            )
            _, psd_pred = signal.welch(
                pred_segments,
                fs=fs,
                window=window,
                nperseg=npt,
                noverlap=0,
                scaling=scaling,
                axis=-1,
            )
            _, psd_err = signal.welch(
                err_segments,
                fs=fs,
                window=window,
                nperseg=npt,
                noverlap=0,
                scaling=scaling,
                axis=-1,
            )
            _, csd = signal.csd(
                ref_segments,
                pred_segments,
                fs=fs,
                window=window,
                nperseg=npt,
                noverlap=0,
                scaling=scaling,
                axis=-1,
            )
            coherence = np.abs(csd) ** 2 / (psd_ref * psd_pred)
            segment_lons = _segment_lons(lon_values, bounds)
            segment_lats = np.asarray(
                [np.nanmedian(lat_values[start:stop]) for start, stop in bounds],
                dtype=float,
            )

    with np.errstate(divide="ignore", invalid="ignore"):
        score = np.where(
            np.isfinite(psd_ref) & (psd_ref > 0.0),
            1.0 - psd_err / psd_ref,
            np.nan,
        )
        wavelength = 1.0 / freqs

    coords = {
        "segment": np.arange(segment_lons.size),
        "wavenumber": ("wavenumber", freqs, {"units": "cycles/km"}),
        "wavelength": ("wavenumber", wavelength, {"units": "km"}),
        "segment_lon": ("segment", segment_lons, {"units": "degrees_east"}),
        "segment_lat": ("segment", segment_lats, {"units": "degrees_north"}),
    }
    return xr.Dataset(
        data_vars={
            "psd_ref": (("segment", "wavenumber"), psd_ref),
            "psd_pred": (("segment", "wavenumber"), psd_pred),
            "psd_err": (("segment", "wavenumber"), psd_err),
            "psd_score": (("segment", "wavenumber"), score),
            "coherence": (("segment", "wavenumber"), coherence),
        },
        coords=coords,
        attrs={"spacing_km": spacing},
    )

psd_score_by_region ¶

psd_score_by_region(ds_segments: Dataset, *, lat_centers: ArrayLike = _DEFAULT_LAT_CENTERS, lon_centers: ArrayLike = _DEFAULT_LON_CENTERS, delta_lat: float = 10.0, delta_lon: float = 10.0, min_segments: int = 2) -> Dataset

Average per-segment PSDs into overlapping circular-longitude boxes.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/segmented_psd.py

def psd_score_by_region(
    ds_segments: xr.Dataset,
    *,
    lat_centers: ArrayLike = _DEFAULT_LAT_CENTERS,
    lon_centers: ArrayLike = _DEFAULT_LON_CENTERS,
    delta_lat: float = 10.0,
    delta_lon: float = 10.0,
    min_segments: int = 2,
) -> xr.Dataset:
    """Average per-segment PSDs into overlapping circular-longitude boxes."""
    lat_values = np.asarray(ds_segments["segment_lat"].values, dtype=float)
    lon_values = np.mod(
        np.asarray(ds_segments["segment_lon"].values, dtype=float), 360.0
    )
    lat_out = np.asarray(lat_centers, dtype=float)
    # Normalize to [0, 360) so callers can pass either [-180, 180] or
    # [0, 360] convention; the circular-distance computation downstream
    # assumes a common modulus.
    lon_out = np.mod(np.asarray(lon_centers, dtype=float), 360.0)
    wavenumber = np.asarray(ds_segments["wavenumber"].values, dtype=float)

    data_names = [
        name
        for name in ("psd_ref", "psd_pred", "psd_err", "coherence")
        if name in ds_segments
    ]
    out = {
        name: np.full(
            (lat_out.size, lon_out.size, wavenumber.size), np.nan, dtype=float
        )
        for name in data_names
    }
    counts = np.zeros((lat_out.size, lon_out.size), dtype=int)

    for i, lat_center in enumerate(lat_out):
        lat_mask = np.abs(lat_values - lat_center) <= delta_lat / 2.0
        for j, lon_center in enumerate(lon_out):
            # Wrap to [-180, 180] to get shortest circular longitude distance.
            lon_delta = np.abs((lon_values - lon_center + 180.0) % 360.0 - 180.0)
            mask = lat_mask & (lon_delta <= delta_lon / 2.0)
            counts[i, j] = int(np.count_nonzero(mask))
            if counts[i, j] < min_segments:
                continue
            for name in data_names:
                out[name][i, j] = np.nanmean(ds_segments[name].values[mask], axis=0)

    coords = {
        "lat": ("lat", lat_out, {"units": "degrees_north"}),
        "lon": ("lon", lon_out, {"units": "degrees_east"}),
        "wavenumber": ("wavenumber", wavenumber, ds_segments["wavenumber"].attrs),
    }
    if "wavelength" in ds_segments:
        coords["wavelength"] = (
            "wavenumber",
            ds_segments["wavelength"].values,
            ds_segments["wavelength"].attrs,
        )

    data_vars: dict[str, Any] = {
        name: (("lat", "lon", "wavenumber"), values) for name, values in out.items()
    }
    # psd_score is only well-defined when both reference and error PSDs
    # are present; otherwise emit NaN rather than raising on a missing key.
    if "psd_ref" in out and "psd_err" in out:
        with np.errstate(divide="ignore", invalid="ignore"):
            score = np.where(
                np.isfinite(out["psd_ref"]) & (out["psd_ref"] > 0.0),
                1.0 - out["psd_err"] / out["psd_ref"],
                np.nan,
            )
        data_vars["psd_score"] = (("lat", "lon", "wavenumber"), score)
    data_vars["n_segments"] = (("lat", "lon"), counts)
    return xr.Dataset(data_vars=data_vars, coords=coords)

band_limited_rmse ¶

band_limited_rmse(prediction: Dataset, reference: Dataset, *, variable: str, bands: Mapping[str, tuple[float, float]], dims: Sequence[str], coord_spacing: Mapping[str, float] | None = None) -> Dataset

RMSE of prediction - reference filtered to each band.

Convenience wrapper around :func:evaluate_by_frequency_band with metric=RMSE(variable, dims=dims).

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/spectral.py

def band_limited_rmse(
    prediction: xr.Dataset,
    reference: xr.Dataset,
    *,
    variable: str,
    bands: Mapping[str, tuple[float, float]],
    dims: Sequence[str],
    coord_spacing: Mapping[str, float] | None = None,
) -> xr.Dataset:
    """RMSE of ``prediction - reference`` filtered to each band.

    Convenience wrapper around :func:`evaluate_by_frequency_band` with
    ``metric=RMSE(variable, dims=dims)``.
    """
    from xrops.metrics._src.pixel import RMSE

    return evaluate_by_frequency_band(
        prediction,
        reference,
        variable=variable,
        bands=bands,
        dims=dims,
        metric=RMSE(variable, dims=tuple(dims)),
        coord_spacing=coord_spacing,
    )

evaluate_by_frequency_band ¶

evaluate_by_frequency_band(prediction: Dataset, reference: Dataset, *, variable: str, bands: Mapping[str, tuple[float, float]], dims: Sequence[str], metric: Operator | None = None, coord_spacing: Mapping[str, float] | None = None) -> Dataset

Filter prediction and reference to each band and apply metric.

Parameters:

Name	Type	Description	Default
`prediction`	`Dataset`	Prediction dataset.	required
`reference`	`Dataset`	Reference dataset.	required
`variable`	`str`	Variable to score; must exist in both inputs.	required
`bands`	`Mapping[str, tuple[float, float]]`	`{name: (low, high)}` in physical coord units of `dims`. `low` is included, `high` is excluded; both must be `>= 0`.	required
`dims`	`Sequence[str]`	Dimensions to FFT along. Coords must carry `units` metadata (or per-dim `coord_spacing` overrides must be supplied).	required
`metric`	`Operator \| None`	Inner :class:`Operator` evaluated on the band-passed `(pred, ref)` pair. Defaults to `RMSE(variable, dims=dims)`.	`None`
`coord_spacing`	`Mapping[str, float] \| None`	Optional per-dim spacing overrides (in the same physical units as `bands`); bypasses unit auto-detection.	`None`

Returns:

Type	Description
`Dataset`	Dataset stacked along a `"band"` dim with auxiliary
`Dataset`	`"band_low"` / `"band_high"` coords. Bands whose `low`
`Dataset`	exceeds the Nyquist magnitude produce NaN entries plus a
`Dataset`	class:`UserWarning` (no exception).

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/spectral.py

def evaluate_by_frequency_band(
    prediction: xr.Dataset,
    reference: xr.Dataset,
    *,
    variable: str,
    bands: Mapping[str, tuple[float, float]],
    dims: Sequence[str],
    metric: Operator | None = None,
    coord_spacing: Mapping[str, float] | None = None,
) -> xr.Dataset:
    """Filter ``prediction`` and ``reference`` to each band and apply ``metric``.

    Args:
        prediction: Prediction dataset.
        reference: Reference dataset.
        variable: Variable to score; must exist in both inputs.
        bands: ``{name: (low, high)}`` in physical coord units of
            ``dims``. ``low`` is included, ``high`` is excluded; both
            must be ``>= 0``.
        dims: Dimensions to FFT along. Coords must carry ``units``
            metadata (or per-dim ``coord_spacing`` overrides must be
            supplied).
        metric: Inner :class:`Operator` evaluated on the band-passed
            ``(pred, ref)`` pair. Defaults to ``RMSE(variable, dims=dims)``.
        coord_spacing: Optional per-dim spacing overrides (in the same
            physical units as ``bands``); bypasses unit auto-detection.

    Returns:
        Dataset stacked along a ``"band"`` dim with auxiliary
        ``"band_low"`` / ``"band_high"`` coords. Bands whose ``low``
        exceeds the Nyquist magnitude produce NaN entries plus a
        :class:`UserWarning` (no exception).
    """
    if variable not in prediction.data_vars:
        raise KeyError(f"prediction missing variable {variable!r}")
    if variable not in reference.data_vars:
        raise KeyError(f"reference missing variable {variable!r}")
    dim_list = list(dims)
    for d in dim_list:
        if d not in prediction.dims:
            raise ValueError(f"prediction is missing dim {d!r}")
        if d not in reference.dims:
            raise ValueError(f"reference is missing dim {d!r}")
    # Bands are physical-frequency, so prediction and reference must
    # share the dim coordinates — otherwise a single ``spacing`` derived
    # from one of them would band-pass the other with the wrong cutoff.
    for d in dim_list:
        if d in prediction.coords and d in reference.coords:
            p = np.asarray(prediction[d].values)
            r = np.asarray(reference[d].values)
            if p.shape != r.shape or not np.allclose(p, r):
                raise ValueError(
                    f"prediction and reference disagree on coord {d!r}; "
                    "regrid one onto the other (e.g. with "
                    "`xrops.interpolate.RegridLike`) before scoring."
                )
    _validate_bands(bands)

    if metric is None:
        from xrops.metrics._src.pixel import RMSE

        metric = RMSE(variable, dims=tuple(dim_list))
    elif not isinstance(metric, Operator):
        raise TypeError(
            "metric must be an Operator instance (e.g. RMSE(...)) so its "
            "configuration is introspectable."
        )

    spacing = _physical_spacing(prediction, dim_list, coord_spacing=coord_spacing)
    nyquist_mag = float(np.sqrt(sum((0.5 / spacing[d]) ** 2 for d in dim_list)))

    pred_da = prediction[variable]
    ref_da = reference[variable]

    pieces: list[xr.DataArray | xr.Dataset] = []
    band_names = list(bands.keys())
    for name in band_names:
        lo, hi = float(bands[name][0]), float(bands[name][1])
        if lo > nyquist_mag:
            warnings.warn(
                f"band {name!r} low={lo} exceeds Nyquist magnitude "
                f"{nyquist_mag:.4g} — emitting NaN.",
                UserWarning,
                stacklevel=2,
            )
            template = metric(prediction, reference)
            pieces.append(xr.full_like(template, np.nan, dtype=float))
            continue
        pred_band = _bandpass(pred_da, dim_list, spacing, lo, hi)
        ref_band = _bandpass(ref_da, dim_list, spacing, lo, hi)
        pred_ds = prediction.assign({variable: pred_band})
        ref_ds = reference.assign({variable: ref_band})
        pieces.append(metric(pred_ds, ref_ds))

    band_index = xr.DataArray(band_names, dims=("band",), name="band")
    out = xr.concat(pieces, dim=band_index)
    band_low = xr.DataArray(
        [float(bands[k][0]) for k in band_names], dims=("band",), name="band_low"
    )
    band_high = xr.DataArray(
        [float(bands[k][1]) for k in band_names], dims=("band",), name="band_high"
    )
    out = out.assign_coords(band_low=band_low, band_high=band_high)
    if isinstance(out, xr.DataArray):
        return out.to_dataset(name=out.name or "score")
    return out

find_intercept_1D ¶

find_intercept_1D(x: Float[ndarray, 'n'], y: Float[ndarray, 'n'], level: float = 0.5, kind: str = 'slinear', **kwargs: Any) -> float

Invert a 1-D monotone-ish curve at y = level and return x.

Uses :class:scipy.interpolate.interp1d on (y, x). Duplicate y values (common for plateaued PSD scores) are collapsed to their first occurrence in the sorted order before interpolating, because interp1d requires a strictly monotone x-axis. Extrapolates silently when level falls outside the range of the deduplicated y.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/spectral.py

def find_intercept_1D(
    x: Float[np.ndarray, "n"],
    y: Float[np.ndarray, "n"],
    level: float = 0.5,
    kind: str = "slinear",
    **kwargs: Any,
) -> float:
    """Invert a 1-D monotone-ish curve at ``y = level`` and return ``x``.

    Uses :class:`scipy.interpolate.interp1d` on ``(y, x)``. Duplicate
    ``y`` values (common for plateaued PSD scores) are collapsed to
    their first occurrence in the sorted order before interpolating,
    because ``interp1d`` requires a strictly monotone x-axis.
    Extrapolates silently when ``level`` falls outside the range of the
    deduplicated ``y``.
    """
    y_arr = np.asarray(y, dtype=float)
    x_arr = np.asarray(x, dtype=float)
    order = np.argsort(y_arr)
    y_sorted = y_arr[order]
    x_sorted = x_arr[order]

    _, first_idx = np.unique(y_sorted, return_index=True)
    first_idx.sort()
    y_unique = y_sorted[first_idx]
    x_unique = x_sorted[first_idx]

    if y_unique.size < 2:
        return float(x_unique.item()) if x_unique.size else float("nan")

    f = interp1d(
        y_unique,
        x_unique,
        fill_value=kwargs.pop("fill_value", "extrapolate"),
        kind=kind,
        **kwargs,
    )
    try:
        return float(np.asarray(f(level)).item())
    except ValueError:
        warnings.warn(
            f"level={level} outside range of y — returning edge value.",
            stacklevel=2,
        )
        y_min, y_max = float(y_unique.min()), float(y_unique.max())
        edge = y_min if level < y_min else y_max
        return float(np.asarray(f(edge)).item())

find_intercept_2D ¶

find_intercept_2D(score: DataArray, level: float = 0.5, space_dim: str = 'freq_lon', time_dim: str = 'freq_time') -> list[DataArray]

Extract the threshold contour of a 2-D (space, time) score field.

Counterpart to :func:find_intercept_1D for the 2-D space-time score: returns the boundary where score == level as data so it can be quoted in a paper or compared across methods quantitatively rather than only overlaid on a figure. Wraps :func:skimage.measure.find_contours. Disconnected boundary pieces are returned as separate segments. Coordinate values along each segment are interpolated linearly from the input grid. Args: score: 2-D :class:xr.DataArray indexed by space_dim and time_dim. level: Score threshold to extract. Default 0.5. space_dim: Name of the spatial-frequency dim. Default "freq_lon". time_dim: Name of the temporal-frequency dim. Default "freq_time". Returns: List of :class:xr.DataArray, one per contour segment, each with dims ("point", "axis") where axis has coordinates (space_dim, time_dim) and contains the segment's coordinate polyline.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/spectral.py

def find_intercept_2D(
    score: xr.DataArray,
    level: float = 0.5,
    space_dim: str = "freq_lon",
    time_dim: str = "freq_time",
) -> list[xr.DataArray]:
    """Extract the threshold contour of a 2-D ``(space, time)`` score field.

    Counterpart to :func:`find_intercept_1D` for the 2-D space-time
    score: returns the boundary where ``score == level`` as data so it
    can be quoted in a paper or compared across methods quantitatively
    rather than only overlaid on a figure.
    Wraps :func:`skimage.measure.find_contours`. Disconnected boundary
    pieces are returned as separate segments. Coordinate values along
    each segment are interpolated linearly from the input grid.
    Args:
        score: 2-D :class:`xr.DataArray` indexed by ``space_dim`` and
            ``time_dim``.
        level: Score threshold to extract. Default ``0.5``.
        space_dim: Name of the spatial-frequency dim. Default
            ``"freq_lon"``.
        time_dim: Name of the temporal-frequency dim. Default
            ``"freq_time"``.
    Returns:
        List of :class:`xr.DataArray`, one per contour segment, each
        with dims ``("point", "axis")`` where ``axis`` has coordinates
        ``(space_dim, time_dim)`` and contains the segment's coordinate
        polyline.
    """
    measure = _require_optional(
        "skimage.measure",
        extra="image",
        feature="find_intercept_2D",
        package="scikit-image",
    )
    find_contours = measure.find_contours

    if score.ndim != 2:
        raise ValueError(
            f"find_intercept_2D expects a 2-D DataArray; got dims {score.dims}."
        )
    if space_dim not in score.dims or time_dim not in score.dims:
        raise ValueError(
            f"score must have dims {space_dim!r} and {time_dim!r}; got {score.dims}."
        )

    # find_contours operates in (row, col) index space; transpose so
    # rows == time_dim (y-axis) and cols == space_dim (x-axis), which
    # matches the conventional space-time score plot orientation.
    s = score.transpose(time_dim, space_dim)
    arr = np.asarray(s.values, dtype=float)
    space_coord = np.asarray(s[space_dim].values, dtype=float)
    time_coord = np.asarray(s[time_dim].values, dtype=float)

    contours = find_contours(arr, level=level)
    segments: list[xr.DataArray] = []
    for c in contours:
        # c[:, 0] is fractional row (time), c[:, 1] is fractional col (space)
        rows = c[:, 0]
        cols = c[:, 1]
        sx = np.interp(cols, np.arange(space_coord.size), space_coord)
        ty = np.interp(rows, np.arange(time_coord.size), time_coord)
        seg = xr.DataArray(
            np.stack([sx, ty], axis=-1),
            dims=("point", "axis"),
            coords={"axis": [space_dim, time_dim]},
        )
        segments.append(seg)
    return segments

psd_error ¶

psd_error(pred: DataArray, ref: DataArray, *, psd_dims: Sequence[str], avg_dims: Sequence[str] | None = None, isotropic: bool = False, **kwargs: Any) -> Dataset

PSD of the prediction error pred - ref.

Parameters:

Name	Type	Description	Default
`pred`	`DataArray`	Prediction DataArray.	required
`ref`	`DataArray`	Reference DataArray.	required
`psd_dims`	`Sequence[str]`	Dimensions over which to take the PSD.	required
`avg_dims`	`Sequence[str] \| None`	Optional dims to conditionally average out after the PSD (e.g. average over `lat` after computing the lon/time spectrum).	`None`
`isotropic`	`bool`	If `True`, use the isotropic power spectrum.	`False`
`**kwargs`	`Any`	Forwarded to the underlying PSD function.	`{}`

Returns:

Type	Description
`Dataset`	Dataset with a single `"error"` variable containing the
`Dataset`	PSD of the error.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/spectral.py

def psd_error(
    pred: xr.DataArray,
    ref: xr.DataArray,
    *,
    psd_dims: Sequence[str],
    avg_dims: Sequence[str] | None = None,
    isotropic: bool = False,
    **kwargs: Any,
) -> xr.Dataset:
    """PSD of the prediction error ``pred - ref``.

    Args:
        pred: Prediction DataArray.
        ref: Reference DataArray.
        psd_dims: Dimensions over which to take the PSD.
        avg_dims: Optional dims to conditionally average out after the
            PSD (e.g. average over ``lat`` after computing the lon/time
            spectrum).
        isotropic: If ``True``, use the isotropic power spectrum.
        **kwargs: Forwarded to the underlying PSD function.

    Returns:
        Dataset with a single ``"error"`` variable containing the
        PSD of the error.
    """
    diff = (pred - ref).rename("error")
    err = power_spectrum(diff, dim=list(psd_dims), isotropic=isotropic, **kwargs)
    err_ds = err.rename("error").to_dataset()
    if avg_dims is not None:
        err_ds = drop_negative_frequencies(err_ds, dims=avg_dims, drop=True)
    return err_ds

psd_score ¶

psd_score(pred: DataArray, ref: DataArray, *, psd_dims: Sequence[str], avg_dims: Sequence[str] | None = None, isotropic: bool = False, **kwargs: Any) -> Dataset

Normalized PSD score 1 - PSD(err) / PSD(ref).

A score of 1 means the error has no power at that scale; 0 means the error has as much power as the reference signal.

Parameters:

Name	Type	Description	Default
`pred`	`DataArray`	Prediction DataArray.	required
`ref`	`DataArray`	Reference DataArray.	required
`psd_dims`	`Sequence[str]`	Dimensions over which to take the PSD.	required
`avg_dims`	`Sequence[str] \| None`	Optional conditional-average dims applied after PSD.	`None`
`isotropic`	`bool`	If `True`, use the isotropic power spectrum.	`False`
`**kwargs`	`Any`	Forwarded to the underlying PSD function.	`{}`

Returns:

Type	Description
`Dataset`	Dataset with a single `"score"` variable.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/spectral.py

def psd_score(
    pred: xr.DataArray,
    ref: xr.DataArray,
    *,
    psd_dims: Sequence[str],
    avg_dims: Sequence[str] | None = None,
    isotropic: bool = False,
    **kwargs: Any,
) -> xr.Dataset:
    """Normalized PSD score ``1 - PSD(err) / PSD(ref)``.

    A score of ``1`` means the error has no power at that scale; ``0``
    means the error has as much power as the reference signal.

    Args:
        pred: Prediction DataArray.
        ref: Reference DataArray.
        psd_dims: Dimensions over which to take the PSD.
        avg_dims: Optional conditional-average dims applied after PSD.
        isotropic: If ``True``, use the isotropic power spectrum.
        **kwargs: Forwarded to the underlying PSD function.

    Returns:
        Dataset with a single ``"score"`` variable.
    """
    err = psd_error(
        pred,
        ref,
        psd_dims=psd_dims,
        avg_dims=avg_dims,
        isotropic=isotropic,
        **kwargs,
    )
    ref_psd = power_spectrum(ref, dim=list(psd_dims), isotropic=isotropic, **kwargs)
    ref_name = ref.name if ref.name is not None else "ref"
    ref_ds = ref_psd.rename(ref_name).to_dataset()
    if avg_dims is not None:
        ref_ds = drop_negative_frequencies(ref_ds, dims=avg_dims, drop=True)
    score = 1.0 - err["error"] / ref_ds[ref_name]
    return score.to_dataset(name="score")

resolved_scale ¶

resolved_scale(score: DataArray | Dataset, frequency: str, level: float = 0.5) -> float

Wavelength (1 / frequency) at which a PSD score crosses level.

Parameters:

Name	Type	Description	Default
`score`	`DataArray \| Dataset`	DataArray of PSD-score values along `frequency`, or a Dataset containing a `"score"` variable.	required
`frequency`	`str`	Name of the frequency coordinate.	required
`level`	`float`	Score threshold (default 0.5 — the commonly-used "resolved scale" threshold).	`0.5`

Returns:

Type	Description
`float`	Scalar wavelength at which the score first crosses `level`.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/spectral.py

def resolved_scale(
    score: xr.DataArray | xr.Dataset,
    frequency: str,
    level: float = 0.5,
) -> float:
    """Wavelength (``1 / frequency``) at which a PSD score crosses ``level``.

    Args:
        score: DataArray of PSD-score values along ``frequency``, or a
            Dataset containing a ``"score"`` variable.
        frequency: Name of the frequency coordinate.
        level: Score threshold (default 0.5 — the commonly-used
            "resolved scale" threshold).

    Returns:
        Scalar wavelength at which the score first crosses ``level``.
    """
    score_da = score["score"] if isinstance(score, xr.Dataset) else score
    freqs = np.asarray(score_da[frequency].values)
    vals = np.asarray(score_da.values)
    positive = freqs > 0
    freqs = freqs[positive]
    vals = vals[positive]
    wavelengths = 1.0 / freqs
    return find_intercept_1D(x=wavelengths, y=vals, level=level)

resolved_scale_2d ¶

resolved_scale_2d(score: DataArray | Dataset, *, level: float = 0.5, space_dim: str = 'freq_lon', time_dim: str = 'freq_time') -> dict[str, float]

Summarize a 2-D contour crossing as min/max space-time wavelengths.

Parameters:

Name	Type	Description	Default
`score`	`DataArray \| Dataset`	2-D score field, or a Dataset containing `"score"`.	required
`level`	`float`	Contour level to summarize.	`0.5`
`space_dim`	`str`	Spatial-frequency dimension name.	`'freq_lon'`
`time_dim`	`str`	Temporal-frequency dimension name.	`'freq_time'`

Returns:

Type	Description
`dict[str, float]`	Dictionary with `lambda_space_min`, `lambda_time_min`,
`dict[str, float]`	`lambda_space_max`, and `lambda_time_max`. If no contour
`dict[str, float]`	exists at `level`, all entries are `NaN`.

Examples:

>>> resolved_scale_2d(score, level=0.5)

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/spectral.py

def resolved_scale_2d(
    score: xr.DataArray | xr.Dataset,
    *,
    level: float = 0.5,
    space_dim: str = "freq_lon",
    time_dim: str = "freq_time",
) -> dict[str, float]:
    """Summarize a 2-D contour crossing as min/max space-time wavelengths.

    Args:
        score: 2-D score field, or a Dataset containing ``"score"``.
        level: Contour level to summarize.
        space_dim: Spatial-frequency dimension name.
        time_dim: Temporal-frequency dimension name.

    Returns:
        Dictionary with ``lambda_space_min``, ``lambda_time_min``,
        ``lambda_space_max``, and ``lambda_time_max``. If no contour
        exists at ``level``, all entries are ``NaN``.

    Examples:
        ```pycon
        >>> resolved_scale_2d(score, level=0.5)
        ```
    """
    score_da = score["score"] if isinstance(score, xr.Dataset) else score
    segments = find_intercept_2D(
        score_da,
        level=level,
        space_dim=space_dim,
        time_dim=time_dim,
    )
    if not segments:
        return {
            "lambda_space_min": float("nan"),
            "lambda_time_min": float("nan"),
            "lambda_space_max": float("nan"),
            "lambda_time_max": float("nan"),
        }

    space_freq = np.concatenate(
        [
            np.asarray(segment.sel(axis=space_dim).values, dtype=float)
            for segment in segments
        ]
    )
    time_freq = np.concatenate(
        [
            np.asarray(segment.sel(axis=time_dim).values, dtype=float)
            for segment in segments
        ]
    )
    space_freq = space_freq[space_freq > 0.0]
    time_freq = time_freq[time_freq > 0.0]
    return {
        "lambda_space_min": (
            float(1.0 / np.max(space_freq)) if space_freq.size else float("nan")
        ),
        "lambda_time_min": (
            float(1.0 / np.max(time_freq)) if time_freq.size else float("nan")
        ),
        "lambda_space_max": (
            float(1.0 / np.min(space_freq)) if space_freq.size else float("nan")
        ),
        "lambda_time_max": (
            float(1.0 / np.min(time_freq)) if time_freq.size else float("nan")
        ),
    }

wavelet_psd_score ¶

wavelet_psd_score(pred: DataArray, ref: DataArray, scales: DataArray, *, dim: tuple[str, str] = ('y', 'x'), x0: float = 50000.0, ntheta: int = 16, k0: float = 1.0, isotropic: bool = True) -> Dataset

Localized wavelet PSD score 1 - WPSD(err) / WPSD(ref).

Parameters:

Name	Type	Description	Default
`pred`	`DataArray`	Prediction DataArray on a locally Cartesian grid.	required
`ref`	`DataArray`	Reference DataArray on the same grid.	required
`scales`	`DataArray`	Positive Morlet scales. The first dimension is used as the scale axis.	required
`dim`	`tuple[str, str]`	Spatial dimensions as `(y, x)` on a locally Cartesian grid.	`('y', 'x')`
`x0`	`float`	Reference length scale in the same units as the spatial coords.	`50000.0`
`ntheta`	`int`	Number of evenly spaced Morlet angles.	`16`
`k0`	`float`	Dimensionless central Morlet wavenumber.	`1.0`
`isotropic`	`bool`	If `True`, average directional power over angle.	`True`

Returns:

Type	Description
`Dataset`	Dataset containing `score` with untrusted COI samples masked to NaN.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/spectral.py

def wavelet_psd_score(
    pred: xr.DataArray,
    ref: xr.DataArray,
    scales: xr.DataArray,
    *,
    dim: tuple[str, str] = ("y", "x"),
    x0: float = 50e3,
    ntheta: int = 16,
    k0: float = 1.0,
    isotropic: bool = True,
) -> xr.Dataset:
    """Localized wavelet PSD score ``1 - WPSD(err) / WPSD(ref)``.

    Args:
        pred: Prediction DataArray on a locally Cartesian grid.
        ref: Reference DataArray on the same grid.
        scales: Positive Morlet scales. The first dimension is used as the
            scale axis.
        dim: Spatial dimensions as ``(y, x)`` on a locally Cartesian grid.
        x0: Reference length scale in the same units as the spatial coords.
        ntheta: Number of evenly spaced Morlet angles.
        k0: Dimensionless central Morlet wavenumber.
        isotropic: If ``True``, average directional power over angle.

    Returns:
        Dataset containing ``score`` with untrusted COI samples masked to NaN.
    """
    err = (pred - ref).rename("error")
    err_psd = wvlt_power_spectrum(
        err,
        scales,
        dim=dim,
        x0=x0,
        ntheta=ntheta,
        k0=k0,
        isotropic=isotropic,
    )
    ref_psd = wvlt_power_spectrum(
        ref,
        scales,
        dim=dim,
        x0=x0,
        ntheta=ntheta,
        k0=k0,
        isotropic=isotropic,
    )
    trust = err_psd["coi_mask"] & ref_psd["coi_mask"]
    score = (1.0 - err_psd / ref_psd).where(trust)
    return score.rename("score").to_dataset()

wavelet_resolved_scale_map ¶

wavelet_resolved_scale_map(truth: DataArray, pred: DataArray, scales: DataArray, *, dim: tuple[str, str] = ('y', 'x'), x0: float = 50000.0, ntheta: int = 16, k0: float = 1.0, threshold: float = 0.5, wavelength_scale: float = 0.001, wavelength_units: str = 'km') -> DataArray

Return the wavelength where local wavelet PSD skill crosses a threshold.

Parameters:

Name	Type	Description	Default
`truth`	`DataArray`	Reference field on a locally Cartesian grid.	required
`pred`	`DataArray`	Prediction field on the same grid.	required
`scales`	`DataArray`	Positive Morlet scales used to compute the local spectra.	required
`dim`	`tuple[str, str]`	Spatial dimensions as `(y, x)`.	`('y', 'x')`
`x0`	`float`	Reference length scale in the same units as the spatial coords.	`50000.0`
`ntheta`	`int`	Number of evenly spaced Morlet angles.	`16`
`k0`	`float`	Dimensionless central Morlet wavenumber.	`1.0`
`threshold`	`float`	Score level defining the resolved scale.	`0.5`
`wavelength_scale`	`float`	Multiplier applied to `1 / wavenumber` to convert into `wavelength_units`. Defaults to `1e-3` for metres → km (matches the historical behaviour); pass `1.0` if your coords are already in km.	`0.001`
`wavelength_units`	`str`	Unit string written to the output attrs.	`'km'`

Returns:

Type	Description
`DataArray`	Two-dimensional resolved wavelength in `wavelength_units`.
`DataArray`	Locations with fewer than two trusted scale samples, or where
`DataArray`	the threshold isn't bracketed by the score column, are NaN.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/spectral.py

def wavelet_resolved_scale_map(
    truth: xr.DataArray,
    pred: xr.DataArray,
    scales: xr.DataArray,
    *,
    dim: tuple[str, str] = ("y", "x"),
    x0: float = 50e3,
    ntheta: int = 16,
    k0: float = 1.0,
    threshold: float = 0.5,
    wavelength_scale: float = 1e-3,
    wavelength_units: str = "km",
) -> xr.DataArray:
    """Return the wavelength where local wavelet PSD skill crosses a threshold.

    Args:
        truth: Reference field on a locally Cartesian grid.
        pred: Prediction field on the same grid.
        scales: Positive Morlet scales used to compute the local spectra.
        dim: Spatial dimensions as ``(y, x)``.
        x0: Reference length scale in the same units as the spatial coords.
        ntheta: Number of evenly spaced Morlet angles.
        k0: Dimensionless central Morlet wavenumber.
        threshold: Score level defining the resolved scale.
        wavelength_scale: Multiplier applied to ``1 / wavenumber`` to
            convert into ``wavelength_units``. Defaults to ``1e-3`` for
            metres → km (matches the historical behaviour); pass
            ``1.0`` if your coords are already in km.
        wavelength_units: Unit string written to the output attrs.

    Returns:
        Two-dimensional resolved wavelength in ``wavelength_units``.
        Locations with fewer than two trusted scale samples, or where
        the threshold isn't bracketed by the score column, are NaN.
    """
    score = wavelet_psd_score(
        pred,
        truth,
        scales,
        dim=dim,
        x0=x0,
        ntheta=ntheta,
        k0=k0,
        isotropic=True,
    )["score"]
    scale_dim = scales.dims[0]
    wavelengths = (
        1.0
        / np.asarray(
            scale_to_wavenumber(scales, x0=x0, k0=k0).values,
            dtype=float,
        )
        * wavelength_scale
    )
    out = xr.apply_ufunc(
        _resolved_scale_column,
        score,
        input_core_dims=[[scale_dim]],
        output_core_dims=[[]],
        kwargs={"wavelengths_km": wavelengths, "threshold": threshold},
        vectorize=True,
        dask="forbidden",
        output_dtypes=[float],
    )
    out.name = "wavelet_resolved_scale"
    out.attrs["units"] = wavelength_units
    out.attrs["threshold"] = threshold
    return out

Multiscale — region & band conditioned¶

Stratify skill by geographic region or by frequency band.

EvaluateByRegion ¶

EvaluateByRegion(metric: Operator, *, regions: Any)

Bases: Operator

Apply an inner metric per region.

Parameters:

Name	Type	Description	Default
`metric`	`Operator`	Inner :class:`Operator` evaluated per region.	required
`regions`	`Any`	One of: * a `regionmask.Regions` instance (requires the optional `regionmask` dep), * an integer-encoded mask :class:`xr.DataArray`, * a `dict[str, xr.DataArray[bool]]`.	required

Pre-normalize once and reuse via :func:normalize_regions if you intend to evaluate multiple metrics over the same regions:

```pycon
>>> mask, names = normalize_regions(natural_earth, ds)
>>> EvaluateByRegion(RMSE(...), regions=mask)
```

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/multiscale.py

def __init__(self, metric: Operator, *, regions: Any) -> None:
    if not isinstance(metric, Operator):
        raise TypeError(
            f"EvaluateByRegion requires an Operator, got {type(metric).__name__!r}."
        )
    self.metric = metric
    self.regions = regions

RegionScores ¶

RegionScores(*, var_ref: str, var_pred: str, regions: DataArray | Regions, lon: str = 'lon', lat: str = 'lat', metrics: Sequence[str] = ('rmse', 'bias', 'correlation', 'explained_variance'), region_dim: str = 'region')

Bases: Operator

Layer-1 wrapper for :func:scores_by_region.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/residuals.py

def __init__(
    self,
    *,
    var_ref: str,
    var_pred: str,
    regions: xr.DataArray | regionmask.Regions,
    lon: str = "lon",
    lat: str = "lat",
    metrics: Sequence[str] = ("rmse", "bias", "correlation", "explained_variance"),
    region_dim: str = "region",
) -> None:
    self.var_ref = var_ref
    self.var_pred = var_pred
    self.regions = regions
    self.lon = lon
    self.lat = lat
    self.metrics = tuple(metrics)
    self.region_dim = region_dim

FrequencyBandSkill ¶

FrequencyBandSkill(variable: str, dims: Sequence[str], bands: Mapping[str, tuple[float, float]], *, metric: Operator | None = None, coord_spacing: Mapping[str, float] | None = None)

Bases: Operator

Apply an inner metric to band-passed (pred, ref) pairs.

Parameters:

Name	Type	Description	Default
`variable`	`str`	Variable to score.	required
`dims`	`Sequence[str]`	Dimensions to FFT along (passed through to inner metric as well when `metric` is omitted).	required
`bands`	`Mapping[str, tuple[float, float]]`	`{name: (low, high)}` in physical coord units.	required
`metric`	`Operator \| None`	Inner :class:`Operator` (default `RMSE(variable, dims=dims)`).	`None`
`coord_spacing`	`Mapping[str, float] \| None`	Optional per-dim spacing override (in the same physical units as `bands`).	`None`

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/spectral.py

def __init__(
    self,
    variable: str,
    dims: Sequence[str],
    bands: Mapping[str, tuple[float, float]],
    *,
    metric: Operator | None = None,
    coord_spacing: Mapping[str, float] | None = None,
) -> None:
    if metric is None:
        from xrops.metrics._src.pixel import RMSE

        metric = RMSE(variable, dims=tuple(dims))
    elif not isinstance(metric, Operator):
        raise TypeError(
            "metric must be an Operator instance (e.g. RMSE(...)) so its "
            "configuration is introspectable."
        )
    _validate_bands(bands)
    self.variable = variable
    self.dims = tuple(dims)
    self.bands = {k: (float(v[0]), float(v[1])) for k, v in bands.items()}
    self.metric = metric
    self.coord_spacing = (
        None
        if coord_spacing is None
        else {k: float(v) for k, v in coord_spacing.items()}
    )

BandLimitedRMSE ¶

BandLimitedRMSE(variable: str, dims: Sequence[str], bands: Mapping[str, tuple[float, float]], *, coord_spacing: Mapping[str, float] | None = None)

Bases: FrequencyBandSkill

Fixed-metric convenience subclass: FrequencyBandSkill with RMSE.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/spectral.py

def __init__(
    self,
    variable: str,
    dims: Sequence[str],
    bands: Mapping[str, tuple[float, float]],
    *,
    coord_spacing: Mapping[str, float] | None = None,
) -> None:
    from xrops.metrics._src.pixel import RMSE

    super().__init__(
        variable,
        dims,
        bands,
        metric=RMSE(variable, dims=tuple(dims)),
        coord_spacing=coord_spacing,
    )

SkillByLeadTime ¶

SkillByLeadTime(metric: Operator, *, lead_dim: str = 'lead_time')

Bases: Operator

Apply an inner metric per lead-time slice.

Parameters:

Name	Type	Description	Default
`metric`	`Operator`	Inner :class:`Operator` (e.g. `RMSE("ssh", dims="time")`) evaluated on each `lead_dim` slice independently.	required
`lead_dim`	`str`	Name of the lead-time dimension on the prediction dataset. Defaults to `"lead_time"`.	`'lead_time'`

Example

>>> from xrops.metrics import RMSE, SkillByLeadTime
>>> op = SkillByLeadTime(RMSE("ssh", dims=("lat", "lon")))
>>> skill = op(pred_ds, ref_ds)  # DataArray indexed by lead_time

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/forecast.py

def __init__(self, metric: Operator, *, lead_dim: str = "lead_time") -> None:
    if not isinstance(metric, Operator):
        raise TypeError(
            f"SkillByLeadTime requires an Operator, got "
            f"{type(metric).__name__!r}; pass an instantiated metric like "
            "RMSE(...) so its config is introspectable."
        )
    self.metric = metric
    self.lead_dim = lead_dim

Physical — balance residuals¶

Conservation- and balance-law residuals used as physical evaluation diagnostics (geostrophic balance, divergence, potential-vorticity conservation, density inversions).

The derivative-based metrics default to lon/lat in degrees with latitude-derived Coriolis. Idealised f/β-plane model output in metres passes the horizontal dims, the geometry and an explicit Coriolis parameter (same f= forms as xrops.ocn):

from xrops.metrics import divergence_error, geostrophic_imbalance

cart = {"dims": ("x", "y"), "geometry": "cartesian"}
geostrophic_imbalance(run, ssh_var="eta", **cart, f=(1e-4, 1.6e-11))  # 0 = balanced
divergence_error(run, **cart, reduce="rms")                           # 1/s

geostrophic_imbalance is the scale-invariant ratio \(\mathrm{rms}(f\,\hat k\times\mathbf u + g\nabla\eta)\,/\,\mathrm{rms}(f\,\hat k\times\mathbf u)\): 0 for exact geostrophy, 1 for purely inertial flow.

GeostrophicBalanceError ¶

GeostrophicBalanceError(*, ssh_var: str = 'ssh', u_var: str = 'u', v_var: str = 'v', lat: str = 'lat', lon: str = 'lon', g: float = GRAVITY, dims: tuple[str, str] | None = None, geometry: Geometry = 'spherical', f: Coriolis = None)

Bases: Operator

Operator wrapper for :func:geostrophic_balance_error.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/physical.py

def __init__(
    self,
    *,
    ssh_var: str = "ssh",
    u_var: str = "u",
    v_var: str = "v",
    lat: str = "lat",
    lon: str = "lon",
    g: float = grad.GRAVITY,
    dims: tuple[str, str] | None = None,
    geometry: Geometry = "spherical",
    f: Coriolis = None,
) -> None:
    self.ssh_var = ssh_var
    self.u_var = u_var
    self.v_var = v_var
    self.lat = lat
    self.lon = lon
    self.g = g
    self.dims = dims
    self.geometry = geometry
    self.f = f

GeostrophicImbalance ¶

GeostrophicImbalance(*, ssh_var: str = 'ssh', u_var: str = 'u', v_var: str = 'v', dims: tuple[str, str] = ('lon', 'lat'), geometry: Geometry = 'spherical', f: Coriolis = None, g: float = GRAVITY, eps: float = 1e-30)

Bases: Operator

Operator wrapper for :func:geostrophic_imbalance.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/physical.py

def __init__(
    self,
    *,
    ssh_var: str = "ssh",
    u_var: str = "u",
    v_var: str = "v",
    dims: tuple[str, str] = ("lon", "lat"),
    geometry: Geometry = "spherical",
    f: Coriolis = None,
    g: float = grad.GRAVITY,
    eps: float = 1e-30,
) -> None:
    self.ssh_var = ssh_var
    self.u_var = u_var
    self.v_var = v_var
    self.dims = dims
    self.geometry = geometry
    self.f = f
    self.g = g
    self.eps = eps

DivergenceError ¶

DivergenceError(*, u_var: str = 'u', v_var: str = 'v', lat: str = 'lat', lon: str = 'lon', dims: tuple[str, str] | None = None, geometry: Geometry = 'spherical', reduce: Literal['rms'] | None = None)

Bases: Operator

Operator wrapper for :func:divergence_error.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/physical.py

def __init__(
    self,
    *,
    u_var: str = "u",
    v_var: str = "v",
    lat: str = "lat",
    lon: str = "lon",
    dims: tuple[str, str] | None = None,
    geometry: Geometry = "spherical",
    reduce: Literal["rms"] | None = None,
) -> None:
    self.u_var = u_var
    self.v_var = v_var
    self.lat = lat
    self.lon = lon
    self.dims = dims
    self.geometry = geometry
    self.reduce = reduce

PVConservationError ¶

PVConservationError(*, pv_var: str = 'pv', traj_dim: str = 'trajectory', time_dim: str = 'time')

Bases: Operator

Operator wrapper for :func:pv_conservation_error.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/physical.py

def __init__(
    self,
    *,
    pv_var: str = "pv",
    traj_dim: str = "trajectory",
    time_dim: str = "time",
) -> None:
    self.pv_var = pv_var
    self.traj_dim = traj_dim
    self.time_dim = time_dim

DensityInversionFraction ¶

DensityInversionFraction(*, density_var: str = 'rho', depth_dim: str = 'depth')

Bases: Operator

Operator wrapper for :func:density_inversion_fraction.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/physical.py

def __init__(self, *, density_var: str = "rho", depth_dim: str = "depth") -> None:
    self.density_var = density_var
    self.depth_dim = depth_dim

Masked — coverage aware¶

Wrap any metric so it ignores masked / missing cells and reports coverage.

MaskedMetric ¶

MaskedMetric(metric: Operator, mask: DataArray | None = None)

Bases: Operator

Wrap any inner metric so it operates only on masked-valid pixels.

Parameters:

Name	Type	Description	Default
`metric`	`Operator`	Inner two-input :class:`Operator` (e.g. `RMSE(...)`).	required
`mask`	`DataArray \| None`	Optional default mask. If `None`, a per-call mask must be supplied via `__call__(pred, ref, mask=...)`.	`None`

Example

>>> from xrops.metrics import RMSE, MaskedMetric
>>> op = MaskedMetric(RMSE("ssh", ("lat", "lon")), mask=ocean)
>>> op(pred_ds, ref_ds)
>>>
>>> op2 = MaskedMetric(RMSE("sst", ("lat", "lon")))
>>> op2(pred_ds, ref_ds, mask=cloud_free_today)

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/masked.py

def __init__(self, metric: Operator, mask: xr.DataArray | None = None) -> None:
    if not isinstance(metric, Operator):
        raise TypeError(
            f"MaskedMetric requires an Operator, got {type(metric).__name__!r}."
        )
    self.metric = metric
    self.mask = mask

BinnedResiduals2D ¶

BinnedResiduals2D(*, var_ref: str, var_pred: str, lon_bins: Sequence[float], lat_bins: Sequence[float], lon: str = 'lon', lat: str = 'lat', statistics: Sequence[str] = ('mean', 'std', 'count', 'rmse'))

Bases: Operator

Layer-1 wrapper for :func:bin_residuals_2d.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/residuals.py

def __init__(
    self,
    *,
    var_ref: str,
    var_pred: str,
    lon_bins: Sequence[float],
    lat_bins: Sequence[float],
    lon: str = "lon",
    lat: str = "lat",
    statistics: Sequence[str] = ("mean", "std", "count", "rmse"),
) -> None:
    self.var_ref = var_ref
    self.var_pred = var_pred
    self.lon_bins = list(lon_bins)
    self.lat_bins = list(lat_bins)
    self.lon = lon
    self.lat = lat
    self.statistics = tuple(statistics)

Distributional distances¶

Distance between predicted and reference distributions.

CRPS ¶

CRPS(variable: str, *, ensemble_dim: str = 'member', dims: Sequence[str] | None = None)

Bases: _DistributionalOp

Two-input CRPS operator.

Inputs: (ds_ensemble, ds_reference) where ds_ensemble has an ensemble_dim axis (one or more members) and ds_reference is deterministic. A one-member ensemble scores as the absolute error |x - y|; see :func:crps_ensemble.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/distributional.py

def __init__(
    self,
    variable: str,
    *,
    ensemble_dim: str = "member",
    dims: Sequence[str] | None = None,
) -> None:
    self.variable = variable
    self.ensemble_dim = ensemble_dim
    self.dims = None if dims is None else list(dims)

Wasserstein1 ¶

Wasserstein1(variable: str, *, sample_dim_a: str = 'member', sample_dim_b: str = 'member', dims: Sequence[str] | None = None)

Bases: _TwoSampleOp

1-D Wasserstein-distance operator.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/distributional.py

def __init__(
    self,
    variable: str,
    *,
    sample_dim_a: str = "member",
    sample_dim_b: str = "member",
    dims: Sequence[str] | None = None,
) -> None:
    self.variable = variable
    self.sample_dim_a = sample_dim_a
    self.sample_dim_b = sample_dim_b
    self.dims = None if dims is None else list(dims)

EnergyDistance ¶

EnergyDistance(variable: str, *, sample_dim_a: str = 'member', sample_dim_b: str = 'member', dims: Sequence[str] | None = None)

Bases: _TwoSampleOp

Two-sample energy-distance operator.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/distributional.py

def __init__(
    self,
    variable: str,
    *,
    sample_dim_a: str = "member",
    sample_dim_b: str = "member",
    dims: Sequence[str] | None = None,
) -> None:
    self.variable = variable
    self.sample_dim_a = sample_dim_a
    self.sample_dim_b = sample_dim_b
    self.dims = None if dims is None else list(dims)

Structural & probabilistic¶

SSIM ¶

SSIM(variable: str, dims: str | Sequence[str], *, window: int | None = None)

Bases: Operator

Structural-similarity-index operator (scikit-image).

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/structural.py

def __init__(
    self,
    variable: str,
    dims: str | Sequence[str],
    *,
    window: int | None = None,
) -> None:
    self.variable = variable
    self.dims = list(_normalize_dims(dims))
    self.window = window

GradientDifference ¶

GradientDifference(variable: str, dims: str | Sequence[str])

Bases: Operator

RMS gradient-difference operator.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/structural.py

def __init__(self, variable: str, dims: str | Sequence[str]) -> None:
    self.variable = variable
    self.dims = list(_normalize_dims(dims))

PhaseShiftError ¶

PhaseShiftError(variable: str, dims: str | Sequence[str], *, periodic: bool = False)

Bases: Operator

FFT cross-correlation phase-shift operator.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/structural.py

def __init__(
    self,
    variable: str,
    dims: str | Sequence[str],
    *,
    periodic: bool = False,
) -> None:
    self.variable = variable
    self.dims = list(_normalize_dims(dims))
    self.periodic = periodic

EnsembleCoverage ¶

EnsembleCoverage(variable: str, *, q: tuple[float, float] = (0.05, 0.95), ensemble_dim: str = 'member')

Bases: _ProbabilisticOp

Fraction of references inside the q ensemble envelope.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/probabilistic.py

def __init__(
    self,
    variable: str,
    *,
    q: tuple[float, float] = (0.05, 0.95),
    ensemble_dim: str = "member",
) -> None:
    super().__init__(variable, ensemble_dim=ensemble_dim)
    self.q = tuple(q)

SpreadSkillRatio ¶

SpreadSkillRatio(variable: str, *, ensemble_dim: str = 'member', dims: Sequence[str] | None = None)

Bases: _ProbabilisticOp

SpreadSkillRatio(variable, ensemble_dim, dims).

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/probabilistic.py

def __init__(
    self,
    variable: str,
    *,
    ensemble_dim: str = "member",
    dims: Sequence[str] | None = None,
) -> None:
    super().__init__(variable, ensemble_dim=ensemble_dim)
    self.dims = None if dims is None else list(dims)

RankHistogram ¶

RankHistogram(variable: str, *, ensemble_dim: str = 'member')

Bases: _ProbabilisticOp

Talagrand rank histogram.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/probabilistic.py

def __init__(
    self,
    variable: str,
    *,
    ensemble_dim: str = "member",
) -> None:
    self.variable = variable
    self.ensemble_dim = ensemble_dim

ReliabilityCurve ¶

ReliabilityCurve(variable: str, *, probability_bins: Float[ndarray, 'edges'] | None = None)

Bases: Operator

Reliability-diagram operator (probability vs event).

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/probabilistic.py

def __init__(
    self,
    variable: str,
    *,
    probability_bins: Float[np.ndarray, "edges"] | None = None,
) -> None:
    self.variable = variable
    self.probability_bins = (
        None if probability_bins is None else np.asarray(probability_bins)
    )

Object / event verification¶

InstanceMatcher ¶

InstanceMatcher(*, iou_threshold: float = 0.5, scores: ndarray | None = None)

Bases: Operator

Greedy IoU matching between predicted and reference mask stacks.

Parameters:

Name	Type	Description	Default
`iou_threshold`	`float`	Minimum IoU for a (pred, ref) pair to count as a match. Defaults to `0.5` (Pascal-VOC convention; the paper uses `0.1` for the much sparser MethaneSAT plume regime).	`0.5`
`scores`	`ndarray \| None`	Optional length-`N_pred` array of per-prediction scores. When `None`, predictions are ranked by area.	`None`

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/instance.py

def __init__(
    self,
    *,
    iou_threshold: float = 0.5,
    scores: np.ndarray | None = None,
) -> None:
    if not 0.0 < iou_threshold <= 1.0:
        raise ValueError("iou_threshold must be in (0, 1]")
    self.iou_threshold = float(iou_threshold)
    self.scores = None if scores is None else np.asarray(scores, dtype=float)

InstanceF1AtIoU ¶

InstanceF1AtIoU(*, iou_threshold: float = 0.5, scores: ndarray | None = None)

Bases: Operator

Instance precision / recall / F1 at a single IoU threshold.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/instance.py

def __init__(
    self,
    *,
    iou_threshold: float = 0.5,
    scores: np.ndarray | None = None,
) -> None:
    if not 0.0 < iou_threshold <= 1.0:
        raise ValueError("iou_threshold must be in (0, 1]")
    self.iou_threshold = float(iou_threshold)
    self.scores = None if scores is None else np.asarray(scores, dtype=float)

AveragePrecisionMatched ¶

AveragePrecisionMatched(*, iou_threshold: float = 0.5, scores: ndarray | None = None)

Bases: Operator

VOC-style Average Precision at a single IoU threshold.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/instance.py

def __init__(
    self,
    *,
    iou_threshold: float = 0.5,
    scores: np.ndarray | None = None,
) -> None:
    if not 0.0 < iou_threshold <= 1.0:
        raise ValueError("iou_threshold must be in (0, 1]")
    self.iou_threshold = float(iou_threshold)
    self.scores = None if scores is None else np.asarray(scores, dtype=float)

CentroidDisplacement ¶

CentroidDisplacement(dims: tuple[str, str] = ('lat', 'lon'))

Bases: Operator

Centroid-displacement operator on labelled-object datasets.

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/structural.py

def __init__(self, dims: tuple[str, str] = ("lat", "lon")) -> None:
    self.dims = tuple(dims)

MaskIoU ¶

Bases: Operator

Pairwise mask-IoU matrix between two (N, H, W) stacks.

Forecast comparison¶

DieboldMariano ¶

DieboldMariano(variable: str, *, dim: str | Sequence[str] | None = None, h: int = 1, alternative: str = 'two-sided', power: float = 2.0, hln_correction: bool = True)

Bases: Operator

Layer-1 Diebold–Mariano test on a forecast-error variable.

Selects variable from two input Datasets of forecast errors (residuals) and runs :func:dm_test, returning a Dataset with statistic and p_value.

Parameters:

Name	Type	Description	Default
`variable`	`str`	Error variable to select from each input Dataset.	required
`dim`	`str \| Sequence[str] \| None`	Dimension(s) to reduce. `None` flattens the full array.	`None`
`h`	`int`	Forecast horizon.	`1`
`alternative`	`str`	`"two-sided"`, `"less"`, or `"greater"`.	`'two-sided'`
`power`	`float`	Loss exponent applied to `\|e\|`.	`2.0`
`hln_correction`	`bool`	Apply the Harvey–Leybourne–Newbold correction.	`True`

Example

>>> op = DieboldMariano("error", dim="time")
>>> result = op(ds_errors_a, ds_errors_b)
>>> result["statistic"], result["p_value"]

Source code in packages/xrtoolz-ops/src/xrops/metrics/_src/dm.py

def __init__(
    self,
    variable: str,
    *,
    dim: str | Sequence[str] | None = None,
    h: int = 1,
    alternative: str = "two-sided",
    power: float = 2.0,
    hln_correction: bool = True,
) -> None:
    self.variable = variable
    self.dim = dim
    self.h = h
    self.alternative = alternative
    self.power = power
    self.hln_correction = hln_correction