2215 ensemble verification

src/CSET/cset_workflow/meta/verification/rose-meta.conf

@@ -73,3 +73,10 @@ description=Create a timeseries of the RMSE for the specified surface fields.
 type=python_boolean
 compulsory=true
 sort-key=scores2
+
+[template variables=SCORES_CRPS_FOR_ENSEMBLE]
+ns=Verification/Scores
+description=Create a timeseries of the CRPS for comparing ensemble members to control.
+type=python_boolean
+compulsory=true
+sort-key=scores3

src/CSET/cset_workflow/rose-suite.conf.example

@@ -127,6 +127,7 @@ RAIN_PRESENCE_SPATIAL_DIFFERENCE=False
 RAIN_PRESENCE_SPATIAL_PLOT=False
 SCORES_RMSE_SPATIAL=False
 SCORES_RMSE_TIMESERIES=False
+SCORES_CRPS_FOR_ENSEMBLE=False
 SCREEN_LEVEL_TEMPERATURE_PROBABILITIES=False
 !!SCREEN_LEVEL_TEMPERATURE_SPATIAL_PROBABILITY_WITHOUT_CONTROL_MEMBER=False
 SELECT_SUBAREA=False

src/CSET/loaders/verification.py

@@ -63,3 +63,17 @@ def load(conf: Config):
                 model_ids=[base_model["id"], model["id"]],
                 aggregation=False,
             )
+    if conf.SCORES_CRPS_FOR_ENSEMBLE:
+        for model, field in itertools.product(models, conf.SURFACE_FIELDS):
+            yield RawRecipe(
+                recipe="crps_for_ensemble.yaml",
+                variables={
+                    "VARNAME": field,
+                    "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None,
+                    "SUBAREA_EXTENT": conf.SUBAREA_EXTENT
+                    if conf.SELECT_SUBAREA
+                    else None,
+                },
+                model_ids=[model["id"]],
+                aggregation=False,
+            )

src/CSET/operators/scoreswrappers.py

@@ -21,13 +21,15 @@
 import numpy as np
 import scores
 import scores.continuous
+import scores.probability
 import xarray as xr
 from iris.cube import Cube, CubeList
 from iris.util import reverse
 
 from CSET._common import is_increasing
 from CSET.operators._utils import fully_equalise_attributes, get_cube_yxcoordname
 from CSET.operators.misc import _extract_common_time_points
+from CSET.operators.read import _realization_callback
 from CSET.operators.regrid import regrid_onto_cube
 
 
@@ -192,3 +194,57 @@
     )
     RMSE.rename(f"RMSE_of_{base.name()}")
     return RMSE
+
+
+def scores_crps_for_ensemble(cubes: Cube | CubeList):
+    r"""Calculate the CRPS for an ensemble.
+
+    Acts as a wrapper around the crps_for_ensemble from ``scores`` ([scoresa]_, [scoresb]_).
+
+    Parameters
+    ----------
+    cubes: iris.cube.Cube
+        A Cube containing ensembles data
+
+    Returns
+    -------
+    crps: iris.cube.Cube
+        A cube containing the crps between the ensemble members and the control
+
+    References
+    ----------
+    .. [scoresa] Leeuwenburg, T., Loveday, N., Ebert, E. E., Cook, H.,
+        Khanarmuei, M., Taggart, R. J., Ramanathan, N., Carroll, M., Chong, S.,
+        Griffiths, A., & Sharples, J. (2024) "scores: A Python package for
+        verifying and evaluating models and predictions with xarray". Journal
+        of Open Source Software, vol. 9, 6889. doi: 10.21105/joss.06889
+
+    .. [scoresb] Leeuwenburg, T., Loveday, N., Ramanathan, N., Chong, S.,
+        Taggart, R. J., Shrestha, D., Khanarmuei, M., Cook, H., Bluett, L., Ebert,
+        E. E., Carroll, M., Trotta, B., Bishop, S., Squire, D. T., Griffiths, A.,
+        Pagano, T. C., Fisher, A. J., Mandelbaum, T., Jinghan, F., … Smallwood, J.
+        (2026) "scores: Metrics for the verification, evaluation and optimisation of
+        forecasts, predictions or models (2.5.0)". Zenodo. doi: 10.5281/zenodo.18638494
+    """
+    ctrl = cubes[np.where(cubes.coord("realization").points[:] == 1)[0][0], :]
+    ens_mem = cubes[np.where(cubes.coord("realization").points[:] != 1)[0][:], :]
+
+    ctrl.data
+    ens_mem.data
+
+    ctrl = xr.DataArray.from_iris(ctrl)
+    ens_mem = xr.DataArray.from_iris(ens_mem)
+
+    crps = xr.DataArray.to_iris(
+        scores.probability.crps_for_ensemble(
+            ens_mem,
+            ctrl,
+            ensemble_member_dim="realization",
+            method="ecdf",
+            preserve_dims="time",
+        )
+    )
+
+    crps.rename(f"CRPS_of_{cubes[0].name()}")
+    _realization_callback(crps)
+    return crps

src/CSET/recipes/verification/crps_for_ensemble.yaml

@@ -0,0 +1,30 @@
+category: Scores
+title: "$VARNAME\nCRPS timeseries between ensemble members and control"
+description: |
+  Extracts and plots the  Continuous Ranked Probability Score (CRPS) for Ensembles in $VARNAME
+  for all times. The CRPS is calculated based on that used in the
+  package [`scores`](https://scores.readthedocs.io/en/stable/tutorials/CRPS_for_Ensembles.html).
+
+
+steps:
+  - operator: read.read_cube
+    file_paths: $INPUT_PATHS
+    constraint:
+      operator: constraints.combine_constraints
+      varname_constraint:
+        operator: constraints.generate_var_constraint
+        varname: $VARNAME
+      cell_methods_constraint:
+        operator: constraints.generate_cell_methods_constraint
+        cell_methods: []
+        varname: $VARNAME
+
+
+  - operator: scoreswrappers.scores_crps_for_ensemble
+
+  - operator: plot.plot_line_series
+    series_coordinate: time
+
+
+  - operator: write.write_cube_to_nc
+    overwrite: True