Coverage for src/CSET/operators/scoreswrappers.py: 91%
42 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-18 14:09 +0000
1# © Crown copyright, Met Office (2022-2026) and CSET contributors.
2#
3# Licensed under the Apache License, Version 2.0 (the "License");
4# you may not use this file except in compliance with the License.
5# You may obtain a copy of the License at
6#
7# http://www.apache.org/licenses/LICENSE-2.0
8#
9# Unless required by applicable law or agreed to in writing, software
10# distributed under the License is distributed on an "AS IS" BASIS,
11# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12# See the License for the specific language governing permissions and
13# limitations under the License.
15"""A module containing wrappers for the scores module."""
17import logging
19import iris
20import iris.exceptions
21import scores
22import scores.continuous
23import xarray as xr
24from iris.cube import Cube, CubeList
25from iris.util import reverse
27from CSET._common import is_increasing
28from CSET.operators._utils import fully_equalise_attributes, get_cube_yxcoordname
29from CSET.operators.misc import _extract_common_time_points
30from CSET.operators.regrid import regrid_onto_cube
33def _sort_cubes_for_verification(cubes: CubeList):
34 """Prepare cubes ready for verification in scores.
36 Parameters
37 ----------
38 cubes: iris.cube.CubeList
39 A CubeList of exact 2 cubes, one from each model.
41 Returns
42 -------
43 base: iris.cube.Cube
44 The cube from the "analysis" in the same format as the other model.
45 other: iris.cube.Cube
46 The cube from the model in the same format as the base model.
48 Raises
49 ------
50 ValueError: "cubes should contain exactly 2 cubes."
51 If any other number of cubes are present.
53 Notes
54 -----
55 This operator is used for sorting the data into the correct format. It
56 is likely going to need to be refactored out of CSET and perhaps moved into
57 `CSET._utils` given common code between here and `misc.difference`.
58 """
59 # Set cubes into correct format using code from difference operator
60 if len(cubes) != 2:
61 raise ValueError("cubes should contain exactly 2 cubes.")
62 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
63 other: Cube = cubes.extract_cube(
64 iris.Constraint(
65 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
66 )
67 )
69 # If cubes contain a pressure coordinate, ensure it is increasing.
70 for cube in cubes:
71 try:
72 if len(cube.coord("pressure").points) > 2: 72 ↛ 70line 72 didn't jump to line 70 because the condition on line 72 was always true
73 if not is_increasing(cube.coord("pressure").points): 73 ↛ 74line 73 didn't jump to line 74 because the condition on line 73 was never true
74 reverse(cube, "pressure")
76 except iris.exceptions.CoordinateNotFoundError:
77 pass
79 # Get spatial coord names.
80 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
81 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
83 # Ensure cubes to compare are on common differencing grid.
84 # This is triggered if either
85 # i) latitude and longitude shapes are not the same. Note grid points
86 # are not compared directly as these can differ through rounding
87 # errors.
88 # ii) or variables are known to often sit on different grid staggering
89 # in different models (e.g. cell center vs cell edge), as is the case
90 # for UM and LFRic comparisons.
91 # In future greater choice of regridding method might be applied depending
92 # on variable type. Linear regridding can in general be appropriate for smooth
93 # variables. Care should be taken with interpretation of differences
94 # given this dependency on regridding.
95 if (
96 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
97 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
98 ) or (
99 base.long_name
100 in [
101 "eastward_wind_at_10m",
102 "northward_wind_at_10m",
103 "northward_wind_at_cell_centres",
104 "eastward_wind_at_cell_centres",
105 "zonal_wind_at_pressure_levels",
106 "meridional_wind_at_pressure_levels",
107 "potential_vorticity_at_pressure_levels",
108 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
109 ]
110 ):
111 logging.debug(
112 "Linear regridding base cube to other grid to compute differences"
113 )
114 base = regrid_onto_cube(base, other, method="Linear")
116 # Figure out if we are comparing between UM and LFRic; flip array if so.
117 base_lat_direction = is_increasing(base.coord(base_lat_name).points)
118 other_lat_direction = is_increasing(other.coord(other_lat_name).points)
119 if base_lat_direction != other_lat_direction: 119 ↛ 120line 119 didn't jump to line 120 because the condition on line 119 was never true
120 reverse(other, other_lat_name)
122 # Extract just common time points.
123 base, other = _extract_common_time_points(base, other)
125 # Equalise attributes so we can merge.
126 fully_equalise_attributes(CubeList([base, other]))
127 logging.debug("Base: %s\nOther: %s", base, other)
129 return base, other
132def scores_rmse(cubes: CubeList, preserved_coordinates: list[str] | str | None = None):
133 r"""Calculate the Root Mean Square Error (RMSE) using scores.
135 Acts as a wrapper around the RMSE calculation from ``scores`` ([scoresa]_, [scoresb]_).
136 It is calculated as
138 .. math:: RMSE = \sqrt{\frac{1}{N} \Sigma(forecast - observations)^2}
140 Parameters
141 ----------
142 cubes: iris.cube.CubeList
143 A CubeList containing exactly two cubes: a base and an "other" model,
144 this can be an analysis and the model.
145 preserved_coordinates: list[str] | str | None, default is None.
146 The coordinates that you wish to preserve in the calculaiton of the
147 RMSE. For example if you want a map of each time you can preserve
148 ["time","grid_latitude", "grid_longitude"] or if you want a time series
149 you can preserve ["time"], if you want to collapse to a single value
150 use `None`. The default is `None`.
152 Returns
153 -------
154 RMSE: iris.cube.Cube
155 A cube containing the RMSE between the base and other cube.
157 References
158 ----------
159 .. [scoresa] Leeuwenburg, T., Loveday, N., Ebert, E. E., Cook, H.,
160 Khanarmuei, M., Taggart, R. J., Ramanathan, N., Carroll, M., Chong, S.,
161 Griffiths, A., & Sharples, J. (2024) "scores: A Python package for
162 verifying and evaluating models and predictions with xarray". Journal
163 of Open Source Software, vol. 9, 6889. doi: 10.21105/joss.06889
165 .. [scoresb] Leeuwenburg, T., Loveday, N., Ramanathan, N., Chong, S.,
166 Taggart, R. J., Shrestha, D., Khanarmuei, M., Cook, H., Bluett, L., Ebert,
167 E. E., Carroll, M., Trotta, B., Bishop, S., Squire, D. T., Griffiths, A.,
168 Pagano, T. C., Fisher, A. J., Mandelbaum, T., Jinghan, F., … Smallwood, J.
169 (2026) "scores: Metrics for the verification, evaluation and optimisation of
170 forecasts, predictions or models (2.5.0)". Zenodo. doi: 10.5281/zenodo.18638494
171 """
172 base, other = _sort_cubes_for_verification(cubes)
173 # Scores operators on xarray data arrays, so we transform the iris cube into an array,
174 # apply scores, and then transform it back.
175 RMSE = xr.DataArray.to_iris(
176 scores.continuous.rmse(
177 xr.DataArray.from_iris(other),
178 xr.DataArray.from_iris(base),
179 preserve_dims=preserved_coordinates,
180 )
181 )
182 RMSE.rename(f"RMSE_of_{base.name()}")
183 return RMSE