From dc4e5ef20be1f56c2ff3788f06bac731b138bd9d Mon Sep 17 00:00:00 2001 From: James Ross Date: Mon, 17 Aug 2026 13:32:07 +0100 Subject: [PATCH 1/3] adding pearson corr # Conflicts: # src/CSET/operators/scoreswrappers.py # tests/operators/test_scoreswrappers.py --- src/CSET/operators/scoreswrappers.py | 61 +++++++++++++- ...ores_model_vs_obs_correlation_pearson.yaml | 79 +++++++++++++++++++ tests/operators/test_scoreswrappers.py | 44 +++++++++++ 3 files changed, 182 insertions(+), 2 deletions(-) create mode 100644 src/CSET/recipes/verification/surface_scores_model_vs_obs_correlation_pearson.yaml diff --git a/src/CSET/operators/scoreswrappers.py b/src/CSET/operators/scoreswrappers.py index ff55ba9b3..47bd3133f 100644 --- a/src/CSET/operators/scoreswrappers.py +++ b/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.Cube + A cube 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) diff --git a/src/CSET/recipes/verification/surface_scores_model_vs_obs_correlation_pearson.yaml b/src/CSET/recipes/verification/surface_scores_model_vs_obs_correlation_pearson.yaml new file mode 100644 index 000000000..28b8a5aec --- /dev/null +++ b/src/CSET/recipes/verification/surface_scores_model_vs_obs_correlation_pearson.yaml @@ -0,0 +1,79 @@ +category: Surface Model vs Observation +title: "$VARNAME scores rmse spatial plot for model vs observation points $SUBAREA_NAME" +description: | + Extracts and plot a spatial plot of the RMSE of $VARNAME for all + times based on mean of matched model and observation points. + + The RMSE is calculated based on that used in the + package [`scores`](https://scores.readthedocs.io/en/stable/api.html#scores.continuous.rmse). + + This recipe can preserve the time, latitude and longitude coordinates, in which case + a sequence of spatial plots of the RMSE are produced according to the averaging + method i.e. SEQ, MAX and MIN. + + Alternatively, the RMSE 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 RMSE. + + A larger RMSE implies a greater error than a smaller RMSE. An + RMSE of zero indicates the two fields match. The RMSE is calculated + on the grid point and thus a spatial view of the RMSE provides useful + information about where the differences are, or if placement errors + are domininating the score (usually indicated by dipoles). + + For RMSE values computed across an entire casse study, the computation + is done within the scores RMSE module, rather than trying to mean across + each timestep. This preserves the nonlinear nature of the RMSE metric. + + + +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 diff --git a/tests/operators/test_scoreswrappers.py b/tests/operators/test_scoreswrappers.py index 64e642c17..ca8517f60 100644 --- a/tests/operators/test_scoreswrappers.py +++ b/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( From dbd01b2fc6d993bcc202f1478c4268511e9bad6d Mon Sep 17 00:00:00 2001 From: James Ross Date: Mon, 17 Aug 2026 14:21:20 +0100 Subject: [PATCH 2/3] cube -> cubelist --- src/CSET/operators/scoreswrappers.py | 4 +- ...ores_model_vs_obs_correlation_pearson.yaml | 79 ------------------- ...res_model_vs_obs_correlation_pearsonr.yaml | 60 ++++++++++++++ tests/operators/test_scoreswrappers.py | 45 ----------- 4 files changed, 62 insertions(+), 126 deletions(-) delete mode 100644 src/CSET/recipes/verification/surface_scores_model_vs_obs_correlation_pearson.yaml create mode 100644 src/CSET/recipes/verification/timeseries_surface_scores_model_vs_obs_correlation_pearsonr.yaml diff --git a/src/CSET/operators/scoreswrappers.py b/src/CSET/operators/scoreswrappers.py index 47bd3133f..1dfad30c2 100644 --- a/src/CSET/operators/scoreswrappers.py +++ b/src/CSET/operators/scoreswrappers.py @@ -347,8 +347,8 @@ def scores_correlation_pearsonr_model_obs( Returns ------- - scores_cube: iris.cube.Cube - A cube containing the PC between the models and observation cube. + scores_cube: iris.cube.CubeList + A cube list containing the PC between the models and observation cube. """ pearsonr_cubes = CubeList() model_list = CubeList() diff --git a/src/CSET/recipes/verification/surface_scores_model_vs_obs_correlation_pearson.yaml b/src/CSET/recipes/verification/surface_scores_model_vs_obs_correlation_pearson.yaml deleted file mode 100644 index 28b8a5aec..000000000 --- a/src/CSET/recipes/verification/surface_scores_model_vs_obs_correlation_pearson.yaml +++ /dev/null @@ -1,79 +0,0 @@ -category: Surface Model vs Observation -title: "$VARNAME scores rmse spatial plot for model vs observation points $SUBAREA_NAME" -description: | - Extracts and plot a spatial plot of the RMSE of $VARNAME for all - times based on mean of matched model and observation points. - - The RMSE is calculated based on that used in the - package [`scores`](https://scores.readthedocs.io/en/stable/api.html#scores.continuous.rmse). - - This recipe can preserve the time, latitude and longitude coordinates, in which case - a sequence of spatial plots of the RMSE are produced according to the averaging - method i.e. SEQ, MAX and MIN. - - Alternatively, the RMSE 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 RMSE. - - A larger RMSE implies a greater error than a smaller RMSE. An - RMSE of zero indicates the two fields match. The RMSE is calculated - on the grid point and thus a spatial view of the RMSE provides useful - information about where the differences are, or if placement errors - are domininating the score (usually indicated by dipoles). - - For RMSE values computed across an entire casse study, the computation - is done within the scores RMSE module, rather than trying to mean across - each timestep. This preserves the nonlinear nature of the RMSE metric. - - - -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 diff --git a/src/CSET/recipes/verification/timeseries_surface_scores_model_vs_obs_correlation_pearsonr.yaml b/src/CSET/recipes/verification/timeseries_surface_scores_model_vs_obs_correlation_pearsonr.yaml new file mode 100644 index 000000000..98c4de369 --- /dev/null +++ b/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 diff --git a/tests/operators/test_scoreswrappers.py b/tests/operators/test_scoreswrappers.py index ca8517f60..cb0a750cd 100644 --- a/tests/operators/test_scoreswrappers.py +++ b/tests/operators/test_scoreswrappers.py @@ -395,51 +395,6 @@ def test_model_obs_pearson_correlation_preserve_in_timelatlon(dummy_cubelist_mod ) -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( From 0c9acdaa7165b3cf98ab2b21cdd697ffacb7b72a Mon Sep 17 00:00:00 2001 From: James Ross Date: Fri, 21 Aug 2026 10:24:45 +0100 Subject: [PATCH 3/3] add tests that i accidently removed --- tests/operators/test_scoreswrappers.py | 45 ++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/tests/operators/test_scoreswrappers.py b/tests/operators/test_scoreswrappers.py index cb0a750cd..ca8517f60 100644 --- a/tests/operators/test_scoreswrappers.py +++ b/tests/operators/test_scoreswrappers.py @@ -395,6 +395,51 @@ def test_model_obs_pearson_correlation_preserve_in_timelatlon(dummy_cubelist_mod ) +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(