Scores [Categorical] Equitable Threat Score Model vs Obs

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

@@ -188,7 +188,7 @@ sort-key=sc-ts-3
 
 [template variables=SCORES_TIMESERIES_PC]
 ns=Verification/ScoresTimeseries
-title=Pearson's Correlation
+title=Pearsons Correlation
 description=SCORES COUNTINOUS METRIC.
             Plots a timeseries plot of the Pearson’s Correlation Coefficient between two models.
             The computation is implemented by the scores package.
@@ -320,16 +320,11 @@ type=python_boolean
 compulsory=true
 sort-key=sc-sp-1_model_vs_obs
 
-[Verification/Scores Categorical]
-[template variables=SCORES_CATEGORICAL_POD]
-ns=Verification/Scores Categorical
-description=Compute Scores Probability of Detection.
-type=python_boolean
-compulsory=true
-trigger=template variables=SCORES_CATEGORICAL_POD_ENTRIES: True
-sort-key=scoresPOD1
 
-[template variables=SCORES_CATEGORICAL_POD_ENTRIES]
+#####
+
+# Categorical Scores metrics against observations.
+[template variables=SCORES_CATEGORICAL_ENTRIES]
 ns=Verification/Scores Categorical
 description=List of variables, operators and thresholds
             i.e. ['<varname>,<operator>,<value>']
@@ -343,4 +338,20 @@ help=Supported operators:
      ['air_temperature,lt,280']
 type=python_list
 compulsory=true
+sort-key=scoresPOD1
+
+[Verification/Scores Categorical]
+[template variables=SCORES_CATEGORICAL_POD]
+ns=Verification/Scores Categorical
+description=Compute Scores Probability of Detection.
+type=python_boolean
+compulsory=true
 sort-key=scoresPOD2
+
+[Verification/Scores Categorical]
+[template variables=SCORES_CATEGORICAL_ETS]
+ns=Verification/Scores Categorical
+description=Compute Scores Equitable Threat Score.
+type=python_boolean
+compulsory=true
+sort-key=scoresPOD3

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

@@ -128,8 +128,9 @@ POINT_OBS=False
 !!POINT_OBS_USE_WMO_STATION_NUMBERS=False
 !!POINT_OBS_WMO_BLOCK_STTN_NUMBERS=[]
 PRESSURE_LEVELS=[]
+SCORES_CATEGORICAL_ETS=False
 SCORES_CATEGORICAL_POD=False
-!!SCORES_CATEGORICAL_POD_ENTRIES=[]
+SCORES_CATEGORICAL_ENTRIES=[]
 PRESSURE_LEVEL_FIELDS=[]
 !!PROB_TEMPERATURE_CONDITION=[]
 !!PROB_TEMPERATURE_THRESHOLD=[]

src/CSET/loaders/verification.py

@@ -69,6 +69,8 @@ def _get_scores_timeseries_categorical(conf):
     scores_timeseries_categorical = []
     if conf.SCORES_CATEGORICAL_POD or conf.SCORES_ALL:
         scores_timeseries_categorical.append("pod")
+    if conf.SCORES_CATEGORICAL_ETS or conf.SCORES_ALL:
+        scores_timeseries_categorical.append("ets")
     return scores_timeseries_categorical
 
 
@@ -291,7 +293,7 @@ def load(conf: Config):
     if scores_timeseries_categorical:
         # Produce timeseries plots of scores categorical metrics for each model.
         for field_and_method, scores_method in itertools.product(
-            conf.SCORES_CATEGORICAL_POD_ENTRIES, scores_timeseries_categorical
+            conf.SCORES_CATEGORICAL_ENTRIES, scores_timeseries_categorical
         ):
             try:
                 var, op, value = field_and_method.split(",")
@@ -308,8 +310,8 @@ def load(conf: Config):
                 variables={
                     "VARNAME": var,
                     "MODEL_NAME": ["OBS"] + [model["name"] for model in models],
-                    "POD_THRESHOLD": value,
-                    "POD_OPERATOR": op,
+                    "THRESHOLD": value,
+                    "OPERATOR": op,
                     "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "",
                     "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None,
                     "SUBAREA_EXTENT": conf.SUBAREA_EXTENT

src/CSET/operators/scoreswrappers.py

@@ -781,3 +781,114 @@ def scores_pod_model_obs(
         scores_results.append(scores_cube)
 
     return scores_results
+
+
+def scores_ets_model_obs(
+    cubes: CubeList,
+    preserved_coordinates: list[str] | str | None,
+    threshold: str,
+    op_func: str,
+):
+    r"""
+    Compute the Equitable Threat Score (ETS) score using Scores ([scoresa]_ [scoresb]_).
+
+    Parameters
+    ----------
+    cubes: iris.cube.CubeList
+        An iris cubelist containing model(s) and an observation cube.
+    preserved_coordinates: list | str | None
+        An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
+        then preserving coordinate 'time' will produce the equitable threat score for each timeslice (shape time). If None,
+        then it will return a single value score for all times/point locations.
+    threshold: str
+        A str containing the threshold to use to generate the binary masks, which subsequently gets turned to a float (but passed as str around the recipe templating).
+    op_func: str
+        A string either containing 'lt' for less than or 'gt for greater than, to determine how the threshold is applied to the data
+        to generate the mask.
+
+    Returns
+    -------
+    cube: iris.cube
+        An iris cube, containing the probability of detection score for further plotting.
+
+    Notes
+    -----
+    The Equitable Threat Score (ETS) evaluates the accuracy of forecasts for events that meet a specified threshold,
+    hile accounting for correct forecasts that could occur purely by chance. Unlike the Probability of Detection (POD),
+    ETS considers hits, misses, and false alarms, providing a more balanced assessment of forecast skill.
+
+    For example, if the threshold is 290 K and op_func is gt (greater than), an observation of 292 K and a forecast of 295 K
+    would be counted as a hit. ETS adjusts the total number of hits by removing the number of hits expected due to random chance.
+
+    It is calculated as:
+
+    .. math::
+
+    ETS = \frac{hits - hits_{random}}
+    {hits + misses + false\ alarms - hits_{random}}
+
+    where
+
+    hits_{random} = \frac{(hits + misses)(hits + false\ alarms)}{total count}
+
+    ETS ranges from -1/3 to 1, where 1 indicates a perfect forecast, 0 indicates no skill beyond random chance, and negative values indicate worse than
+    random chance.
+    """
+    # Split out model(s) and obs
+    models = CubeList()
+    for c in cubes:
+        if "observed" in c.long_name:
+            observed = c
+        else:
+            models.append(c)
+
+    # Setup cubelist to store results
+    scores_results = iris.cube.CubeList()
+
+    # Setup operators greater than, less than.
+    ops = {
+        "gt": operator.gt,
+        "lt": operator.lt,
+    }
+
+    try:
+        op = ops[op_func]
+    except KeyError as err:
+        raise ValueError(f"Operator {op_func} not supported.") from err
+
+    for model in models:
+        # Convert obs cubes to xarray and resolve preserved dimensions.
+        other_xr = xr.DataArray.from_iris(model)
+        base_xr = xr.DataArray.from_iris(observed)
+        preserve_dims = _resolve_preserve_dims(
+            observed, other_xr, preserved_coordinates
+        )
+
+        # Create event operator object using threshold and operator direction.
+        event_operator = scores.categorical.ThresholdEventOperator(
+            default_event_threshold=float(threshold), default_op_fn=op
+        )
+
+        # Generate binary fields using the event operator.
+        forecast_binary, observed_binary = event_operator.make_event_tables(
+            other_xr, base_xr
+        )
+
+        # Create binary contigency manager, as per Scores API, using transform to preserve preserve_dims
+        contingency_manager = scores.categorical.BinaryContingencyManager(
+            forecast_binary, observed_binary
+        ).transform(preserve_dims=preserve_dims)
+
+        # Get ETS from the contigency manager, and convert back to an iris cube.
+        scores_cube = xr.DataArray.to_iris(contingency_manager.equitable_threat_score())
+
+        # Rename cube so it plots correctly alongside correcting cube units.
+        scores_cube.rename(
+            f"Equitable_Threat_Score_{op_func}_{threshold}_{observed.name()}"
+        )
+        scores_cube.units = "1"
+        scores_cube.attributes["model_name"] = model.attributes["model_name"]
+
+        scores_results.append(scores_cube)
+
+    return scores_results

src/CSET/recipes/verification/surface_categorical_model_obs_ets.yaml

@@ -0,0 +1,82 @@
+category: Scores Categorical
+title: Timeseries of ETS for "$VARNAME over $SUBAREA_NAME, using threshold $OPERATOR $THRESHOLD"
+description: |
+
+    Extracts and plots the ETS for $VARNAME for each model against observations as a timeseries.
+
+    The Equitable Threat Score (ETS) in [`scores`](https://scores.readthedocs.io/en/stable/api.html#scores.categorical.BasicContingencyManager.equitable_threat_score)
+    evaluates the accuracy of forecasts for events that meet a specified threshold,
+    hile accounting for correct forecasts that could occur purely by chance. Unlike the Probability of Detection (POD),
+    ETS considers hits, misses, and false alarms, providing a more balanced assessment of forecast skill.
+
+    For example, if the threshold is 290 K and op_func is gt (greater than), an observation of 292 K and a forecast of 295 K
+    would be counted as a hit. ETS adjusts the total number of hits by removing the number of hits expected due to random chance.
+
+    It is calculated as:
+
+    ETS = frac{hits - hits_{random}}
+    {hits + misses + false\ alarms - hits_{random}}
+
+
+    where
+
+    hits_{random} = frac{(hits + misses)(hits + false\ alarms)}{total count}
+
+    ETS ranges from -1/3 to 1, where 1 indicates a perfect forecast, 0 indicates no skill beyond random chance, and negative values indicate worse than
+    random chance.
+
+    References
+    https://scores.readthedocs.io/en/stable/api.html#scores.categorical.BasicContingencyManager.equitable_threat_score
+
+steps:
+  - operator: read.read_cubes
+    file_paths: $INPUT_PATHS
+    model_names: $MODEL_NAME
+    constraint:
+      operator: constraints.combine_constraints
+      varname_constraint:
+        operator: constraints.generate_var_constraint
+        varname: ['observed_$VARNAME', '$VARNAME']
+      cell_methods_constraint:
+        operator: constraints.generate_cell_methods_constraint
+        cell_methods: []
+      pressure_level_constraint:
+        operator: constraints.generate_level_constraint
+        coordinate: "pressure"
+        levels: []
+    subarea_type: $SUBAREA_TYPE
+    subarea_extent: $SUBAREA_EXTENT
+
+  - operator: misc.extract_common_points
+    coordinate: time
+
+  - operator: misc.combine_cubes_into_cubelist
+    first:
+        operator: filters.filter_cubes
+        constraint:
+            operator: constraints.generate_var_constraint
+            varname: observed_$VARNAME
+    second:
+       operator: regrid.interpolate_to_point_cube
+       fld:
+           operator: filters.filter_multiple_cubes
+           constraint:
+               operator: constraints.combine_constraints
+               var_constraint:
+                   operator: constraints.generate_var_constraint
+                   varname: $VARNAME
+       point_cube:
+           operator: filters.filter_cubes
+           constraint:
+               operator: constraints.generate_var_constraint
+               varname: observed_$VARNAME
+
+  - operator: scoreswrappers.scores_ets_model_obs
+    preserved_coordinates: "time"
+    threshold: $THRESHOLD
+    op_func: $OPERATOR
+
+  - operator: plot.plot_line_series
+
+  - operator: write.write_cube_to_nc
+    overwrite: True

src/CSET/recipes/verification/surface_categorical_model_obs_pod.yaml

@@ -1,10 +1,11 @@
 category: Scores Categorical
-title: Timeseries of POD for "$VARNAME over $SUBAREA_NAME, using threshold $POD_OPERATOR $POD_THRESHOLD"
+title: Timeseries of POD for "$VARNAME over $SUBAREA_NAME, using threshold $OPERATOR $THRESHOLD"
 description: |
 
-  Create a timeseries of the probability of detection.
+  Extracts and plots the POD for $VARNAME for each model against observations as a timeseries.
 
-  The probability of detection calculates the proportion of observed events that meet a threshold that were correctly forecast by the model.
+  The probability of detection (POD)  in [`scores`](https://scores.readthedocs.io/en/stable/api.html#scores.categorical.BasicContingencyManager.probability_of_detection)
+  calculates the proportion of observed events that meet a threshold that were correctly forecast by the model.
   For example, if threshold is 290K and op_func is gt (greater than), and at some station a temperature was recorded as 292K and the model produced
   295k, that would be a positive hit. It does not take into account how far above/below a threshold a model forecasts.
 
@@ -62,8 +63,8 @@ steps:
 
   - operator: scoreswrappers.scores_pod_model_obs
     preserved_coordinates: "time"
-    threshold: $POD_THRESHOLD
-    op_func: $POD_OPERATOR
+    threshold: $THRESHOLD
+    op_func: $OPERATOR
 
   - operator: plot.plot_line_series