Skip to content

Spatial & Ensemble Metrics

Spatial verification and ensemble analysis metrics.

Spatial and Ensemble Metrics for Atmospheric Sciences (Aero Protocol Compliant)

BSS(obs, mod, threshold, dim=None, axis=None)

Brier Skill Score (BSS) for probabilistic forecasts (Aero Protocol).

Typical Use Cases

  • Evaluating the accuracy of probabilistic binary forecasts relative to climatology.
  • Common in meteorological verification for event occurrence.

Parameters

obs : xarray.DataArray or numpy.ndarray Observed binary outcomes (0 or 1) or continuous values (will be binarized). mod : xarray.DataArray or numpy.ndarray Forecast probabilities (0 to 1) or continuous values (will be binarized). threshold : float Threshold for converting values to binary events. dim : str or iterable of str, optional Dimension(s) along which to compute the score (xarray only). axis : int or iterable of int, optional Axis or axes along which to compute the score (numpy only).

Returns

xarray.DataArray, numpy.ndarray, or float Brier Skill Score.

Source code in src/monet_stats/spatial_ensemble_metrics.py 
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
def BSS(
    obs: Union[xr.DataArray, np.ndarray],
    mod: Union[xr.DataArray, np.ndarray],
    threshold: float,
    dim: Optional[Union[str, Iterable[str]]] = None,
    axis: Optional[Union[int, Iterable[int]]] = None,
) -> Union[xr.DataArray, np.ndarray, float]:
    """
    Brier Skill Score (BSS) for probabilistic forecasts (Aero Protocol).

    Typical Use Cases
    -----------------
    - Evaluating the accuracy of probabilistic binary forecasts relative to climatology.
    - Common in meteorological verification for event occurrence.

    Parameters
    ----------
    obs : xarray.DataArray or numpy.ndarray
        Observed binary outcomes (0 or 1) or continuous values (will be binarized).
    mod : xarray.DataArray or numpy.ndarray
        Forecast probabilities (0 to 1) or continuous values (will be binarized).
    threshold : float
        Threshold for converting values to binary events.
    dim : str or iterable of str, optional
        Dimension(s) along which to compute the score (xarray only).
    axis : int or iterable of int, optional
        Axis or axes along which to compute the score (numpy only).

    Returns
    -------
    xarray.DataArray, numpy.ndarray, or float
        Brier Skill Score.
    """
    if isinstance(obs, xr.DataArray) and isinstance(mod, xr.DataArray):
        obs, mod = xr.align(obs, mod, join="inner")
        # Binarize if not already
        o_bin = (obs >= threshold).astype(float)
        m_prob = (mod >= threshold).astype(float)

        bs = ((m_prob - o_bin) ** 2).mean(dim=dim)
        obs_clim = o_bin.mean(dim=dim)
        bs_ref = ((obs_clim - o_bin) ** 2).mean(dim=dim)

        res = xr.where(bs_ref != 0, 1.0 - (bs / bs_ref), 0.0)
        return _update_history(res, "Brier Skill Score (BSS)")

    o = np.asarray(obs)
    m = np.asarray(mod)
    o_bin = (o >= threshold).astype(float)
    m_prob = (m >= threshold).astype(float)

    bs = np.nanmean((m_prob - o_bin) ** 2, axis=axis)
    obs_clim = np.nanmean(o_bin, axis=axis)
    if axis is not None:
        # Keep dimensions for subtraction
        obs_clim_kd = np.nanmean(o_bin, axis=axis, keepdims=True)
    else:
        obs_clim_kd = obs_clim

    bs_ref = np.nanmean((obs_clim_kd - o_bin) ** 2, axis=axis)

    with np.errstate(divide="ignore", invalid="ignore"):
        res = np.where(bs_ref != 0, 1.0 - (bs / bs_ref), 0.0)
        return res.item() if np.ndim(res) == 0 else res

CRPS(ensemble, obs, axis=0)

Continuous Ranked Probability Score (CRPS) for ensemble forecasts.

Supports lazy evaluation via Xarray/Dask.

Parameters

ensemble : xarray.DataArray or numpy.ndarray Ensemble forecasts. If DataArray, should have an ensemble dimension. obs : xarray.DataArray or numpy.ndarray Observed values. axis : int or str, optional Axis or dimension corresponding to ensemble members. Default is 0.

Returns

xarray.DataArray or numpy.ndarray CRPS values.

Examples

import numpy as np ens = np.array([[1, 2], [2, 3], [3, 4]]) obs = np.array([2, 3]) CRPS(ens, obs, axis=0) array([0.22222222, 0.22222222])

Source code in src/monet_stats/spatial_ensemble_metrics.py 
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
def CRPS(
    ensemble: Union[xr.DataArray, np.ndarray],
    obs: Union[xr.DataArray, np.ndarray],
    axis: Union[int, str] = 0,
) -> Union[xr.DataArray, np.ndarray]:
    """
    Continuous Ranked Probability Score (CRPS) for ensemble forecasts.

    Supports lazy evaluation via Xarray/Dask.

    Parameters
    ----------
    ensemble : xarray.DataArray or numpy.ndarray
        Ensemble forecasts. If DataArray, should have an ensemble dimension.
    obs : xarray.DataArray or numpy.ndarray
        Observed values.
    axis : int or str, optional
        Axis or dimension corresponding to ensemble members. Default is 0.

    Returns
    -------
    xarray.DataArray or numpy.ndarray
        CRPS values.

    Examples
    --------
    >>> import numpy as np
    >>> ens = np.array([[1, 2], [2, 3], [3, 4]])
    >>> obs = np.array([2, 3])
    >>> CRPS(ens, obs, axis=0)
    array([0.22222222, 0.22222222])
    """

    def _crps_numpy(ens, observation, ens_axis=0):
        ens_sorted = np.sort(ens, axis=ens_axis)
        n = ens.shape[ens_axis]
        # Compute empirical CDFs
        cdf_ens = np.arange(1, n + 1) / n
        shape = [1] * ens.ndim
        shape[ens_axis] = n
        cdf_ens = np.reshape(cdf_ens, shape)
        # Broadcast obs for comparison
        obs_broadcast = np.expand_dims(observation, ens_axis)
        cdf_obs = (ens_sorted >= obs_broadcast).astype(float)
        return np.sum((cdf_ens - cdf_obs) ** 2, axis=ens_axis)

    if isinstance(ensemble, xr.DataArray) and isinstance(obs, xr.DataArray):
        # Determine core dimension
        if isinstance(axis, int):
            ens_dim = ensemble.dims[axis]
        else:
            ens_dim = axis

        res = xr.apply_ufunc(
            _crps_numpy,
            ensemble,
            obs,
            input_core_dims=[[ens_dim], []],
            output_core_dims=[[]],
            kwargs={"ens_axis": -1},
            dask="parallelized",
            output_dtypes=[float],
            dask_gufunc_kwargs={"allow_rechunk": True},
        )
        return _update_history(res, "Continuous Ranked Probability Score (CRPS)")

    return _crps_numpy(np.asarray(ensemble), np.asarray(obs), ens_axis=axis)

EDS(obs, mod, threshold, dim=None, axis=None)

Extreme Dependency Score (EDS) for rare event detection (Aero Protocol).

Typical Use Cases

  • Assessing model performance for rare extreme events (e.g., heavy precipitation).
  • Used when traditional scores like CSI or ETS go to zero as the event becomes rarer.

Parameters

obs : xarray.DataArray or numpy.ndarray Observed field. mod : xarray.DataArray or numpy.ndarray Model field. threshold : float Event threshold to define the extreme event. dim : str or iterable of str, optional Dimension(s) along which to compute the score (xarray only). If None, reduces over all dimensions. axis : int or iterable of int, optional Axis or axes along which to compute the score (numpy only).

Returns

xarray.DataArray, numpy.ndarray, or float Extreme Dependency Score.

Examples

import numpy as np obs = np.zeros((10, 10)); obs[5, 5] = 1 mod = np.zeros((10, 10)); mod[5, 5] = 1 EDS(obs, mod, threshold=0.5) 1.0

Source code in src/monet_stats/spatial_ensemble_metrics.py 
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
def EDS(
    obs: Union[xr.DataArray, np.ndarray],
    mod: Union[xr.DataArray, np.ndarray],
    threshold: float,
    dim: Optional[Union[str, Iterable[str]]] = None,
    axis: Optional[Union[int, Iterable[int]]] = None,
) -> Union[xr.DataArray, np.ndarray, float]:
    """
    Extreme Dependency Score (EDS) for rare event detection (Aero Protocol).

    Typical Use Cases
    -----------------
    - Assessing model performance for rare extreme events (e.g., heavy precipitation).
    - Used when traditional scores like CSI or ETS go to zero as the event becomes rarer.

    Parameters
    ----------
    obs : xarray.DataArray or numpy.ndarray
        Observed field.
    mod : xarray.DataArray or numpy.ndarray
        Model field.
    threshold : float
        Event threshold to define the extreme event.
    dim : str or iterable of str, optional
        Dimension(s) along which to compute the score (xarray only).
        If None, reduces over all dimensions.
    axis : int or iterable of int, optional
        Axis or axes along which to compute the score (numpy only).

    Returns
    -------
    xarray.DataArray, numpy.ndarray, or float
        Extreme Dependency Score.

    Examples
    --------
    >>> import numpy as np
    >>> obs = np.zeros((10, 10)); obs[5, 5] = 1
    >>> mod = np.zeros((10, 10)); mod[5, 5] = 1
    >>> EDS(obs, mod, threshold=0.5)
    1.0
    """
    if isinstance(obs, xr.DataArray) and isinstance(mod, xr.DataArray):
        obs, mod = xr.align(obs, mod, join="inner")
        obs_bin = obs >= threshold
        mod_bin = mod >= threshold

        mask = obs.notnull() & mod.notnull()
        obs_bin = obs_bin & mask
        mod_bin = mod_bin & mask

        hits = (obs_bin & mod_bin).sum(dim=dim)
        n_obs = obs_bin.sum(dim=dim)
        n_mod = mod_bin.sum(dim=dim)
        n = mask.sum(dim=dim)

        # EDS = log(hits/n) / log(n_obs*n_mod / n^2)
        # Avoid log(0) and division by zero
        # p = n_obs / n, q = n_mod / n
        with np.errstate(divide="ignore", invalid="ignore"):
            eds = xr.where(
                (hits > 0) & (n_obs > 0) & (n_mod > 0), np.log(hits / n) / np.log((n_obs / n) * (n_mod / n)), np.nan
            )

        return _update_history(eds, "Extreme Dependency Score (EDS)")

    # NumPy path
    obs_arr = np.asarray(obs)
    mod_arr = np.asarray(mod)
    mask = ~np.isnan(obs_arr) & ~np.isnan(mod_arr)

    obs_bin = (obs_arr >= threshold) & mask
    mod_bin = (mod_arr >= threshold) & mask

    hits = np.sum(obs_bin & mod_bin, axis=axis)
    n_obs = np.sum(obs_bin, axis=axis)
    n_mod = np.sum(mod_bin, axis=axis)
    n = np.sum(mask, axis=axis)

    with np.errstate(divide="ignore", invalid="ignore"):
        res = np.where(
            (hits > 0) & (n_obs > 0) & (n_mod > 0), np.log(hits / n) / np.log((n_obs / n) * (n_mod / n)), np.nan
        )
        return res.item() if np.ndim(res) == 0 else res

SAL(obs, mod, threshold=None, lat_dim='lat', lon_dim='lon')

Structure-Amplitude-Location (SAL) score for spatial verification (Aero Protocol).

This implementation is vectorized over non-spatial dimensions using xarray.apply_ufunc and supports lazy evaluation for dimensions other than the spatial ones.

Parameters

obs : xarray.DataArray or numpy.ndarray Observed field. Should be 2D (lat, lon) or multi-dimensional. mod : xarray.DataArray or numpy.ndarray Model field. threshold : float, optional Threshold for object identification. If None, uses mean of observations per slice (Lazy-friendly). lat_dim : str, optional Name of the latitude dimension. Default is 'lat'. lon_dim : str, optional Name of the longitude dimension. Default is 'lon'.

Returns

S, A, L : xarray.DataArray, numpy.ndarray, or float Structure, Amplitude, and Location components.

Examples

import xarray as xr import numpy as np obs = xr.DataArray(np.random.rand(10, 10, 10), dims=['time', 'lat', 'lon']) mod = xr.DataArray(np.random.rand(10, 10, 10), dims=['time', 'lat', 'lon']) S, A, L = SAL(obs, mod)

Source code in src/monet_stats/spatial_ensemble_metrics.py 
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
def SAL(
    obs: Union[xr.DataArray, np.ndarray],
    mod: Union[xr.DataArray, np.ndarray],
    threshold: Optional[float] = None,
    lat_dim: str = "lat",
    lon_dim: str = "lon",
) -> Union[
    Tuple[xr.DataArray, xr.DataArray, xr.DataArray],
    Tuple[np.ndarray, np.ndarray, np.ndarray],
    Tuple[float, float, float],
]:
    """
    Structure-Amplitude-Location (SAL) score for spatial verification (Aero Protocol).

    This implementation is vectorized over non-spatial dimensions using `xarray.apply_ufunc`
    and supports lazy evaluation for dimensions other than the spatial ones.

    Parameters
    ----------
    obs : xarray.DataArray or numpy.ndarray
        Observed field. Should be 2D (lat, lon) or multi-dimensional.
    mod : xarray.DataArray or numpy.ndarray
        Model field.
    threshold : float, optional
        Threshold for object identification. If None, uses mean of observations
        per slice (Lazy-friendly).
    lat_dim : str, optional
        Name of the latitude dimension. Default is 'lat'.
    lon_dim : str, optional
        Name of the longitude dimension. Default is 'lon'.

    Returns
    -------
    S, A, L : xarray.DataArray, numpy.ndarray, or float
        Structure, Amplitude, and Location components.

    Examples
    --------
    >>> import xarray as xr
    >>> import numpy as np
    >>> obs = xr.DataArray(np.random.rand(10, 10, 10), dims=['time', 'lat', 'lon'])
    >>> mod = xr.DataArray(np.random.rand(10, 10, 10), dims=['time', 'lat', 'lon'])
    >>> S, A, L = SAL(obs, mod)
    """
    is_xr = isinstance(obs, xr.DataArray) and isinstance(mod, xr.DataArray)

    if is_xr:
        obs, mod = xr.align(obs, mod, join="inner")
        # Determine spatial dimensions
        if lat_dim not in obs.dims or lon_dim not in obs.dims:
            # Fallback to last two dimensions if lat/lon not found
            spatial_dims = [obs.dims[-2], obs.dims[-1]]
        else:
            spatial_dims = [lat_dim, lon_dim]
    else:
        # For numpy, assume last two dimensions are spatial
        spatial_dims = [-2, -1]

    def _sal_wrapper(o: np.ndarray, m: np.ndarray, t: Optional[float]) -> Tuple[float, float, float]:
        """
        Wrapper for _sal_numpy to be used with xarray.apply_ufunc.

        Parameters
        ----------
        o : numpy.ndarray
            Observed field slice.
        m : numpy.ndarray
            Model field slice.
        t : float, optional
            Threshold.

        Returns
        -------
        Tuple[float, float, float]
            S, A, L components.

        Examples
        --------
        >>> import numpy as np
        >>> o = np.random.rand(10, 10)
        >>> m = np.random.rand(10, 10)
        >>> _sal_wrapper(o, m, 0.5)
        """
        return _sal_numpy(o, m, t)

    res_s, res_a, res_l = xr.apply_ufunc(
        _sal_wrapper,
        obs,
        mod,
        threshold,
        input_core_dims=[spatial_dims, spatial_dims, []],
        output_core_dims=[[], [], []],
        vectorize=True,
        dask="parallelized",
        output_dtypes=[float, float, float],
    )

    if is_xr:
        res_s = _update_history(res_s, "SAL: Structure component")
        res_a = _update_history(res_a, "SAL: Amplitude component")
        res_l = _update_history(res_l, "SAL: Location component")

    return res_s, res_a, res_l

ensemble_mean(ensemble, axis=0)

Calculate the ensemble mean.

Parameters

ensemble : xarray.DataArray or numpy.ndarray Ensemble forecasts. axis : int or str, optional Axis or dimension corresponding to ensemble members. Default is 0.

Returns

xarray.DataArray or numpy.ndarray Ensemble mean.

Source code in src/monet_stats/spatial_ensemble_metrics.py 
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
def ensemble_mean(
    ensemble: Union[xr.DataArray, np.ndarray],
    axis: Union[int, str] = 0,
) -> Union[xr.DataArray, np.ndarray]:
    """
    Calculate the ensemble mean.

    Parameters
    ----------
    ensemble : xarray.DataArray or numpy.ndarray
        Ensemble forecasts.
    axis : int or str, optional
        Axis or dimension corresponding to ensemble members. Default is 0.

    Returns
    -------
    xarray.DataArray or numpy.ndarray
        Ensemble mean.
    """
    if isinstance(ensemble, xr.DataArray):
        dim = axis
        if isinstance(axis, int):
            dim = ensemble.dims[axis]
        res = ensemble.mean(dim=dim)
        return _update_history(res, "Ensemble Mean")
    return np.mean(ensemble, axis=axis)

ensemble_std(ensemble, axis=0)

Calculate the ensemble standard deviation.

Parameters

ensemble : xarray.DataArray or numpy.ndarray Ensemble forecasts. axis : int or str, optional Axis or dimension corresponding to ensemble members. Default is 0.

Returns

xarray.DataArray or numpy.ndarray Ensemble standard deviation.

Source code in src/monet_stats/spatial_ensemble_metrics.py 
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
def ensemble_std(
    ensemble: Union[xr.DataArray, np.ndarray],
    axis: Union[int, str] = 0,
) -> Union[xr.DataArray, np.ndarray]:
    """
    Calculate the ensemble standard deviation.

    Parameters
    ----------
    ensemble : xarray.DataArray or numpy.ndarray
        Ensemble forecasts.
    axis : int or str, optional
        Axis or dimension corresponding to ensemble members. Default is 0.

    Returns
    -------
    xarray.DataArray or numpy.ndarray
        Ensemble standard deviation.
    """
    if isinstance(ensemble, xr.DataArray):
        dim = axis
        if isinstance(axis, int):
            dim = ensemble.dims[axis]
        res = ensemble.std(dim=dim)
        return _update_history(res, "Ensemble Standard Deviation")
    return np.std(ensemble, axis=axis)

rank_histogram(ensemble, obs, axis=0)

Calculate the rank histogram counts.

Parameters

ensemble : xarray.DataArray or numpy.ndarray Ensemble forecasts. obs : xarray.DataArray or numpy.ndarray Observed values. axis : int or str, optional Axis or dimension corresponding to ensemble members. Default is 0.

Returns

xarray.DataArray or numpy.ndarray Rank histogram counts.

Examples

import numpy as np ens = np.array([[1, 2], [2, 3], [3, 4]]) obs = np.array([2, 3]) rank_histogram(ens, obs, axis=0) array([0., 0., 2., 0.])

Source code in src/monet_stats/spatial_ensemble_metrics.py 
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
def rank_histogram(
    ensemble: Union[xr.DataArray, np.ndarray],
    obs: Union[xr.DataArray, np.ndarray],
    axis: Union[int, str] = 0,
) -> Union[xr.DataArray, np.ndarray]:
    """
    Calculate the rank histogram counts.

    Parameters
    ----------
    ensemble : xarray.DataArray or numpy.ndarray
        Ensemble forecasts.
    obs : xarray.DataArray or numpy.ndarray
        Observed values.
    axis : int or str, optional
        Axis or dimension corresponding to ensemble members. Default is 0.

    Returns
    -------
    xarray.DataArray or numpy.ndarray
        Rank histogram counts.

    Examples
    --------
    >>> import numpy as np
    >>> ens = np.array([[1, 2], [2, 3], [3, 4]])
    >>> obs = np.array([2, 3])
    >>> rank_histogram(ens, obs, axis=0)
    array([0., 0., 2., 0.])
    """

    def _rank_numpy(ens, observation, ens_axis=0):
        o_exp = np.expand_dims(observation, ens_axis)
        full_ensemble = np.concatenate([ens, o_exp], axis=ens_axis)
        ranks = np.argsort(full_ensemble, axis=ens_axis)
        obs_rank = np.argmax(ranks == ens.shape[ens_axis], axis=ens_axis)
        n_ens = ens.shape[ens_axis]
        hist, _ = np.histogram(obs_rank, bins=np.arange(n_ens + 2))
        return hist.astype(float)

    if isinstance(ensemble, xr.DataArray) and isinstance(obs, xr.DataArray):
        if isinstance(axis, int):
            ens_dim = ensemble.dims[axis]
        else:
            ens_dim = axis

        def _rank_ufunc(ens, observation):
            o_exp = np.expand_dims(observation, -1)
            full_ensemble = np.concatenate([ens, o_exp], axis=-1)
            ranks = np.argsort(full_ensemble, axis=-1)
            return np.argmax(ranks == ens.shape[-1], axis=-1)

        obs_rank = xr.apply_ufunc(
            _rank_ufunc,
            ensemble,
            obs,
            input_core_dims=[[ens_dim], []],
            output_core_dims=[[]],
            dask="parallelized",
            output_dtypes=[int],
        )

        n_ens = ensemble.sizes[ens_dim]
        bins = np.arange(n_ens + 2)
        if hasattr(obs_rank.data, "dask"):
            import dask.array as da

            hist, _ = da.histogram(obs_rank.data, bins=bins)
            res = xr.DataArray(hist, dims="rank", coords={"rank": np.arange(n_ens + 1)})
        else:
            hist, _ = np.histogram(obs_rank.values, bins=bins)
            res = xr.DataArray(hist, dims="rank", coords={"rank": np.arange(n_ens + 1)})
        return _update_history(res, "Rank Histogram")

    return _rank_numpy(np.asarray(ensemble), np.asarray(obs), ens_axis=axis)

spread_error(ensemble, obs, axis=0, dim=None, reduce_axis=None)

Spread-Error Relationship for ensemble forecasts (Aero Protocol).

Typical Use Cases

  • Assessing if the ensemble spread is a good proxy for the forecast error.
  • Ideally, mean spread should equal RMSE of the ensemble mean.

Parameters

ensemble : xarray.DataArray or numpy.ndarray Ensemble forecasts. obs : xarray.DataArray or numpy.ndarray Observed values. axis : int or str, optional Axis or dimension corresponding to ensemble members. Default is 0. dim : str or iterable of str, optional Dimension(s) along which to compute the mean spread and error (xarray only). If None, reduces over all dimensions. reduce_axis : int or iterable of int, optional Axis or axes along which to compute the mean spread and error (numpy only).

Returns

mean_spread : float, numpy.ndarray, or xarray.DataArray Mean ensemble spread. mean_error : float, numpy.ndarray, or xarray.DataArray Mean absolute error of ensemble mean vs. obs.

Source code in src/monet_stats/spatial_ensemble_metrics.py 
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
def spread_error(
    ensemble: Union[xr.DataArray, np.ndarray],
    obs: Union[xr.DataArray, np.ndarray],
    axis: Union[int, str] = 0,
    dim: Optional[Union[str, Iterable[str]]] = None,
    reduce_axis: Optional[Union[int, Iterable[int]]] = None,
) -> Tuple[Any, Any]:
    """
    Spread-Error Relationship for ensemble forecasts (Aero Protocol).

    Typical Use Cases
    -----------------
    - Assessing if the ensemble spread is a good proxy for the forecast error.
    - Ideally, mean spread should equal RMSE of the ensemble mean.

    Parameters
    ----------
    ensemble : xarray.DataArray or numpy.ndarray
        Ensemble forecasts.
    obs : xarray.DataArray or numpy.ndarray
        Observed values.
    axis : int or str, optional
        Axis or dimension corresponding to ensemble members. Default is 0.
    dim : str or iterable of str, optional
        Dimension(s) along which to compute the mean spread and error (xarray only).
        If None, reduces over all dimensions.
    reduce_axis : int or iterable of int, optional
        Axis or axes along which to compute the mean spread and error (numpy only).

    Returns
    -------
    mean_spread : float, numpy.ndarray, or xarray.DataArray
        Mean ensemble spread.
    mean_error : float, numpy.ndarray, or xarray.DataArray
        Mean absolute error of ensemble mean vs. obs.
    """
    if isinstance(ensemble, xr.DataArray) and isinstance(obs, xr.DataArray):
        # Resolve ensemble dimension
        e_dim = axis
        if isinstance(axis, int):
            e_dim = ensemble.dims[axis]

        spread_field = ensemble.std(dim=e_dim)
        ens_mean_field = ensemble.mean(dim=e_dim)
        error_field = abs(ens_mean_field - obs)

        # Average over specified dimensions (or all remaining)
        m_spread = spread_field.mean(dim=dim)
        m_error = error_field.mean(dim=dim)

        return _update_history(m_spread, "Mean Ensemble Spread"), _update_history(m_error, "Mean Ensemble Error")

    # NumPy path
    ens = np.asarray(ensemble)
    observation = np.asarray(obs)

    # Calculate spread and ensemble mean
    spread_f = np.std(ens, axis=axis)
    ens_m_f = np.mean(ens, axis=axis)
    error_f = np.abs(ens_m_f - observation)

    # Average over specified axes
    m_spread = np.nanmean(spread_f, axis=reduce_axis)
    m_error = np.nanmean(error_f, axis=reduce_axis)

    if np.ndim(m_spread) == 0:
        return float(m_spread), float(m_error)
    return m_spread, m_error