Adding model-obs pearson correlation recipe + operator

src/CSET/operators/scoreswrappers.py

@@ -326,6 +326,53 @@ def scores_additive_bias_model_obs(
     return additive_bias_cubes
 
 
+def scores_correlation_pearsonr_model_obs(
+    cubes: CubeList, preserved_coordinates: list[str] | str | None = None
+):
+    r"""Calculate the Pearson's Correlation (PC) coefficient using scores.
+
+    Acts as a wrapper around the PC calculation from ``scores`` ([scoresa]_, [scoresb]_).
+
+    Parameters
+    ----------
+    cubes: iris.cube.CubeList
+        A CubeList containing exactly two cubes: a base and an "other" model,
+        this can be an analysis and the model.
+    preserved_coordinates: list[str] | str | None, default is None.
+        The coordinates that you wish to preserve in the calculation of the
+        PC. 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 PC between the models and observation cube.
+    """
+    pearsonr_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)
+        pearsonr = scores_correlation_pearsonr(
+            input_cubelist, preserved_coordinates, obs_model_comparison=True
+        )
+        model_name = model.attributes["model_name"]
+        pearsonr.attributes["model_name"] = model_name
+        pearsonr_cubes.append(pearsonr)
+
+    return pearsonr_cubes
+
+
 def scores_rmse(
     cubes: CubeList,
     preserved_coordinates: list[str] | str | None = None,
@@ -590,7 +637,9 @@ def scores_additive_bias(
 
 
 def scores_correlation_pearsonr(
-    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 Pearson's Correlation (PC) coefficient using scores.
 
@@ -613,8 +662,16 @@ def scores_correlation_pearsonr(
     scores_cubelist: iris.cube.CubeList
         A cubelist containing the PC between the base and other cube(s).
     """
-    base, others = _sort_cube_into_base_and_other(cubes)
     scores_cubelist = CubeList()
+    if obs_model_comparison:
+        for cb in cubes:
+            if "observed" in cb.long_name:
+                base = cb
+            else:
+                others = [cb]
+    else:
+        base, others = _sort_cube_into_base_and_other(cubes)
+
     for other in others:
         base, other = _process_cubes_for_verification(base, other)
 

src/CSET/recipes/verification/timeseries_surface_scores_model_vs_obs_correlation_pearsonr.yaml

@@ -0,0 +1,60 @@
+category: Surface Model verus Observation Time Series
+title: "$VARNAME scores pearson correlation time series for model vs observation points $SUBAREA_NAME"
+description: |
+
+  Extracts and plots the Pearson's Correlation coefficient in $VARNAME
+  computed over the domain for each timestep. The PC is calculated based on that used in the
+  package [`scores`](https://scores.readthedocs.io/en/stable/api.html#scores.continuous.correlation.pearsonr).
+  This recipe allows the preservation of the time coordinate to produce a timeseries. Therefore, the PC is
+  collapsed over all other coordinates in the cube and calculated for every timestep.
+
+  The PC coefficient provides information on the linear relationship between two fields.
+  Perfectly correlated fields would yield a PC coefficient of 1, whereas perfectly anti-correlated
+  fields would yield a PC coefficient of -1. A value of 0 would indicate that the fields are non-correlated.
+
+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_correlation_pearsonr_model_obs
+    preserved_coordinates: "time"
+
+  - operator: plot.plot_line_series
+
+  - operator: write.write_cube_to_nc
+    overwrite: True

tests/operators/test_scoreswrappers.py

@@ -351,6 +351,50 @@ def test_model_obs_rmse_preserve_in_timelatlon(dummy_cubelist_model_obs):
         assert cube.attributes["model_name"] == model_name
 
 
+def test_model_obs_pearson_correlation_preserve_in_time(dummy_cubelist_model_obs):
+    """Pearson correlation collapsed over station, preserving only the time dimension."""
+    corr = scoreswrappers.scores_correlation_pearsonr_model_obs(
+        dummy_cubelist_model_obs, "time"
+    )
+    assert isinstance(corr, CubeList)
+    assert len(corr) == 2
+    model_names = ["model_a", "model_b"]
+    for cube, model_name in zip(corr, model_names, strict=True):
+        assert (
+            cube.name() == "Pearson_Correlation_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_pearson_correlation_preserve_in_latlon(dummy_cubelist_model_obs):
+    """Pearson correlation collapsed over time, preserving the station dimension via lat/lon."""
+    corr = scoreswrappers.scores_correlation_pearsonr_model_obs(
+        dummy_cubelist_model_obs, ["longitude", "latitude"]
+    )
+    assert isinstance(corr, CubeList)
+    assert len(corr) == 2
+    model_names = ["model_a", "model_b"]
+    for cube, model_name in zip(corr, model_names, strict=True):
+        assert (
+            cube.name() == "Pearson_Correlation_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_pearson_correlation_preserve_in_timelatlon(dummy_cubelist_model_obs):
+    """Pearson correlation with nothing collapsed, preserving both time and station dimensions."""
+    with pytest.raises(
+        ValueError, match="You cannot preserve all dimensions with pearsonr."
+    ):
+        scoreswrappers.scores_correlation_pearsonr_model_obs(
+            dummy_cubelist_model_obs, ["time", "longitude", "latitude"]
+        )
+
+
 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(