Coverage for src/CSET/operators/scoreswrappers.py: 78%
96 statements
« prev ^ index » next coverage.py v7.15.0, created at 2026-07-09 15:26 +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 numpy as np
22import scores
23import scores.continuous
24import xarray as xr
25from iris.cube import Cube, CubeList
26from iris.util import reverse
28from CSET._common import is_increasing
29from CSET.operators._utils import fully_equalise_attributes, get_cube_yxcoordname
30from CSET.operators.misc import _extract_common_time_points
31from CSET.operators.regrid import regrid_onto_cube
34def _sort_cubes_for_verification(cubes: CubeList):
35 """Prepare cubes ready for verification in scores.
37 Parameters
38 ----------
39 cubes: iris.cube.CubeList
40 A CubeList of exact 2 cubes, one from each model.
42 Returns
43 -------
44 base: iris.cube.Cube
45 The cube from the "analysis" in the same format as the other model.
46 other: iris.cube.Cube
47 The cube from the model in the same format as the base model.
49 Raises
50 ------
51 ValueError: "cubes should contain exactly 2 cubes."
52 If any other number of cubes are present.
54 Notes
55 -----
56 This operator is used for sorting the data into the correct format. It
57 is likely going to need to be refactored out of CSET and perhaps moved into
58 `CSET._utils` given common code between here and `misc.difference`.
59 """
60 # Set cubes into correct format using code from difference operator
61 if len(cubes) != 2:
62 raise ValueError("cubes should contain exactly 2 cubes.")
63 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
64 other: Cube = cubes.extract_cube(
65 iris.Constraint(
66 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
67 )
68 )
70 # If cubes contain a pressure coordinate, ensure it is increasing.
71 for cube in cubes:
72 try:
73 if len(cube.coord("pressure").points) > 2: 73 ↛ 71line 73 didn't jump to line 71 because the condition on line 73 was always true
74 if not is_increasing(cube.coord("pressure").points): 74 ↛ 75line 74 didn't jump to line 75 because the condition on line 74 was never true
75 reverse(cube, "pressure")
77 except iris.exceptions.CoordinateNotFoundError:
78 pass
80 # Extract just common time points.
81 base, other = _extract_common_time_points(base, other)
83 # Get spatial coord names.
84 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
85 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
87 # Ensure cubes to compare are on common differencing grid.
88 # This is triggered if either
89 # i) latitude and longitude shapes are not the same. Note grid points
90 # are not compared directly as these can differ through rounding
91 # errors.
92 # ii) or variables are known to often sit on different grid staggering
93 # in different models (e.g. cell center vs cell edge), as is the case
94 # for UM and LFRic comparisons.
95 # In future greater choice of regridding method might be applied depending
96 # on variable type. Linear regridding can in general be appropriate for smooth
97 # variables. Care should be taken with interpretation of differences
98 # given this dependency on regridding.
99 if (
100 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
101 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
102 ) or (
103 base.long_name
104 in [
105 "eastward_wind_at_10m",
106 "northward_wind_at_10m",
107 "northward_wind_at_cell_centres",
108 "eastward_wind_at_cell_centres",
109 "zonal_wind_at_pressure_levels",
110 "meridional_wind_at_pressure_levels",
111 "potential_vorticity_at_pressure_levels",
112 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
113 ]
114 ):
115 logging.debug(
116 "Linear regridding base cube to other grid to compute differences"
117 )
118 base = regrid_onto_cube(base, other, method="Linear")
120 # Figure out if we are comparing between UM and LFRic; flip array if so.
121 base_lat_direction = is_increasing(base.coord(base_lat_name).points)
122 other_lat_direction = is_increasing(other.coord(other_lat_name).points)
123 if base_lat_direction != other_lat_direction: 123 ↛ 125line 123 didn't jump to line 125 because the condition on line 123 was never true
124 # Copy base cube for correct coordinate information.
125 other_tmp = base.copy()
126 # Flip the data and place in the copied cube.
127 other_tmp.data = np.flip(
128 other.data, other.coord(other_lat_name).cube_dims(other)
129 )
130 # Use original name and units from the other cube.
131 other_tmp.rename(other.name())
132 other_tmp.units = other.units
133 # Replace the cube.
134 other = other_tmp
136 # Equalise attributes so we can merge.
137 fully_equalise_attributes(CubeList([base, other]))
138 logging.debug("Base: %s\nOther: %s", base, other)
140 return base, other
143def _resolve_preserve_dims(
144 cube: Cube,
145 data_array: xr.DataArray,
146 preserved_coordinates: list[str] | str | None,
147) -> list[str] | None:
148 """Resolve preserve coordinates to xarray dimension names.
150 The ``scores`` package expects preserve dimensions to match xarray
151 dimension names. In Iris data, commonly used coordinates such as ``time``
152 may be auxiliary coordinates attached to a differently named dimension
153 (e.g. ``dim0``). This helper maps coordinate names to their underlying
154 dimension names.
155 """
156 if preserved_coordinates is None:
157 return None
159 coord_names = (
160 [preserved_coordinates]
161 if isinstance(preserved_coordinates, str)
162 else preserved_coordinates
163 )
164 preserve_dims: list[str] = []
166 for coord_name in coord_names:
167 # Already an xarray dimension name.
168 if coord_name in data_array.dims: 168 ↛ 174line 168 didn't jump to line 174 because the condition on line 168 was always true
169 if coord_name not in preserve_dims: 169 ↛ 171line 169 didn't jump to line 171 because the condition on line 169 was always true
170 preserve_dims.append(coord_name)
171 continue
173 # Otherwise, map coordinate name to dimension index/indices.
174 try:
175 dim_indices = cube.coord_dims(coord_name)
176 except iris.exceptions.CoordinateNotFoundError:
177 # Keep original name so scores raises a clear error for unknown keys.
178 if coord_name not in preserve_dims:
179 preserve_dims.append(coord_name)
180 continue
182 for dim_index in dim_indices:
183 dim_name = data_array.dims[dim_index]
184 if dim_name not in preserve_dims:
185 preserve_dims.append(dim_name)
187 return preserve_dims
190def scores_rmse(cubes: CubeList, preserved_coordinates: list[str] | str | None = None):
191 r"""Calculate the Root Mean Square Error (RMSE) using scores.
193 Acts as a wrapper around the RMSE calculation from ``scores`` ([scoresa]_, [scoresb]_).
194 It is calculated as
196 .. math:: RMSE = \sqrt{\frac{1}{N} \Sigma(forecast - observations)^2}
198 Parameters
199 ----------
200 cubes: iris.cube.CubeList
201 A CubeList containing exactly two cubes: a base and an "other" model,
202 this can be an analysis and the model.
203 preserved_coordinates: list[str] | str | None, default is None.
204 The coordinates that you wish to preserve in the calculaiton of the
205 RMSE. For example if you want a map of each time you can preserve
206 ["time","grid_latitude", "grid_longitude"] or if you want a time series
207 you can preserve ["time"], if you want to collapse to a single value
208 use `None`. The default is `None`.
210 Returns
211 -------
212 scores_cube: iris.cube.Cube
213 A cube containing the RMSE between the base and other cube.
215 References
216 ----------
217 .. [scoresa] Leeuwenburg, T., Loveday, N., Ebert, E. E., Cook, H.,
218 Khanarmuei, M., Taggart, R. J., Ramanathan, N., Carroll, M., Chong, S.,
219 Griffiths, A., & Sharples, J. (2024) "scores: A Python package for
220 verifying and evaluating models and predictions with xarray". Journal
221 of Open Source Software, vol. 9, 6889. doi: 10.21105/joss.06889
223 .. [scoresb] Leeuwenburg, T., Loveday, N., Ramanathan, N., Chong, S.,
224 Taggart, R. J., Shrestha, D., Khanarmuei, M., Cook, H., Bluett, L., Ebert,
225 E. E., Carroll, M., Trotta, B., Bishop, S., Squire, D. T., Griffiths, A.,
226 Pagano, T. C., Fisher, A. J., Mandelbaum, T., Jinghan, F., … Smallwood, J.
227 (2026) "scores: Metrics for the verification, evaluation and optimisation of
228 forecasts, predictions or models (2.5.0)". Zenodo. doi: 10.5281/zenodo.18638494
229 """
230 base, other = _sort_cubes_for_verification(cubes)
232 # SB
233 other_xr = xr.DataArray.from_iris(other)
234 base_xr = xr.DataArray.from_iris(base)
235 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
237 # if preserved_coordinates == ["grid_latitude", "grid_longitude"]:
238 # time_coord = base.coord[0]("time")
239 # Scores operators on xarray data arrays, so we transform the iris cube into an array,
240 # apply scores, and then transform it back.
241 scores_cube = xr.DataArray.to_iris(
242 scores.continuous.rmse(
243 other_xr,
244 base_xr,
245 preserve_dims=preserve_dims,
246 )
247 )
249 # If time is aggregated out, attach a scalar time coordinate with bounds
250 # so plotting can display the aggregated period in the title.
251 try:
252 if not scores_cube.coords("time"):
253 base_time = base.coord("time")
254 time_vals = (
255 base_time.bounds.flatten()
256 if base_time.has_bounds()
257 else base_time.points
258 )
259 t_start = float(time_vals[0])
260 t_end = float(time_vals[-1])
261 t_mid = 0.5 * (t_start + t_end)
263 scores_cube.add_aux_coord(
264 iris.coords.AuxCoord(
265 t_mid,
266 standard_name=base_time.standard_name,
267 long_name=base_time.long_name,
268 var_name=base_time.var_name,
269 units=base_time.units,
270 bounds=np.array([t_start, t_end]),
271 attributes=base_time.attributes.copy(),
272 )
273 )
274 except iris.exceptions.CoordinateNotFoundError:
275 pass
277 scores_cube.rename(f"RMSE_of_{base.name()}")
278 # if preserved_coordinates == ["grid_latitude", "grid_longitude"]:
279 # scores_cube.add_aux_coord(time_coord)
280 return scores_cube
283def scores_mae(cubes: CubeList, preserved_coordinates: list[str] | str | None = None):
284 """Calculate the Mean Absolute Error (MAE) using scores."""
285 base, other = _sort_cubes_for_verification(cubes)
286 # Scores operators on xarray data arrays, so we transform the iris cube into an array,
287 # apply scores, and then transform it back.
288 scores_cube = xr.DataArray.to_iris(
289 scores.continuous.mae(
290 xr.DataArray.from_iris(other),
291 xr.DataArray.from_iris(base),
292 preserve_dims=preserved_coordinates,
293 )
294 )
295 scores_cube.rename(f"MAE_of_{base.name()}")
296 return scores_cube
299def scores_additive_bias(
300 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
301):
302 """Calculate the Additive Bias (Mean Error) using scores."""
303 base, other = _sort_cubes_for_verification(cubes)
304 # Scores operators on xarray data arrays, so we transform the iris cube into an array,
305 # apply scores, and then transform it back.
306 scores_cube = xr.DataArray.to_iris(
307 scores.continuous.additive_bias(
308 xr.DataArray.from_iris(other),
309 xr.DataArray.from_iris(base),
310 preserve_dims=preserved_coordinates,
311 )
312 )
313 scores_cube.rename(f"Additive_Bias_of_{base.name()}")
314 return scores_cube
317def scores_correlation_pearsonr(
318 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
319):
320 """Calculate the Pearson's Correlation using scores."""
321 base, other = _sort_cubes_for_verification(cubes)
322 # Scores operators on xarray data arrays, so we transform the iris cube into an array,
323 # apply scores, and then transform it back.
324 scores_cube = xr.DataArray.to_iris(
325 scores.continuous.correlation.pearsonr(
326 xr.DataArray.from_iris(other),
327 xr.DataArray.from_iris(base),
328 preserve_dims=preserved_coordinates,
329 )
330 )
331 scores_cube.rename(f"Pearson_Correlation_of_{base.name()}")
332 return scores_cube