diff --git a/src/CSET/operators/scoreswrappers.py b/src/CSET/operators/scoreswrappers.py index 681f7f5be..3383ade79 100644 --- a/src/CSET/operators/scoreswrappers.py +++ b/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) diff --git a/src/CSET/recipes/verification/surface_scores_model_vs_obs_additive_bias.yaml b/src/CSET/recipes/verification/surface_scores_model_vs_obs_additive_bias.yaml new file mode 100644 index 000000000..e10d7dad7 --- /dev/null +++ b/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 diff --git a/src/CSET/recipes/verification/timeseries_surface_scores_model_vs_obs_additive_bias.yaml b/src/CSET/recipes/verification/timeseries_surface_scores_model_vs_obs_additive_bias.yaml new file mode 100644 index 000000000..2bf8306c3 --- /dev/null +++ b/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 diff --git a/tests/operators/test_scoreswrappers.py b/tests/operators/test_scoreswrappers.py index a4430a27a..926e00279 100644 --- a/tests/operators/test_scoreswrappers.py +++ b/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(