Adding model vs obs additive bias

src/CSET/operators/scoreswrappers.py

@@ -245,6 +245,52 @@ def scores_rmse_model_obs(
     return rmse_cubes
 
 
+def scores_additive_bias_model_obs(
+    cubes: CubeList, preserved_coordinates: list[str] | str | None = None
+):
+    r"""Calculate the Additive Bias (Mean Error) using scores.
+
+    Acts as a wrapper around the ME calculation from ``scores`` ([scoresa]_, [scoresb]_).
+
+    Parameters
+    ----------
+    cubes: iris.cube.CubeList
+        A CubeList containing an observation cube and at least one model cube.
+    preserved_coordinates: list[str] | str | None, default is None.
+        The coordinates that you wish to preserve in the calculaiton of the
+        ME. For example if you want a map of each time you can preserve
+        ["time","latitude", "longitude"] or if you want a time series
+        you can preserve ["time"], if you want to collapse to a single value
+        use `None`. The default is `None`.
+
+    Returns
+    -------
+    scores_cube: iris.cube.CubeList
+        A cube list containing the ME between the models and observation cube.
+    """
+    additive_bias_cubes = CubeList()
+    model_list = CubeList()
+
+    for cb in cubes:
+        if "observed" in cb.long_name:
+            observed = cb
+        else:
+            model_list.append(cb)
+
+    for model in model_list:
+        input_cubelist = CubeList()
+        input_cubelist.append(observed)
+        input_cubelist.append(model)
+        additive_bias = scores_additive_bias(
+            input_cubelist, preserved_coordinates, obs_model_comparison=True
+        )
+        model_name = model.attributes["model_name"]
+        additive_bias.attributes["model_name"] = model_name
+        additive_bias_cubes.append(additive_bias)
+
+    return additive_bias_cubes
+
+
 def scores_rmse(
     cubes: CubeList,
     preserved_coordinates: list[str] | str | None = None,
@@ -404,7 +450,9 @@ def scores_mae(cubes: CubeList, preserved_coordinates: list[str] | str | None =
 
 
 def scores_additive_bias(
-    cubes: CubeList, preserved_coordinates: list[str] | str | None = None
+    cubes: CubeList,
+    preserved_coordinates: list[str] | str | None = None,
+    obs_model_comparison: bool = False,
 ):
     r"""Calculate the Additive Bias (Mean Error) using scores.
 
@@ -427,7 +475,14 @@ def scores_additive_bias(
     scores_cube: iris.cube.Cube
         A cube containing the ME between the base and other cube.
     """
-    base, other = _sort_cubes_for_verification(cubes)
+    if obs_model_comparison:
+        for cb in cubes:
+            if "observed" in cb.long_name:
+                base = cb
+            else:
+                other = cb
+    else:
+        base, other = _sort_cubes_for_verification(cubes)
 
     # Copy the coordinates of the input cubes.
     other_xr = xr.DataArray.from_iris(other)

src/CSET/recipes/verification/surface_scores_model_vs_obs_additive_bias.yaml

@@ -0,0 +1,75 @@
+category: Surface Model vs Observation
+title: "$VARNAME scores additive bias spatial plot for model vs observation points $SUBAREA_NAME"
+description: |
+  Extracts and plots the Mean Error (also known as the Additive Bias) in $METHOD of $VARNAME for model versus observation.
+  The ME is calculated based on that used in the
+  package [`scores`](https://scores.readthedocs.io/en/stable/api.html#scores.continuous.mean_error).
+  This recipe preserves the time, latitude and longitude coordinates.
+  This means a sequence of spatial plots of the ME can be produced.
+  Currently only working for deterministic models.
+
+  This recipe can preserve the time, latitude and longitude coordinates, in which case
+  a sequence of spatial plots of the ME are produced according to the averaging
+  method i.e. SEQ, MAX and MIN.
+
+  Alternatively, the ME can be computed over all time points in the case study, i.e.
+  preserving just the latitude and longitude coordinates by selecting the method CASE.
+  This will produce a single output plot of the total ME.
+
+  A larger ME implies a greater error than a smaller ME. Although a near zero ME could be indicative
+  of a good match between two fields, this could also arise if there are large compensating
+  positive and negative  errors. The ME is calculated on the grid point and thus a spatial view of
+  the ME provides useful information about whether a field has a bias for over or under-forecasting.
+
+
+
+steps:
+  - operator: read.read_cubes
+    file_paths: $INPUT_PATHS
+    model_names: $MODEL_NAME
+    subarea_type: $SUBAREA_TYPE
+    subarea_extent: $SUBAREA_EXTENT
+    constraint:
+      operator: constraints.generate_var_constraint
+      varname: ['observed_$VARNAME', '$VARNAME']
+
+  - operator: filters.filter_multiple_cubes
+    constraint:
+        operator: constraints.generate_cell_methods_constraint
+        cell_methods: []
+
+  - 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_rmse_model_obs
+    preserved_coordinates: $PRESERVED_COORDS
+
+  - operator: collapse.collapse
+    coordinate: ["time"]
+    method: $METHOD
+
+  - operator: plot.spatial_pcolormesh_plot
+
+  - operator: write.write_cube_to_nc
+    overwrite: True

src/CSET/recipes/verification/timeseries_surface_scores_model_vs_obs_additive_bias.yaml

@@ -0,0 +1,61 @@
+category: Surface Model verus Observation Time Series
+title: "$VARNAME scores additive bias time series for model vs observation points $SUBAREA_NAME"
+description: |
+  Extracts and plots the Mean Error (also known as the Additive Bias) in $VARNAME
+  computed over the domain for each timestep. The ME is calculated based on that used in the
+  package [`scores`](https://scores.readthedocs.io/en/stable/api.html#scores.continuous.mean_error).
+  This recipe allows the preservation of the time coordinate to produce a timeseries. Therefore, the ME is
+  collapsed over all other coordinates in the cube and calculated for every timestep.
+
+  A larger ME implies a greater error than a smaller ME. Although a near zero ME could be indicative
+  of a good match between two fields, this could also arise if there are large compensating positive and
+  negative  errors. The ME is calculated on the grid point and thus a spatial view of the ME provides useful
+  information about whether a field has a bias for over or under-forecasting.
+
+
+steps:
+  - operator: read.read_cubes
+    file_paths: $INPUT_PATHS
+    model_names: $MODEL_NAME
+    subarea_type: $SUBAREA_TYPE
+    subarea_extent: $SUBAREA_EXTENT
+    constraint:
+      operator: constraints.generate_var_constraint
+      varname: ['observed_$VARNAME', '$VARNAME']
+
+  - operator: filters.filter_multiple_cubes
+    constraint:
+        operator: constraints.generate_cell_methods_constraint
+        cell_methods: []
+
+  - 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_additive_bias_model_obs
+    preserved_coordinates: "time"
+
+  - operator: plot.plot_line_series
+
+  - operator: write.write_cube_to_nc
+    overwrite: True

tests/operators/test_scoreswrappers.py

@@ -344,6 +344,51 @@ def test_model_obs_rmse_preserve_in_timelatlon(dummy_cubelist_model_obs):
         assert cube.attributes["model_name"] == model_name
 
 
+def test_model_obs_additive_bias_preserve_in_time(dummy_cubelist_model_obs):
+    """Additive bias collapsed over station, preserving only the time dimension."""
+    bias = scoreswrappers.scores_additive_bias_model_obs(
+        dummy_cubelist_model_obs, "time"
+    )
+    assert isinstance(bias, CubeList)
+    assert len(bias) == 2
+    model_names = ["model_a", "model_b"]
+    for cube, model_name in zip(bias, model_names, strict=True):
+        assert cube.name() == "Additive_Bias_of_observed_temperature_at_screen_level"
+        assert cube.units == "K"
+        assert cube.shape == (36,)
+        assert cube.attributes["model_name"] == model_name
+
+
+def test_model_obs_additive_bias_preserve_in_latlon(dummy_cubelist_model_obs):
+    """Additive bias collapsed over time, preserving the station dimension via lat/lon."""
+    bias = scoreswrappers.scores_additive_bias_model_obs(
+        dummy_cubelist_model_obs, ["longitude", "latitude"]
+    )
+    assert isinstance(bias, CubeList)
+    assert len(bias) == 2
+    model_names = ["model_a", "model_b"]
+    for cube, model_name in zip(bias, model_names, strict=True):
+        assert cube.name() == "Additive_Bias_of_observed_temperature_at_screen_level"
+        assert cube.units == "K"
+        assert cube.shape == (28,)
+        assert cube.attributes["model_name"] == model_name
+
+
+def test_model_obs_additive_bias_preserve_in_timelatlon(dummy_cubelist_model_obs):
+    """Additive bias with nothing collapsed, preserving both time and station dimensions."""
+    bias = scoreswrappers.scores_additive_bias_model_obs(
+        dummy_cubelist_model_obs, ["time", "longitude", "latitude"]
+    )
+    assert isinstance(bias, CubeList)
+    assert len(bias) == 2
+    model_names = ["model_a", "model_b"]
+    for cube, model_name in zip(bias, model_names, strict=True):
+        assert cube.name() == "Additive_Bias_of_observed_temperature_at_screen_level"
+        assert cube.units == "K"
+        assert cube.shape == (36, 28)
+        assert cube.attributes["model_name"] == model_name
+
+
 def test_pod_gt_2x2_manual_case(make_cube_categorical_testing):
     """Test basic 2x2 case for manual validation."""
     obs = make_cube_categorical_testing(