From 0f62a74224bbd97d36cae9f5bedb452a2f9d60a3 Mon Sep 17 00:00:00 2001 From: James Ross Date: Mon, 17 Aug 2026 11:56:37 +0100 Subject: [PATCH] adding model obs MAE --- src/CSET/operators/scoreswrappers.py | 61 +++++++++++++- .../surface_difference_scores_MAE.yaml | 2 +- .../surface_scores_model_vs_obs_MAE.yaml | 82 +++++++++++++++++++ ...eries_surface_scores_model_vs_obs_MAE.yaml | 65 +++++++++++++++ tests/operators/test_scoreswrappers.py | 43 ++++++++++ 5 files changed, 250 insertions(+), 3 deletions(-) create mode 100644 src/CSET/recipes/verification/surface_scores_model_vs_obs_MAE.yaml create mode 100644 src/CSET/recipes/verification/timeseries_surface_scores_model_vs_obs_MAE.yaml diff --git a/src/CSET/operators/scoreswrappers.py b/src/CSET/operators/scoreswrappers.py index 681f7f5be..6e09485a6 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_mae_model_obs( + cubes: CubeList, preserved_coordinates: list[str] | str | None = None +): + r"""Calculate the Mean Absolute Error (MAE) using scores. + + Acts as a wrapper around the MAE 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 + MAE. 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 MAE between the models and observation cube. + """ + mae_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) + mae = scores_mae( + input_cubelist, preserved_coordinates, obs_model_comparison=True + ) + model_name = model.attributes["model_name"] + mae.attributes["model_name"] = model_name + mae_cubes.append(mae) + + return mae_cubes + + def scores_rmse( cubes: CubeList, preserved_coordinates: list[str] | str | None = None, @@ -332,7 +378,11 @@ def scores_rmse( return scores_cube -def scores_mae(cubes: CubeList, preserved_coordinates: list[str] | str | None = None): +def scores_mae( + cubes: CubeList, + preserved_coordinates: list[str] | str | None = None, + obs_model_comparison: bool = False, +): r"""Calculate the Mean Absolute Error (MAE) using scores. Acts as a wrapper around the MAE calculation from ``scores`` ([scoresa]_, [scoresb]_). @@ -354,7 +404,14 @@ def scores_mae(cubes: CubeList, preserved_coordinates: list[str] | str | None = scores_cube: iris.cube.Cube A cube containing the MAE 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_difference_scores_MAE.yaml b/src/CSET/recipes/verification/surface_difference_scores_MAE.yaml index cc77db011..e202a9acf 100644 --- a/src/CSET/recipes/verification/surface_difference_scores_MAE.yaml +++ b/src/CSET/recipes/verification/surface_difference_scores_MAE.yaml @@ -16,7 +16,7 @@ description: | This will produce a single output plot of the total MAE. A larger MAE implies a greater error than a smaller MAE. An - RMSE of zero indicates the two fields most likely match. The MAE is calculated + MAE of zero indicates the two fields most likely match. The MAE is calculated on the grid point and thus a spatial view of the MAE provides useful information about where the differences are, or if placement errors are domininating the score (usually indicated by dipoles). MAE is a fair diff --git a/src/CSET/recipes/verification/surface_scores_model_vs_obs_MAE.yaml b/src/CSET/recipes/verification/surface_scores_model_vs_obs_MAE.yaml new file mode 100644 index 000000000..ae729591b --- /dev/null +++ b/src/CSET/recipes/verification/surface_scores_model_vs_obs_MAE.yaml @@ -0,0 +1,82 @@ +category: Surface Model vs Observation +title: "$VARNAME scores mae spatial plot for model vs observation points $SUBAREA_NAME" +description: | + Extracts and plot a spatial plot of the mean absolute error of $VARNAME for all + times based on mean of matched model and observation points. + + The MAE is calculated based on that used in the + package [`scores`](https://scores.readthedocs.io/en/stable/api.html#scores.continuous.mae). + This recipe preserves the time, latitude and longitude coordinates. + This means a sequence of spatial plots of the MAE can be produced. + + This recipe can preserve the time, latitude and longitude coordinates, in which case + a sequence of spatial plots of the MAE are produced according to the averaging + method i.e. SEQ, MAX and MIN. + + Alternatively, the MAE 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 MAE. + + A larger MAE implies a greater error than a smaller MAE. An + MAE of zero indicates the two fields most likely match. The MAE is calculated + on the grid point and thus a spatial view of the MAE provides useful + information about where the differences are, or if placement errors + are domininating the score (usually indicated by dipoles). MAE is a fair + measure of error when forecasting the median and is not sensitive to outliers. + + For MAE values computed across an entire casse study, the computation + is done within the scores MAE module, rather than trying to mean across + each timestep. This preserves the nonlinear nature of the MAE 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_mae_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_MAE.yaml b/src/CSET/recipes/verification/timeseries_surface_scores_model_vs_obs_MAE.yaml new file mode 100644 index 000000000..64f362298 --- /dev/null +++ b/src/CSET/recipes/verification/timeseries_surface_scores_model_vs_obs_MAE.yaml @@ -0,0 +1,65 @@ +category: Surface Model verus Observation Time Series +title: "$VARNAME scores mae time series for model vs observation points $SUBAREA_NAME" +description: Extracts and plot a time series of $VARNAME for all times based on RMSE between given models 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.mae). + + This recipe allows the preservation of the time coordinate to produce a timeseries. Therefore, the MAE is + collapsed over all other coordinates in the cube and calculated for every timestep. + + + A larger MAE implies a greater error than a smaller MAE. An + MAE of zero indicates the two fields most likely match. The MAE is calculated + on the grid point and thus a spatial view of the MAE provides useful + information about where the differences are, or if placement errors + are domininating the score (usually indicated by dipoles). MAE is a fair + measure of error when forecasting the median and is not sensitive to outliers. + + +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_mae_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..616848d40 100644 --- a/tests/operators/test_scoreswrappers.py +++ b/tests/operators/test_scoreswrappers.py @@ -344,6 +344,49 @@ def test_model_obs_rmse_preserve_in_timelatlon(dummy_cubelist_model_obs): assert cube.attributes["model_name"] == model_name +def test_model_obs_mae_preserve_in_time(dummy_cubelist_model_obs): + """MAE collapsed over station, preserving only the time dimension.""" + mae = scoreswrappers.scores_mae_model_obs(dummy_cubelist_model_obs, "time") + assert isinstance(mae, CubeList) + assert len(mae) == 2 + model_names = ["model_a", "model_b"] + for cube, model_name in zip(mae, model_names, strict=True): + assert cube.name() == "MAE_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_mae_preserve_in_latlon(dummy_cubelist_model_obs): + """MAE collapsed over time, preserving the station dimension via lat/lon.""" + mae = scoreswrappers.scores_mae_model_obs( + dummy_cubelist_model_obs, ["longitude", "latitude"] + ) + assert isinstance(mae, CubeList) + assert len(mae) == 2 + model_names = ["model_a", "model_b"] + for cube, model_name in zip(mae, model_names, strict=True): + assert cube.name() == "MAE_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_mae_preserve_in_timelatlon(dummy_cubelist_model_obs): + """MAE with nothing collapsed, preserving both time and station dimensions.""" + mae = scoreswrappers.scores_mae_model_obs( + dummy_cubelist_model_obs, ["time", "longitude", "latitude"] + ) + assert isinstance(mae, CubeList) + assert len(mae) == 2 + model_names = ["model_a", "model_b"] + for cube, model_name in zip(mae, model_names, strict=True): + assert cube.name() == "MAE_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(