Coverage for src/CSET/operators/scoreswrappers.py: 81%
106 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-16 13:33 +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 scores.probability
25import xarray as xr
26from iris.cube import Cube, CubeList
27from iris.util import reverse
29from CSET._common import is_increasing
30from CSET.operators._utils import fully_equalise_attributes, get_cube_yxcoordname
31from CSET.operators.constraints import (
32 generate_realization_constraint,
33 generate_remove_single_ensemble_member_constraint,
34)
35from CSET.operators.misc import _extract_common_time_points
36from CSET.operators.read import _realization_callback
37from CSET.operators.regrid import regrid_onto_cube
40def _sort_cubes_for_verification(cubes: CubeList):
41 """Prepare cubes ready for verification in scores.
43 Parameters
44 ----------
45 cubes: iris.cube.CubeList
46 A CubeList of exact 2 cubes, one from each model.
48 Returns
49 -------
50 base: iris.cube.Cube
51 The cube from the "analysis" in the same format as the other model.
52 other: iris.cube.Cube
53 The cube from the model in the same format as the base model.
55 Raises
56 ------
57 ValueError: "cubes should contain exactly 2 cubes."
58 If any other number of cubes are present.
60 Notes
61 -----
62 This operator is used for sorting the data into the correct format. It
63 is likely going to need to be refactored out of CSET and perhaps moved into
64 `CSET._utils` given common code between here and `misc.difference`.
65 """
66 # Set cubes into correct format using code from difference operator
67 if len(cubes) != 2:
68 raise ValueError("cubes should contain exactly 2 cubes.")
69 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
70 other: Cube = cubes.extract_cube(
71 iris.Constraint(
72 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
73 )
74 )
76 # If cubes contain a pressure coordinate, ensure it is increasing.
77 for cube in cubes:
78 try:
79 if len(cube.coord("pressure").points) > 2: 79 ↛ 77line 79 didn't jump to line 77 because the condition on line 79 was always true
80 if not is_increasing(cube.coord("pressure").points): 80 ↛ 81line 80 didn't jump to line 81 because the condition on line 80 was never true
81 reverse(cube, "pressure")
83 except iris.exceptions.CoordinateNotFoundError:
84 pass
86 # Extract just common time points.
87 base, other = _extract_common_time_points(base, other)
89 # Get spatial coord names.
90 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
91 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
93 # Ensure cubes to compare are on common differencing grid.
94 # This is triggered if either
95 # i) latitude and longitude shapes are not the same. Note grid points
96 # are not compared directly as these can differ through rounding
97 # errors.
98 # ii) or variables are known to often sit on different grid staggering
99 # in different models (e.g. cell center vs cell edge), as is the case
100 # for UM and LFRic comparisons.
101 # In future greater choice of regridding method might be applied depending
102 # on variable type. Linear regridding can in general be appropriate for smooth
103 # variables. Care should be taken with interpretation of differences
104 # given this dependency on regridding.
105 if (
106 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
107 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
108 ) or (
109 base.long_name
110 in [
111 "eastward_wind_at_10m",
112 "northward_wind_at_10m",
113 "northward_wind_at_cell_centres",
114 "eastward_wind_at_cell_centres",
115 "zonal_wind_at_pressure_levels",
116 "meridional_wind_at_pressure_levels",
117 "potential_vorticity_at_pressure_levels",
118 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
119 ]
120 ):
121 logging.debug(
122 "Linear regridding base cube to other grid to compute differences"
123 )
124 base = regrid_onto_cube(base, other, method="Linear")
126 # Figure out if we are comparing between UM and LFRic; flip array if so.
127 base_lat_direction = is_increasing(base.coord(base_lat_name).points)
128 other_lat_direction = is_increasing(other.coord(other_lat_name).points)
129 if base_lat_direction != other_lat_direction: 129 ↛ 131line 129 didn't jump to line 131 because the condition on line 129 was never true
130 # Copy base cube for correct coordinate information.
131 other_tmp = base.copy()
132 # Flip the data and place in the copied cube.
133 other_tmp.data = np.flip(
134 other.data, other.coord(other_lat_name).cube_dims(other)
135 )
136 # Use original name and units from the other cube.
137 other_tmp.rename(other.name())
138 other_tmp.units = other.units
139 # Replace the cube.
140 other = other_tmp
142 # Equalise attributes so we can merge.
143 fully_equalise_attributes(CubeList([base, other]))
144 logging.debug("Base: %s\nOther: %s", base, other)
146 return base, other
149def _resolve_preserve_dims(
150 cube: Cube,
151 data_array: xr.DataArray,
152 preserved_coordinates: list[str] | str | None,
153) -> list[str] | None:
154 """Resolve preserve coordinates to xarray dimension names.
156 The ``scores`` package expects preserve dimensions to match xarray
157 dimension names. In Iris data, commonly used coordinates such as ``time``
158 may be auxiliary coordinates attached to a differently named dimension
159 (e.g. ``dim0``). This helper maps coordinate names to their underlying
160 dimension names.
161 """
162 if preserved_coordinates is None:
163 return None
165 coord_names = (
166 [preserved_coordinates]
167 if isinstance(preserved_coordinates, str)
168 else preserved_coordinates
169 )
170 preserve_dims: list[str] = []
172 for coord_name in coord_names:
173 # Already an xarray dimension name.
174 if coord_name in data_array.dims: 174 ↛ 180line 174 didn't jump to line 180 because the condition on line 174 was always true
175 if coord_name not in preserve_dims: 175 ↛ 177line 175 didn't jump to line 177 because the condition on line 175 was always true
176 preserve_dims.append(coord_name)
177 continue
179 # Otherwise, map coordinate name to dimension index/indices.
180 try:
181 dim_indices = cube.coord_dims(coord_name)
182 except iris.exceptions.CoordinateNotFoundError:
183 # Keep original name so scores raises a clear error for unknown keys.
184 if coord_name not in preserve_dims:
185 preserve_dims.append(coord_name)
186 continue
188 for dim_index in dim_indices:
189 dim_name = data_array.dims[dim_index]
190 if dim_name not in preserve_dims:
191 preserve_dims.append(dim_name)
193 return preserve_dims
196def scores_rmse(cubes: CubeList, preserved_coordinates: list[str] | str | None = None):
197 r"""Calculate the Root Mean Square Error (RMSE) using scores.
199 Acts as a wrapper around the RMSE calculation from ``scores`` ([scoresa]_, [scoresb]_).
200 It is calculated as
202 .. math:: RMSE = \sqrt{\frac{1}{N} \Sigma(forecast - observations)^2}
204 Parameters
205 ----------
206 cubes: iris.cube.CubeList
207 A CubeList containing exactly two cubes: a base and an "other" model,
208 this can be an analysis and the model.
209 preserved_coordinates: list[str] | str | None, default is None.
210 The coordinates (or xarray dimension names) that you wish to preserve in the calculaiton of the
211 RMSE. For example if you want a map of each time you can preserve
212 ["time","grid_latitude", "grid_longitude"] or if you want a time series
213 you can preserve ["time"], if you want to collapse to a single value
214 use `None`. The default is `None`.
216 Returns
217 -------
218 RMSE: iris.cube.Cube
219 A cube containing the RMSE between the base and other cube.
221 References
222 ----------
223 .. [scoresa] Leeuwenburg, T., Loveday, N., Ebert, E. E., Cook, H.,
224 Khanarmuei, M., Taggart, R. J., Ramanathan, N., Carroll, M., Chong, S.,
225 Griffiths, A., & Sharples, J. (2024) "scores: A Python package for
226 verifying and evaluating models and predictions with xarray". Journal
227 of Open Source Software, vol. 9, 6889. doi: 10.21105/joss.06889
229 .. [scoresb] Leeuwenburg, T., Loveday, N., Ramanathan, N., Chong, S.,
230 Taggart, R. J., Shrestha, D., Khanarmuei, M., Cook, H., Bluett, L., Ebert,
231 E. E., Carroll, M., Trotta, B., Bishop, S., Squire, D. T., Griffiths, A.,
232 Pagano, T. C., Fisher, A. J., Mandelbaum, T., Jinghan, F., … Smallwood, J.
233 (2026) "scores: Metrics for the verification, evaluation and optimisation of
234 forecasts, predictions or models (2.5.0)". Zenodo. doi: 10.5281/zenodo.18638494
235 """
236 base, other = _sort_cubes_for_verification(cubes)
238 other_xr = xr.DataArray.from_iris(other)
239 base_xr = xr.DataArray.from_iris(base)
240 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
242 # Scores operators on xarray data arrays, so we transform the iris cube into an array,
243 # apply scores, and then transform it back.
244 RMSE = xr.DataArray.to_iris(
245 scores.continuous.rmse(
246 other_xr,
247 base_xr,
248 preserve_dims=preserve_dims,
249 )
250 )
252 # If time is aggregated out, attach a scalar time coordinate with bounds
253 # so plotting can display the aggregated period in the title.
254 try:
255 if not RMSE.coords("time"):
256 base_time = base.coord("time")
257 time_vals = (
258 base_time.bounds.flatten()
259 if base_time.has_bounds()
260 else base_time.points
261 )
262 t_start = float(time_vals[0])
263 t_end = float(time_vals[-1])
264 t_mid = 0.5 * (t_start + t_end)
266 RMSE.add_aux_coord(
267 iris.coords.AuxCoord(
268 t_mid,
269 standard_name=base_time.standard_name,
270 long_name=base_time.long_name,
271 var_name=base_time.var_name,
272 units=base_time.units,
273 bounds=np.array([t_start, t_end]),
274 attributes=base_time.attributes.copy(),
275 )
276 )
277 except iris.exceptions.CoordinateNotFoundError:
278 pass
280 RMSE.rename(f"RMSE_of_{base.name()}")
281 return RMSE
284def scores_crps_for_ensemble(
285 cubes: Cube | CubeList, method: str = "ecdf", control_member: int = 0
286) -> iris.Constraint:
287 r"""Calculate the CRPS for an ensemble.
289 Acts as a wrapper around the crps_for_ensemble from ``scores`` ([scores_a]_, [scores_b]_).
291 Lower CRPS values are better (implies experiment distribution is closer to control distribution/observations),
292 larger values are worse (implies distributions are dissimilar).
293 It is applicable across time and spatial scales as the focus is on the distribution of the values.
294 Default method is ecdf. ecdf is exact value from the empirical distributions,
295 whereas fair produces an approximated value based on a random sample of the underlying distribution.
297 See [CRPS] for further information.
299 Parameters
300 ----------
301 cubes: iris.cube.Cube
302 A Cube containing ensembles data
304 Returns
305 -------
306 crps: iris.cube.Cube
307 A cube containing the crps between the ensemble members and the control
309 References
310 ----------
311 .. [scores_a] Leeuwenburg, T., Loveday, N., Ebert, E. E., Cook, H.,
312 Khanarmuei, M., Taggart, R. J., Ramanathan, N., Carroll, M., Chong, S.,
313 Griffiths, A., & Sharples, J. (2024) "scores: A Python package for
314 verifying and evaluating models and predictions with xarray". Journal
315 of Open Source Software, vol. 9, 6889. doi: 10.21105/joss.06889
317 .. [scores_b] Leeuwenburg, T., Loveday, N., Ramanathan, N., Chong, S.,
318 Taggart, R. J., Shrestha, D., Khanarmuei, M., Cook, H., Bluett, L., Ebert,
319 E. E., Carroll, M., Trotta, B., Bishop, S., Squire, D. T., Griffiths, A.,
320 Pagano, T. C., Fisher, A. J., Mandelbaum, T., Jinghan, F., … Smallwood, J.
321 (2026) "scores: Metrics for the verification, evaluation and optimisation of
322 forecasts, predictions or models (2.5.0)". Zenodo. doi: 10.5281/zenodo.18638494
324 .. [CRPS]
325 Hersbach, H., 2000: Decomposition of the Continuous Ranked
326 Probability Score for Ensemble Prediction Systems. Wea.
327 Forecasting, 15, 559–570, https://doi.org/10.1175/1520-0434(2000)015<0559:DOTCRP>2.0.CO;2.
328 """
329 if control_member != 0:
330 logging.warning("control member is usual 0")
332 if control_member not in cubes.coords("realization")[0].points:
333 new_control_member = cubes.coords("realization")[0].points[0]
334 logging.warning(
335 f"control member value {control_member} out of bounds, defaulting to control member={new_control_member}"
336 )
337 control_member = new_control_member
339 if cubes.coord("time").shape[0] == 1:
340 raise ValueError("Cube has only one time point.")
342 if cubes.coord("realization").shape[0] < 3:
343 raise ValueError("Cube should have one control member and at least two members")
345 ctrl = cubes.extract(generate_realization_constraint([control_member]))
346 ens_mem = cubes.extract(
347 generate_remove_single_ensemble_member_constraint(control_member)
348 )
350 # Realising the data in advance provides a large speedup
351 _ = ctrl.data
352 _ = ens_mem.data
353 del _
355 ctrl = xr.DataArray.from_iris(ctrl)
356 ens_mem = xr.DataArray.from_iris(ens_mem)
358 crps = xr.DataArray.to_iris(
359 scores.probability.crps_for_ensemble(
360 ens_mem,
361 ctrl,
362 ensemble_member_dim="realization",
363 method=method,
364 preserve_dims="time",
365 )
366 )
368 crps.rename(f"CRPS_of_{cubes[0].name()}")
369 _realization_callback(crps)
370 return crps