Add Pearson correlation coefficient into verification

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

@@ -73,3 +73,17 @@ description=Create a timeseries of the RMSE for the specified surface fields.
 type=python_boolean
 compulsory=true
 sort-key=scores2
+
+[template variables=SCORES_PEARSONR_SPATIAL]
+ns=Verification/Scores
+description=Create spatially mapped Pearson correlation coeff plots for the specified surface fields.
+type=python_boolean
+compulsory=true
+sort-key=scores3
+
+[template variables=SCORES_PEARSONR_TIMESERIES]
+ns=Verification/Scores
+description=Create a timeseries of the Pearson correlation coeff for the specified surface fields.
+type=python_boolean
+compulsory=true
+sort-key=scores4

src/CSET/loaders/verification.py

@@ -63,3 +63,42 @@ def load(conf: Config):
                 model_ids=[base_model["id"], model["id"]],
                 aggregation=False,
             )
+
+    if conf.SCORES_PEARSONR_SPATIAL:
+        base_model = models[0]
+        for model, field, method in itertools.product(
+            models[1:], conf.SURFACE_FIELDS, conf.SPATIAL_SURFACE_FIELD_METHOD
+        ):
+            yield RawRecipe(
+                recipe="surface_pearsonr_scores.yaml",
+                variables={
+                    "VARNAME": field,
+                    "BASE_MODEL": base_model["name"],
+                    "OTHER_MODEL": model["name"],
+                    "METHOD": method,
+                    "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None,
+                    "SUBAREA_EXTENT": conf.SUBAREA_EXTENT
+                    if conf.SELECT_SUBAREA
+                    else None,
+                },
+                model_ids=[base_model["id"], model["id"]],
+                aggregation=False,
+            )
+
+    if conf.SCORES_PEARSONR_TIMESERIES:
+        base_model = models[0]
+        for model, field in itertools.product(models[1:], conf.SURFACE_FIELDS):
+            yield RawRecipe(
+                recipe="timeseries_surface_pearsonr_scores.yaml",
+                variables={
+                    "VARNAME": field,
+                    "BASE_MODEL": base_model["name"],
+                    "OTHER_MODEL": model["name"],
+                    "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None,
+                    "SUBAREA_EXTENT": conf.SUBAREA_EXTENT
+                    if conf.SELECT_SUBAREA
+                    else None,
+                },
+                model_ids=[base_model["id"], model["id"]],
+                aggregation=False,
+            )

src/CSET/operators/scoreswrappers.py

@@ -192,3 +192,46 @@ def scores_rmse(cubes: CubeList, preserved_coordinates: list[str] | str | None =
     )
     RMSE.rename(f"RMSE_of_{base.name()}")
     return RMSE
+
+
+def scores_pearsonr(
+    cubes: CubeList, preserved_coordinates: list[str] | str | None = None
+):
+    """Calculate Pearson's correlation coefficient using scores.
+
+    Acts as a wrapper around the pearsonr 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 calculaiton of the
+        RMSE. For example if you want a map of each time you can preserve
+        ["time","grid_latitude", "grid_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
+    -------
+    pearsonr: iris.cube.Cube
+        A cube containing the Pearson correlation coefficient between the base and other cube.
+
+    References
+    ----------
+    add references here
+    """
+    base, other = _sort_cubes_for_verification(cubes)
+    # Scores operators on xarray data arrays, so we transform the iris cube into an array,
+    # apply scores, and then transform it back.
+
+    pearsonr = xr.DataArray.to_iris(
+        scores.continuous.correlation.pearsonr(
+            xr.DataArray.from_iris(other),
+            xr.DataArray.from_iris(base),
+            preserve_dims=preserved_coordinates,
+        )
+    )
+    pearsonr.rename(f"Pearsonr_of_{base.name()}")
+    return pearsonr

src/CSET/recipes/verification/surface_pearsonr_scores.yaml

@@ -0,0 +1,41 @@
+category: Scores
+title: "$METHOD $VARNAME\nRMSE between $OTHER_MODEL and $BASE_MODEL"
+description: |
+  Extracts and plots the Pearson Correlation Coefficient in $METHOD of $VARNAME
+  for all times. The Pearson Correlation Coefficient is calculated based on that used in the
+  package [`scores`](https://scores.readthedocs.io/en/stable/api.html#scores.continuous.correlation.pearsonr). This recipe preserves the time, latitude and longitude coordinates.
+  This means a sequence of spatial plots of the Pearson Correlation Coefficient can be produced.
+
+  Description here.....
+
+steps:
+  - operator: read.read_cubes
+    file_paths: $INPUT_PATHS
+    model_names: [$BASE_MODEL, $OTHER_MODEL]
+    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
+      pressure_level_constraint:
+        operator: constraints.generate_level_constraint
+        coordinate: "pressure"
+        levels: []
+    subarea_type: $SUBAREA_TYPE
+    subarea_extent: $SUBAREA_EXTENT
+
+  - operator: scoreswrappers.scores_pearsonr
+    preserved_coordinates: ["time", "grid_latitude", "grid_longitude"]
+
+  - 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_pearsonr_scores.yaml

@@ -0,0 +1,40 @@
+category: Scores
+title: "$VARNAME\nPEARSONR timeseries between $OTHER_MODEL and $BASE_MODEL"
+description: |
+  Extracts and plots the Pearsons correlation coefficient in $VARNAME
+  for all times. The Pearsons correlation coefficient 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
+  Pearsons correlation coefficient is collapsed over all other coordinates in the cube
+  and calculated for every timestep.
+
+  More description here
+
+steps:
+  - operator: read.read_cubes
+    file_paths: $INPUT_PATHS
+    model_names: [$BASE_MODEL, $OTHER_MODEL]
+    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
+      pressure_level_constraint:
+        operator: constraints.generate_level_constraint
+        coordinate: "pressure"
+        levels: []
+    subarea_type: $SUBAREA_TYPE
+    subarea_extent: $SUBAREA_EXTENT
+
+  - operator: scoreswrappers.scores_pearsonr
+    preserved_coordinates: ["time"]
+
+  - operator: plot.plot_line_series
+    series_coordinate: time
+
+  - operator: write.write_cube_to_nc
+    overwrite: True