diff --git a/src/CSET/cset_workflow/meta/verification/rose-meta.conf b/src/CSET/cset_workflow/meta/verification/rose-meta.conf index bdd307a73..5565dd694 100644 --- a/src/CSET/cset_workflow/meta/verification/rose-meta.conf +++ b/src/CSET/cset_workflow/meta/verification/rose-meta.conf @@ -53,38 +53,140 @@ sort-key=met5 ################################### -# Scores -[Verification/Scores] -ns=Verification/Scores -sort-key=scores -title=Scores [template variables=SCORES_ALL] -ns=Verification/Scores +ns=Verification description=If selected, all scores shown below will be calculated. NOTE scores currently only available for model-model intercomparisons. type=python_boolean -trigger=template variables=SCORES_RMSE_SPATIAL: False; - template variables=SCORES_RMSE_TIMESERIES: False; +trigger=template variables=SCORES_SPATIAL_RMSE: False; + template variables=SCORES_SPATIAL_MAE: False; + template variables=SCORES_SPATIAL_AB: False; + template variables=SCORES_TIMESERIES_RMSE: False; + template variables=SCORES_TIMESERIES_MAE: False; + template variables=SCORES_TIMESERIES_AB: False; + template variables=SCORES_TIMESERIES_PC: False; compulsory=true sort-key=scores1 -[template variables=SCORES_RMSE_SPATIAL] -ns=Verification/Scores -description=Create spatially mapped RMSE plots for the specified surface fields. - The analysis method(s) set in SPATIAL_SURFACE_FIELD_METHOD will be used. - NOTE currently only available for model-model intercomparisons. + +# Scores spatial plots +[Verification/ScoresSpatial] +ns=Verification/ScoresSpatial +sort-key=scoresspatial +title=Scores spatial + +[template variables=SPATIAL_SCORES_FIELD_METHOD] +ns=Verification/ScoresSpatial +description=Select analysis method(s) for output mapped plots. Add all options required. + Leave blank or set to "SEQ" (sequence) for plots each diagnostic output time. + For time-collapsed outputs over each analysis period, set to "CASE", "MAX", "MIN" etc. + Selecting "CASE" will compute the spatial scores metrics over the entire case study. +help=Quoted analysis methods. Settings should be based on available iris.analysis methods for collapsing cube dimensions. +type=python_list +compulsory=true +sort-key=sc-sp-0 + +[template variables=SCORES_SPATIAL_RMSE] +ns=Verification/ScoresSpatial +title=Root Mean Square Error +description=SCORES COUNTINOUS METRIC. + Plots a 2D spatial plot of the RMSE between two models. + The computation is implemented by the scores package. + For the time coordinate, the method(s) set in + SPATIAL_SCORES_FIELD_METHOD will be used. +help=This computes the metric using the scores package. Details of this can + be found at https://scores.readthedocs.io/en/stable/included.html#continuous type=python_boolean compulsory=true -sort-key=scores2 +sort-key=sc-sp-1 -[template variables=SCORES_RMSE_TIMESERIES] -ns=Verification/Scores -description=Create a timeseries of the RMSE for the specified surface fields. - NOTE currently only available for model-model intercomparisons. +[template variables=SCORES_SPATIAL_MAE] +ns=Verification/ScoresSpatial +title=Mean Absolute Error +description=SCORES COUNTINOUS METRIC. + Plots a 2D spatial plot of the MAE between two models. + The computation is implemented by the scores package. + For the time coordinate, the method(s) set in + SPATIAL_SCORES_FIELD_METHOD will be used. +help=This computes the metric using the scores package. Details of this can + be found at https://scores.readthedocs.io/en/stable/included.html#continuous +type=python_boolean +compulsory=true +sort-key=sc-sp-2 + +[template variables=SCORES_SPATIAL_AB] +ns=Verification/ScoresSpatial +title=Additive Bias (Mean Error) +description=SCORES COUNTINOUS METRIC. + Plots a 2D spatial plot of the Additive Bias between two models. + The Additive Bias is also known as the Mean Error. + The computation is implemented by the scores package. + For the time coordinate, the method(s) set in + SPATIAL_SCORES_FIELD_METHOD will be used. +help=This computes the metric using the scores package. Details of this can + be found at https://scores.readthedocs.io/en/stable/included.html#continuous +type=python_boolean +compulsory=true +sort-key=sc-sp-3 + + +################################### + +# Scores timeseries plots +[Verification/ScoresTimeseries] +ns=Verification/ScoresTimeseries +sort-key=scorestimeseries +title=Scores timeseries + +[template variables=SCORES_TIMESERIES_RMSE] +ns=Verification/ScoresTimeseries +title=Root Mean Square Error +description=SCORES COUNTINOUS METRIC. + Plots a timeseries plot of the RMSE between two models. + The computation is implemented by the scores package. +help=This computes the metric using the scores package. Details of this can + be found at https://scores.readthedocs.io/en/stable/included.html#continuous +type=python_boolean +compulsory=true +sort-key=sc-ts-1 + +[template variables=SCORES_TIMESERIES_MAE] +ns=Verification/ScoresTimeseries +title=Mean Absolute Error +description=SCORES COUNTINOUS METRIC. + Plots a timeseries plot of the MAE between two models. + The computation is implemented by the scores package. +help=This computes the metric using the scores package. Details of this can + be found at https://scores.readthedocs.io/en/stable/included.html#continuous +type=python_boolean +compulsory=true +sort-key=sc-ts-2 + +[template variables=SCORES_TIMESERIES_AB] +ns=Verification/ScoresTimeseries +title=Additive Bias (Mean Error) +description=SCORES COUNTINOUS METRIC. + Plots a timeseries plot of the Additive Bias between two models. + The Additive Bias is also known as the Mean Error. + The computation is implemented by the scores package. +help=This computes the metric using the scores package. Details of this can + be found at https://scores.readthedocs.io/en/stable/included.html#continuous +type=python_boolean +compulsory=true +sort-key=sc-ts-3 + +[template variables=SCORES_TIMESERIES_PC] +ns=Verification/ScoresTimeseries +title=Pearson's Correlation +description=SCORES COUNTINOUS METRIC. + Plots a timeseries plot of the Pearson’s Correlation Coefficient between two models. + The computation is implemented by the scores package. +help=This computes the metric using the scores package. Details of this can + be found at https://scores.readthedocs.io/en/stable/included.html#continuous type=python_boolean compulsory=true -sort-key=scores3 +sort-key=sc-ts-4 [template variables=SCORES_CRPS_FOR_ENSEMBLE] ns=Verification/Scores diff --git a/src/CSET/loaders/verification.py b/src/CSET/loaders/verification.py index 6d1c05a8a..18623c243 100644 --- a/src/CSET/loaders/verification.py +++ b/src/CSET/loaders/verification.py @@ -19,24 +19,76 @@ from CSET.recipes import Config, RawRecipe, get_models +def _get_scores_spatial_methods(conf): + """Compile list of the required scores spatial plots.""" + scores_spatial_methods = [] + if conf.SCORES_SPATIAL_RMSE or conf.SCORES_ALL: + scores_spatial_methods.append("RMSE") + if conf.SCORES_SPATIAL_AB or conf.SCORES_ALL: + scores_spatial_methods.append("additive_bias") + if conf.SCORES_SPATIAL_MAE or conf.SCORES_ALL: + scores_spatial_methods.append("MAE") + return scores_spatial_methods + + +def _get_scores_timeseries_methods(conf): + """Compile list of the required scores timeseries plots.""" + scores_timeseries_methods = [] + if conf.SCORES_TIMESERIES_RMSE or conf.SCORES_ALL: + scores_timeseries_methods.append("RMSE") + if conf.SCORES_TIMESERIES_AB or conf.SCORES_ALL: + scores_timeseries_methods.append("additive_bias") + if conf.SCORES_TIMESERIES_MAE or conf.SCORES_ALL: + scores_timeseries_methods.append("MAE") + if conf.SCORES_TIMESERIES_PC or conf.SCORES_ALL: + scores_timeseries_methods.append("correlation_pearsonr") + return scores_timeseries_methods + + def load(conf: Config): """Yield recipes from the given workflow configuration.""" # Load a list of model detail dictionaries. models = get_models(conf.asdict()) # Models are listed in order, so model 1 is the first element. - if conf.SCORES_RMSE_SPATIAL or conf.SCORES_ALL: + scores_spatial_methods = _get_scores_spatial_methods(conf) + if scores_spatial_methods: + # Produce 2D spatial plots of scores metrics. base_model = models[0] - for model, field, method in itertools.product( - models[1:], conf.SURFACE_FIELDS, conf.SPATIAL_SURFACE_FIELD_METHOD + for model, field, method, scores_method in itertools.product( + models[1:], + conf.SURFACE_FIELDS, + conf.SPATIAL_SCORES_FIELD_METHOD, + scores_spatial_methods, ): + preserved_coords = ["time", "grid_latitude", "grid_longitude"] + method_null = "" + scores_method_case = "CASE" + scores_coords_case = ["grid_latitude", "grid_longitude"] + if scores_method == "RMSE" and method == scores_method_case: + # Set the preserved coords and collapse method required + # to produce RMSE spatial plot over an entire case study. + preserved_coords = scores_coords_case + method = method_null + if scores_method == "MAE" and method == scores_method_case: + # Set the preserved coords and collapse method required + # to produce MAE spatial plot over an entire case study. + preserved_coords = scores_coords_case + method = method_null + if scores_method == "additive_bias" and method == scores_method_case: + # Set the preserved coords and collapse method required + # to produce ME additive bias spatial plot over an entire case study. + preserved_coords = scores_coords_case + method = method_null yield RawRecipe( - recipe="surface_rmse_scores.yaml", + recipe=f"surface_difference_scores_{scores_method}.yaml", variables={ "VARNAME": field, "BASE_MODEL": base_model["name"], "OTHER_MODEL": model["name"], "METHOD": method, + "PRESERVED_COORDS": preserved_coords, + "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA @@ -46,15 +98,20 @@ def load(conf: Config): aggregation=False, ) - if conf.SCORES_RMSE_TIMESERIES or conf.SCORES_ALL: + scores_timeseries_methods = _get_scores_timeseries_methods(conf) + if scores_timeseries_methods: + # Produce timeseries plots of scores metrics averaged over the domain for each case study. base_model = models[0] - for model, field in itertools.product(models[1:], conf.SURFACE_FIELDS): + for model, field, scores_method in itertools.product( + models[1:], conf.SURFACE_FIELDS, scores_timeseries_methods + ): yield RawRecipe( - recipe="timeseries_surface_rmse_scores.yaml", + recipe=f"timeseries_surface_difference_scores_{scores_method}.yaml", variables={ "VARNAME": field, "BASE_MODEL": base_model["name"], "OTHER_MODEL": model["name"], + "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA diff --git a/src/CSET/operators/scoreswrappers.py b/src/CSET/operators/scoreswrappers.py index 26e6cc449..dad8ddc61 100644 --- a/src/CSET/operators/scoreswrappers.py +++ b/src/CSET/operators/scoreswrappers.py @@ -146,6 +146,53 @@ def _sort_cubes_for_verification(cubes: CubeList): return base, other +def _resolve_preserve_dims( + cube: Cube, + data_array: xr.DataArray, + preserved_coordinates: list[str] | str | None, +) -> list[str] | None: + """Resolve preserve coordinates to xarray dimension names. + + The ``scores`` package expects preserve dimensions to match xarray + dimension names. In Iris data, commonly used coordinates such as ``time`` + may be auxiliary coordinates attached to a differently named dimension + (e.g. ``dim0``). This helper maps coordinate names to their underlying + dimension names and helps to convert from iris to xarray coordinate dimension names. + """ + if preserved_coordinates is None: + return None + + coord_names = ( + [preserved_coordinates] + if isinstance(preserved_coordinates, str) + else preserved_coordinates + ) + preserve_dims: list[str] = [] + + for coord_name in coord_names: + # Already an xarray dimension name. + if coord_name in data_array.dims: + if coord_name not in preserve_dims: + preserve_dims.append(coord_name) + continue + + # Otherwise, map coordinate name to dimension index/indices. + try: + dim_indices = cube.coord_dims(coord_name) + except iris.exceptions.CoordinateNotFoundError: + # Keep original name so scores raises a clear error for unknown keys. + if coord_name not in preserve_dims: + preserve_dims.append(coord_name) + continue + + for dim_index in dim_indices: + dim_name = data_array.dims[dim_index] + if dim_name not in preserve_dims: + preserve_dims.append(dim_name) + + return preserve_dims + + def scores_rmse(cubes: CubeList, preserved_coordinates: list[str] | str | None = None): r"""Calculate the Root Mean Square Error (RMSE) using scores. @@ -168,7 +215,7 @@ def scores_rmse(cubes: CubeList, preserved_coordinates: list[str] | str | None = Returns ------- - RMSE: iris.cube.Cube + scores_cube: iris.cube.Cube A cube containing the RMSE between the base and other cube. References @@ -187,17 +234,315 @@ def scores_rmse(cubes: CubeList, preserved_coordinates: list[str] | str | None = forecasts, predictions or models (2.5.0)". Zenodo. doi: 10.5281/zenodo.18638494 """ base, other = _sort_cubes_for_verification(cubes) - # Scores operators on xarray data arrays, so we transform the iris cube into an array, + + # Copy the coordinates of the input cubes. + other_xr = xr.DataArray.from_iris(other) + base_xr = xr.DataArray.from_iris(base) + preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates) + + # Scores operates on xarray data arrays, so we transform the iris cube into an array, # apply scores, and then transform it back. - RMSE = xr.DataArray.to_iris( + scores_cube = xr.DataArray.to_iris( scores.continuous.rmse( - xr.DataArray.from_iris(other), - xr.DataArray.from_iris(base), - preserve_dims=preserved_coordinates, + other_xr, + base_xr, + preserve_dims=preserve_dims, ) ) - RMSE.rename(f"RMSE_of_{base.name()}") - return RMSE + + # If time is aggregated out, attach a scalar time coordinate with bounds + # so plotting can display the aggregated period in the title. + try: + if not scores_cube.coords("time"): + base_time = base.coord("time") + time_vals = ( + base_time.bounds.flatten() + if base_time.has_bounds() + else base_time.points + ) + t_start = float(time_vals[0]) + t_end = float(time_vals[-1]) + t_mid = 0.5 * (t_start + t_end) + + scores_cube.add_aux_coord( + iris.coords.AuxCoord( + t_mid, + standard_name=base_time.standard_name, + long_name=base_time.long_name, + var_name=base_time.var_name, + units=base_time.units, + bounds=np.array([t_start, t_end]), + attributes=base_time.attributes.copy(), + ) + ) + except iris.exceptions.CoordinateNotFoundError: + pass + + scores_cube.rename(f"RMSE_of_{base.name()}") + # if preserved_coordinates == ["grid_latitude", "grid_longitude"]: + # scores_cube.add_aux_coord(time_coord) + return scores_cube + + +def scores_mae(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 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 + MAE. 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 + ------- + scores_cube: iris.cube.Cube + A cube containing the MAE between the base and other cube. + + References + ---------- + .. [scoresa] Leeuwenburg, T., Loveday, N., Ebert, E. E., Cook, H., + Khanarmuei, M., Taggart, R. J., Ramanathan, N., Carroll, M., Chong, S., + Griffiths, A., & Sharples, J. (2024) "scores: A Python package for + verifying and evaluating models and predictions with xarray". Journal + of Open Source Software, vol. 9, 6889. doi: 10.21105/joss.06889 + + .. [scoresb] Leeuwenburg, T., Loveday, N., Ramanathan, N., Chong, S., + Taggart, R. J., Shrestha, D., Khanarmuei, M., Cook, H., Bluett, L., Ebert, + E. E., Carroll, M., Trotta, B., Bishop, S., Squire, D. T., Griffiths, A., + Pagano, T. C., Fisher, A. J., Mandelbaum, T., Jinghan, F., … Smallwood, J. + (2026) "scores: Metrics for the verification, evaluation and optimisation of + forecasts, predictions or models (2.5.0)". Zenodo. doi: 10.5281/zenodo.18638494 + """ + base, other = _sort_cubes_for_verification(cubes) + + # Copy the coordinates of the input cubes. + other_xr = xr.DataArray.from_iris(other) + base_xr = xr.DataArray.from_iris(base) + preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates) + + # Scores operates on xarray data arrays, so we transform the iris cube into an array, + # apply scores, and then transform it back. + scores_cube = xr.DataArray.to_iris( + scores.continuous.mae( + other_xr, + base_xr, + preserve_dims=preserve_dims, + ) + ) + + # If time is aggregated out, attach a scalar time coordinate with bounds + # so plotting can display the aggregated period in the title. + try: + if not scores_cube.coords("time"): + base_time = base.coord("time") + time_vals = ( + base_time.bounds.flatten() + if base_time.has_bounds() + else base_time.points + ) + t_start = float(time_vals[0]) + t_end = float(time_vals[-1]) + t_mid = 0.5 * (t_start + t_end) + + scores_cube.add_aux_coord( + iris.coords.AuxCoord( + t_mid, + standard_name=base_time.standard_name, + long_name=base_time.long_name, + var_name=base_time.var_name, + units=base_time.units, + bounds=np.array([t_start, t_end]), + attributes=base_time.attributes.copy(), + ) + ) + except iris.exceptions.CoordinateNotFoundError: + pass + + scores_cube.rename(f"MAE_of_{base.name()}") + return scores_cube + + +def scores_additive_bias( + 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 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 + ME. 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 + ------- + scores_cube: iris.cube.Cube + A cube containing the ME between the base and other cube. + + References + ---------- + .. [scoresa] Leeuwenburg, T., Loveday, N., Ebert, E. E., Cook, H., + Khanarmuei, M., Taggart, R. J., Ramanathan, N., Carroll, M., Chong, S., + Griffiths, A., & Sharples, J. (2024) "scores: A Python package for + verifying and evaluating models and predictions with xarray". Journal + of Open Source Software, vol. 9, 6889. doi: 10.21105/joss.06889 + + .. [scoresb] Leeuwenburg, T., Loveday, N., Ramanathan, N., Chong, S., + Taggart, R. J., Shrestha, D., Khanarmuei, M., Cook, H., Bluett, L., Ebert, + E. E., Carroll, M., Trotta, B., Bishop, S., Squire, D. T., Griffiths, A., + Pagano, T. C., Fisher, A. J., Mandelbaum, T., Jinghan, F., … Smallwood, J. + (2026) "scores: Metrics for the verification, evaluation and optimisation of + forecasts, predictions or models (2.5.0)". Zenodo. doi: 10.5281/zenodo.18638494 + """ + base, other = _sort_cubes_for_verification(cubes) + + # Copy the coordinates of the input cubes. + other_xr = xr.DataArray.from_iris(other) + base_xr = xr.DataArray.from_iris(base) + preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates) + + # Scores operates on xarray data arrays, so we transform the iris cube into an array, + # apply scores, and then transform it back. + scores_cube = xr.DataArray.to_iris( + scores.continuous.additive_bias( + other_xr, + base_xr, + preserve_dims=preserve_dims, + ) + ) + + # If time is aggregated out, attach a scalar time coordinate with bounds + # so plotting can display the aggregated period in the title. + try: + if not scores_cube.coords("time"): + base_time = base.coord("time") + time_vals = ( + base_time.bounds.flatten() + if base_time.has_bounds() + else base_time.points + ) + t_start = float(time_vals[0]) + t_end = float(time_vals[-1]) + t_mid = 0.5 * (t_start + t_end) + + scores_cube.add_aux_coord( + iris.coords.AuxCoord( + t_mid, + standard_name=base_time.standard_name, + long_name=base_time.long_name, + var_name=base_time.var_name, + units=base_time.units, + bounds=np.array([t_start, t_end]), + attributes=base_time.attributes.copy(), + ) + ) + except iris.exceptions.CoordinateNotFoundError: + pass + scores_cube.rename(f"Additive_Bias_of_{base.name()}") + return scores_cube + + +def scores_correlation_pearsonr( + 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 calculaiton of the + PC. 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 + ------- + scores_cube: iris.cube.Cube + A cube containing the PC between the base and other cube. + + References + ---------- + .. [scoresa] Leeuwenburg, T., Loveday, N., Ebert, E. E., Cook, H., + Khanarmuei, M., Taggart, R. J., Ramanathan, N., Carroll, M., Chong, S., + Griffiths, A., & Sharples, J. (2024) "scores: A Python package for + verifying and evaluating models and predictions with xarray". Journal + of Open Source Software, vol. 9, 6889. doi: 10.21105/joss.06889 + + .. [scoresb] Leeuwenburg, T., Loveday, N., Ramanathan, N., Chong, S., + Taggart, R. J., Shrestha, D., Khanarmuei, M., Cook, H., Bluett, L., Ebert, + E. E., Carroll, M., Trotta, B., Bishop, S., Squire, D. T., Griffiths, A., + Pagano, T. C., Fisher, A. J., Mandelbaum, T., Jinghan, F., … Smallwood, J. + (2026) "scores: Metrics for the verification, evaluation and optimisation of + forecasts, predictions or models (2.5.0)". Zenodo. doi: 10.5281/zenodo.18638494 + """ + base, other = _sort_cubes_for_verification(cubes) + + # Copy the coordinates of the input cubes. + other_xr = xr.DataArray.from_iris(other) + base_xr = xr.DataArray.from_iris(base) + preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates) + + # Scores operates on xarray data arrays, so we transform the iris cube into an array, + # apply scores, and then transform it back. + scores_cube = xr.DataArray.to_iris( + scores.continuous.correlation.pearsonr( + other_xr, + base_xr, + preserve_dims=preserve_dims, + ) + ) + + # If time is aggregated out, attach a scalar time coordinate with bounds + # so plotting can display the aggregated period in the title. + try: + if not scores_cube.coords("time"): + base_time = base.coord("time") + time_vals = ( + base_time.bounds.flatten() + if base_time.has_bounds() + else base_time.points + ) + t_start = float(time_vals[0]) + t_end = float(time_vals[-1]) + t_mid = 0.5 * (t_start + t_end) + + scores_cube.add_aux_coord( + iris.coords.AuxCoord( + t_mid, + standard_name=base_time.standard_name, + long_name=base_time.long_name, + var_name=base_time.var_name, + units=base_time.units, + bounds=np.array([t_start, t_end]), + attributes=base_time.attributes.copy(), + ) + ) + except iris.exceptions.CoordinateNotFoundError: + pass + + scores_cube.rename(f"Pearson_Correlation_of_{base.name()}") + return scores_cube def scores_crps_for_ensemble( diff --git a/src/CSET/recipes/verification/surface_difference_scores_MAE.yaml b/src/CSET/recipes/verification/surface_difference_scores_MAE.yaml new file mode 100644 index 000000000..cc77db011 --- /dev/null +++ b/src/CSET/recipes/verification/surface_difference_scores_MAE.yaml @@ -0,0 +1,63 @@ +category: Scores +title: "$METHOD $VARNAME\nScores Mean Absolute Error between $OTHER_MODEL and $BASE_MODEL $SUBAREA_NAME" +description: | + Extracts and plots the Mean Absolute Error in $METHOD of $VARNAME. + 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 + RMSE 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: [$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 + + # The preserved coordinates specifies whether or not the + # MAE is to be computed over the entire case study. + - operator: scoreswrappers.scores_mae + preserved_coordinates: $PRESERVED_COORDS + + # This collapse is relevant only if the MAE has not been + # computed across the entire case study. + - 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/surface_difference_scores_RMSE.yaml b/src/CSET/recipes/verification/surface_difference_scores_RMSE.yaml new file mode 100644 index 000000000..9b4642cf4 --- /dev/null +++ b/src/CSET/recipes/verification/surface_difference_scores_RMSE.yaml @@ -0,0 +1,60 @@ +category: Scores +title: "$METHOD $VARNAME\nScores RMSE between $OTHER_MODEL and $BASE_MODEL $SUBAREA_NAME" +description: | + Extracts and plots the Root Mean Square Error in $METHOD of $VARNAME. + 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: [$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 + + # The preserved coordinates specifies whether or not the + # RMSE is to be computed over the entire case study. + - operator: scoreswrappers.scores_rmse + preserved_coordinates: $PRESERVED_COORDS + + # This collapse is relevant only if the RMSE has not been + # computed across the entire case study. + - 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/surface_difference_scores_additive_bias.yaml b/src/CSET/recipes/verification/surface_difference_scores_additive_bias.yaml new file mode 100644 index 000000000..9cc6b22ca --- /dev/null +++ b/src/CSET/recipes/verification/surface_difference_scores_additive_bias.yaml @@ -0,0 +1,58 @@ +category: Scores +title: "$METHOD $VARNAME\nScores Mean Error between $OTHER_MODEL and $BASE_MODEL $SUBAREA_NAME" +description: | + Extracts and plots the Mean Error (also known as the Additive Bias) in $METHOD of $VARNAME. + 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: [$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 + + # The preserved coordinates specifies whether or not the + # ME is to be computed over the entire case study. + - operator: scoreswrappers.scores_additive_bias + preserved_coordinates: $PRESERVED_COORDS + + # This collapse is relevant only if the ME has not been + # computed across the entire case study. + - 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_difference_scores_MAE.yaml b/src/CSET/recipes/verification/timeseries_surface_difference_scores_MAE.yaml new file mode 100644 index 000000000..d5973fbbb --- /dev/null +++ b/src/CSET/recipes/verification/timeseries_surface_difference_scores_MAE.yaml @@ -0,0 +1,44 @@ +category: Scores +title: "$VARNAME\nScores Mean Absolute Error timeseries between $OTHER_MODEL and $BASE_MODEL $SUBAREA_NAME" +description: | + Extracts and plots the Mean Absolute Error in $VARNAME + computed over the domain for each timestep. 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 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: [$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_mae + preserved_coordinates: ["time"] + + - operator: plot.plot_line_series + series_coordinate: time + + - operator: write.write_cube_to_nc + overwrite: True diff --git a/src/CSET/recipes/verification/surface_rmse_scores.yaml b/src/CSET/recipes/verification/timeseries_surface_difference_scores_RMSE.yaml similarity index 63% rename from src/CSET/recipes/verification/surface_rmse_scores.yaml rename to src/CSET/recipes/verification/timeseries_surface_difference_scores_RMSE.yaml index 5e834ea94..7915c332d 100644 --- a/src/CSET/recipes/verification/surface_rmse_scores.yaml +++ b/src/CSET/recipes/verification/timeseries_surface_difference_scores_RMSE.yaml @@ -1,13 +1,14 @@ category: Scores -title: "$METHOD $VARNAME\nRMSE between $OTHER_MODEL and $BASE_MODEL" +title: "$VARNAME\nScores RMSE timeseries between $OTHER_MODEL and $BASE_MODEL $SUBAREA_NAME" description: | - Extracts and plots the Root Mean Square Error in $METHOD of $VARNAME - for all times. 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 preserves the time, latitude and longitude coordinates. - This means a sequence of spatial plots of the RMSE can be produced. + Extracts and plots the Root Mean Square Error in $VARNAME + computed over the domain for each timestep. 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 allows the preservation of the time coordinate to produce a timeseries. Therefore, the RMSE is + collapsed over all other coordinates in the cube and calculated for every timestep. A larger RMSE implies a greater error than a smaller RMSE. An - RMSE of zero implies the two field match. The RMSE is calculated + 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). @@ -33,13 +34,10 @@ steps: subarea_extent: $SUBAREA_EXTENT - operator: scoreswrappers.scores_rmse - preserved_coordinates: ["time", "grid_latitude", "grid_longitude"] + preserved_coordinates: ["time"] - - operator: collapse.collapse - coordinate: ["time"] - method: $METHOD - - - operator: plot.spatial_pcolormesh_plot + - operator: plot.plot_line_series + series_coordinate: time - operator: write.write_cube_to_nc overwrite: True diff --git a/src/CSET/recipes/verification/timeseries_surface_difference_scores_additive_bias.yaml b/src/CSET/recipes/verification/timeseries_surface_difference_scores_additive_bias.yaml new file mode 100644 index 000000000..df9ee7c15 --- /dev/null +++ b/src/CSET/recipes/verification/timeseries_surface_difference_scores_additive_bias.yaml @@ -0,0 +1,42 @@ +category: Scores +title: "$VARNAME\nScores Mean Error timeseries between $OTHER_MODEL and $BASE_MODEL $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: [$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_additive_bias + preserved_coordinates: ["time"] + + - operator: plot.plot_line_series + series_coordinate: time + + - operator: write.write_cube_to_nc + overwrite: True diff --git a/src/CSET/recipes/verification/timeseries_surface_difference_scores_correlation_pearsonr.yaml b/src/CSET/recipes/verification/timeseries_surface_difference_scores_correlation_pearsonr.yaml new file mode 100644 index 000000000..179a735f9 --- /dev/null +++ b/src/CSET/recipes/verification/timeseries_surface_difference_scores_correlation_pearsonr.yaml @@ -0,0 +1,42 @@ +category: Scores +title: "$VARNAME\nScores Pearson's Correlation timeseries between $OTHER_MODEL and $BASE_MODEL $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: [$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_correlation_pearsonr + preserved_coordinates: ["time"] + + # Plot a timeseries of the Pearson Correlation + - operator: plot.plot_line_series + series_coordinate: time + + - operator: write.write_cube_to_nc + overwrite: True diff --git a/tests/operators/test_scoreswrappers.py b/tests/operators/test_scoreswrappers.py index 6f2f17799..1f035c3b1 100644 --- a/tests/operators/test_scoreswrappers.py +++ b/tests/operators/test_scoreswrappers.py @@ -36,6 +36,65 @@ ) +def test_scores_correlation_pearsonr(cube: iris.cube.Cube): + """Test taking the Pearson correlation between two cubes.""" + # Data preparation. + other_cube = cube.copy() + del other_cube.attributes["cset_comparison_base"] + cubes = iris.cube.CubeList([cube, other_cube]) + + # Take difference. + correlation_pearsonr_cube = scoreswrappers.scores_correlation_pearsonr(cubes) + + # As both cubes use the same data, check the Pearson correlation is one. + assert isinstance(correlation_pearsonr_cube, iris.cube.Cube) + assert np.allclose( + correlation_pearsonr_cube.data, + np.ones_like(correlation_pearsonr_cube.data), + atol=1e-9, + ) + assert correlation_pearsonr_cube.standard_name is None + assert ( + correlation_pearsonr_cube.long_name == "Pearson_Correlation_of_air_temperature" + ) + + +def test_scores_additive_bias(cube: iris.cube.Cube): + """Test taking the additive bias between two cubes.""" + # Data preparation. + other_cube = cube.copy() + del other_cube.attributes["cset_comparison_base"] + cubes = iris.cube.CubeList([cube, other_cube]) + + # Take difference. + additive_bias_cube = scoreswrappers.scores_additive_bias(cubes) + + # As both cubes use the same data, check the additive bias is zero. + assert isinstance(additive_bias_cube, iris.cube.Cube) + assert np.allclose( + additive_bias_cube.data, np.zeros_like(additive_bias_cube.data), atol=1e-9 + ) + assert additive_bias_cube.standard_name is None + assert additive_bias_cube.long_name == "Additive_Bias_of_air_temperature" + + +def test_scores_mae(cube: iris.cube.Cube): + """Test taking the mae between two cubes.""" + # Data preparation. + other_cube = cube.copy() + del other_cube.attributes["cset_comparison_base"] + cubes = iris.cube.CubeList([cube, other_cube]) + + # Take difference. + mae_cube = scoreswrappers.scores_mae(cubes) + + # As both cubes use the same data, check the mae is zero. + assert isinstance(mae_cube, iris.cube.Cube) + assert np.allclose(mae_cube.data, np.zeros_like(mae_cube.data), atol=1e-9) + assert mae_cube.standard_name is None + assert mae_cube.long_name == "MAE_of_air_temperature" + + def test_scores_rmse(cube: iris.cube.Cube): """Test taking the rmse between two cubes.""" # Data preparation. @@ -53,7 +112,7 @@ def test_scores_rmse(cube: iris.cube.Cube): assert rmse_cube.long_name == "RMSE_of_air_temperature" -def test_scores_rmse_nonzero(cube: iris.cube.Cube): +def test_scores_rmse_nonzero(): """Test taking the rmse between two different cubes.""" # Data preparation. cube = iris.cube.Cube(