diff --git a/src/CSET/operators/scoreswrappers.py b/src/CSET/operators/scoreswrappers.py index ff55ba9b3..1dfad30c2 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.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) 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 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(