Coverage for src/CSET/operators/scoreswrappers.py: 82%
160 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-07 15:12 +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
39logger = logging.getLogger(__name__)
42def _sort_cubes_for_verification(cubes: CubeList):
43 """Prepare cubes ready for verification in scores.
45 Parameters
46 ----------
47 cubes: iris.cube.CubeList
48 A CubeList of exact 2 cubes, one from each model.
50 Returns
51 -------
52 base: iris.cube.Cube
53 The cube from the "analysis" in the same format as the other model.
54 other: iris.cube.Cube
55 The cube from the model in the same format as the base model.
57 Raises
58 ------
59 ValueError: "cubes should contain exactly 2 cubes."
60 If any other number of cubes are present.
62 Notes
63 -----
64 This operator is used for sorting the data into the correct format. It
65 is likely going to need to be refactored out of CSET and perhaps moved into
66 `CSET._utils` given common code between here and `misc.difference`.
67 """
68 # Set cubes into correct format using code from difference operator
69 if len(cubes) != 2:
70 raise ValueError("cubes should contain exactly 2 cubes.")
71 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
72 other: Cube = cubes.extract_cube(
73 iris.Constraint(
74 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
75 )
76 )
78 # If cubes contain a pressure coordinate, ensure it is increasing.
79 for cube in cubes:
80 try:
81 if len(cube.coord("pressure").points) > 2 and not is_increasing( 81 ↛ 84line 81 didn't jump to line 84 because the condition on line 81 was never true
82 cube.coord("pressure").points
83 ):
84 reverse(cube, "pressure")
86 except iris.exceptions.CoordinateNotFoundError:
87 pass
89 # Extract just common time points.
90 base, other = _extract_common_time_points(base, other)
92 # Get spatial coord names.
93 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
94 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
96 # Ensure cubes to compare are on common differencing grid.
97 # This is triggered if either
98 # i) latitude and longitude shapes are not the same. Note grid points
99 # are not compared directly as these can differ through rounding
100 # errors.
101 # ii) or variables are known to often sit on different grid staggering
102 # in different models (e.g. cell center vs cell edge), as is the case
103 # for UM and LFRic comparisons.
104 # In future greater choice of regridding method might be applied depending
105 # on variable type. Linear regridding can in general be appropriate for smooth
106 # variables. Care should be taken with interpretation of differences
107 # given this dependency on regridding.
108 if (
109 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
110 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
111 ) or (
112 base.long_name
113 in [
114 "eastward_wind_at_10m",
115 "northward_wind_at_10m",
116 "northward_wind_at_cell_centres",
117 "eastward_wind_at_cell_centres",
118 "zonal_wind_at_pressure_levels",
119 "meridional_wind_at_pressure_levels",
120 "potential_vorticity_at_pressure_levels",
121 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
122 ]
123 ):
124 logger.debug("Linear regridding base cube to other grid to compute differences")
125 base = regrid_onto_cube(base, other, method="Linear")
127 # Figure out if we are comparing between UM and LFRic; flip array if so.
128 base_lat_direction = is_increasing(base.coord(base_lat_name).points)
129 other_lat_direction = is_increasing(other.coord(other_lat_name).points)
130 if base_lat_direction != other_lat_direction: 130 ↛ 132line 130 didn't jump to line 132 because the condition on line 130 was never true
131 # Copy base cube for correct coordinate information.
132 other_tmp = base.copy()
133 # Flip the data and place in the copied cube.
134 other_tmp.data = np.flip(
135 other.data, other.coord(other_lat_name).cube_dims(other)
136 )
137 # Use original name and units from the other cube.
138 other_tmp.rename(other.name())
139 other_tmp.units = other.units
140 # Replace the cube.
141 other = other_tmp
143 # Equalise attributes so we can merge.
144 fully_equalise_attributes(CubeList([base, other]))
145 logger.debug("Base: %s\nOther: %s", base, other)
147 return base, other
150def _resolve_preserve_dims(
151 cube: Cube,
152 data_array: xr.DataArray,
153 preserved_coordinates: list[str] | str | None,
154) -> list[str] | None:
155 """Resolve preserve coordinates to xarray dimension names.
157 The ``scores`` package expects preserve dimensions to match xarray
158 dimension names. In Iris data, commonly used coordinates such as ``time``
159 may be auxiliary coordinates attached to a differently named dimension
160 (e.g. ``dim0``). This helper maps coordinate names to their underlying
161 dimension names and helps to convert from iris to xarray coordinate dimension names.
162 """
163 if preserved_coordinates is None:
164 return None
166 coord_names = (
167 [preserved_coordinates]
168 if isinstance(preserved_coordinates, str)
169 else preserved_coordinates
170 )
171 preserve_dims: list[str] = []
173 for coord_name in coord_names:
174 # Already an xarray dimension name.
175 if coord_name in data_array.dims: 175 ↛ 181line 175 didn't jump to line 181 because the condition on line 175 was always true
176 if coord_name not in preserve_dims: 176 ↛ 178line 176 didn't jump to line 178 because the condition on line 176 was always true
177 preserve_dims.append(coord_name)
178 continue
180 # Otherwise, map coordinate name to dimension index/indices.
181 try:
182 dim_indices = cube.coord_dims(coord_name)
183 except iris.exceptions.CoordinateNotFoundError:
184 # Keep original name so scores raises a clear error for unknown keys.
185 if coord_name not in preserve_dims:
186 preserve_dims.append(coord_name)
187 continue
189 for dim_index in dim_indices:
190 dim_name = data_array.dims[dim_index]
191 if dim_name not in preserve_dims:
192 preserve_dims.append(dim_name)
194 return preserve_dims
197def scores_rmse(cubes: CubeList, preserved_coordinates: list[str] | str | None = None):
198 r"""Calculate the Root Mean Square Error (RMSE) using scores.
200 Acts as a wrapper around the RMSE calculation from ``scores`` ([scoresa]_, [scoresb]_).
201 It is calculated as
203 .. math:: RMSE = \sqrt{\frac{1}{N} \Sigma(forecast - observations)^2}
205 Parameters
206 ----------
207 cubes: iris.cube.CubeList
208 A CubeList containing exactly two cubes: a base and an "other" model,
209 this can be an analysis and the model.
210 preserved_coordinates: list[str] | str | None, default is None.
211 The coordinates (or xarray dimension names) that you wish to preserve in the calculaiton of the
212 RMSE. For example if you want a map of each time you can preserve
213 ["time","grid_latitude", "grid_longitude"] or if you want a time series
214 you can preserve ["time"], if you want to collapse to a single value
215 use `None`. The default is `None`.
217 Returns
218 -------
219 scores_cube: iris.cube.Cube
220 A cube containing the RMSE between the base and other cube.
221 """
222 base, other = _sort_cubes_for_verification(cubes)
224 # Copy the coordinates of the input cubes.
225 other_xr = xr.DataArray.from_iris(other)
226 base_xr = xr.DataArray.from_iris(base)
227 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
229 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
230 # apply scores, and then transform it back.
231 scores_cube = xr.DataArray.to_iris(
232 scores.continuous.rmse(
233 other_xr,
234 base_xr,
235 preserve_dims=preserve_dims,
236 )
237 )
239 # If time is aggregated out, attach a scalar time coordinate with bounds
240 # so plotting can display the aggregated period in the title.
241 try:
242 if not scores_cube.coords("time"):
243 base_time = base.coord("time")
244 time_vals = (
245 base_time.bounds.flatten()
246 if base_time.has_bounds()
247 else base_time.points
248 )
249 t_start = float(time_vals[0])
250 t_end = float(time_vals[-1])
251 t_mid = 0.5 * (t_start + t_end)
253 scores_cube.add_aux_coord(
254 iris.coords.AuxCoord(
255 t_mid,
256 standard_name=base_time.standard_name,
257 long_name=base_time.long_name,
258 var_name=base_time.var_name,
259 units=base_time.units,
260 bounds=np.array([t_start, t_end]),
261 attributes=base_time.attributes.copy(),
262 )
263 )
264 except iris.exceptions.CoordinateNotFoundError:
265 pass
267 scores_cube.rename(f"RMSE_of_{base.name()}")
268 # if preserved_coordinates == ["grid_latitude", "grid_longitude"]:
269 # scores_cube.add_aux_coord(time_coord)
270 return scores_cube
273def scores_mae(cubes: CubeList, preserved_coordinates: list[str] | str | None = None):
274 r"""Calculate the Mean Absolute Error (MAE) using scores.
276 Acts as a wrapper around the MAE calculation from ``scores`` ([scoresa]_, [scoresb]_).
278 Parameters
279 ----------
280 cubes: iris.cube.CubeList
281 A CubeList containing exactly two cubes: a base and an "other" model,
282 this can be an analysis and the model.
283 preserved_coordinates: list[str] | str | None, default is None.
284 The coordinates that you wish to preserve in the calculaiton of the
285 MAE. For example if you want a map of each time you can preserve
286 ["time","grid_latitude", "grid_longitude"] or if you want a time series
287 you can preserve ["time"], if you want to collapse to a single value
288 use `None`. The default is `None`.
290 Returns
291 -------
292 scores_cube: iris.cube.Cube
293 A cube containing the MAE between the base and other cube.
294 """
295 base, other = _sort_cubes_for_verification(cubes)
297 # Copy the coordinates of the input cubes.
298 other_xr = xr.DataArray.from_iris(other)
299 base_xr = xr.DataArray.from_iris(base)
300 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
302 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
303 # apply scores, and then transform it back.
304 scores_cube = xr.DataArray.to_iris(
305 scores.continuous.mae(
306 other_xr,
307 base_xr,
308 preserve_dims=preserve_dims,
309 )
310 )
312 # If time is aggregated out, attach a scalar time coordinate with bounds
313 # so plotting can display the aggregated period in the title.
314 try:
315 if not scores_cube.coords("time"): 315 ↛ 340line 315 didn't jump to line 340 because the condition on line 315 was always true
316 base_time = base.coord("time")
317 time_vals = (
318 base_time.bounds.flatten()
319 if base_time.has_bounds()
320 else base_time.points
321 )
322 t_start = float(time_vals[0])
323 t_end = float(time_vals[-1])
324 t_mid = 0.5 * (t_start + t_end)
326 scores_cube.add_aux_coord(
327 iris.coords.AuxCoord(
328 t_mid,
329 standard_name=base_time.standard_name,
330 long_name=base_time.long_name,
331 var_name=base_time.var_name,
332 units=base_time.units,
333 bounds=np.array([t_start, t_end]),
334 attributes=base_time.attributes.copy(),
335 )
336 )
337 except iris.exceptions.CoordinateNotFoundError:
338 pass
340 scores_cube.rename(f"MAE_of_{base.name()}")
341 return scores_cube
344def scores_additive_bias(
345 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
346):
347 r"""Calculate the Additive Bias (Mean Error) using scores.
349 Acts as a wrapper around the ME calculation from ``scores`` ([scoresa]_, [scoresb]_).
351 Parameters
352 ----------
353 cubes: iris.cube.CubeList
354 A CubeList containing exactly two cubes: a base and an "other" model,
355 this can be an analysis and the model.
356 preserved_coordinates: list[str] | str | None, default is None.
357 The coordinates that you wish to preserve in the calculaiton of the
358 ME. For example if you want a map of each time you can preserve
359 ["time","grid_latitude", "grid_longitude"] or if you want a time series
360 you can preserve ["time"], if you want to collapse to a single value
361 use `None`. The default is `None`.
363 Returns
364 -------
365 scores_cube: iris.cube.Cube
366 A cube containing the ME between the base and other cube.
367 """
368 base, other = _sort_cubes_for_verification(cubes)
370 # Copy the coordinates of the input cubes.
371 other_xr = xr.DataArray.from_iris(other)
372 base_xr = xr.DataArray.from_iris(base)
373 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
375 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
376 # apply scores, and then transform it back.
377 scores_cube = xr.DataArray.to_iris(
378 scores.continuous.additive_bias(
379 other_xr,
380 base_xr,
381 preserve_dims=preserve_dims,
382 )
383 )
385 # If time is aggregated out, attach a scalar time coordinate with bounds
386 # so plotting can display the aggregated period in the title.
387 try:
388 if not scores_cube.coords("time"): 388 ↛ 412line 388 didn't jump to line 412 because the condition on line 388 was always true
389 base_time = base.coord("time")
390 time_vals = (
391 base_time.bounds.flatten()
392 if base_time.has_bounds()
393 else base_time.points
394 )
395 t_start = float(time_vals[0])
396 t_end = float(time_vals[-1])
397 t_mid = 0.5 * (t_start + t_end)
399 scores_cube.add_aux_coord(
400 iris.coords.AuxCoord(
401 t_mid,
402 standard_name=base_time.standard_name,
403 long_name=base_time.long_name,
404 var_name=base_time.var_name,
405 units=base_time.units,
406 bounds=np.array([t_start, t_end]),
407 attributes=base_time.attributes.copy(),
408 )
409 )
410 except iris.exceptions.CoordinateNotFoundError:
411 pass
412 scores_cube.rename(f"Additive_Bias_of_{base.name()}")
413 return scores_cube
416def scores_correlation_pearsonr(
417 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
418):
419 r"""Calculate the Pearson's Correlation (PC) coefficient using scores.
421 Acts as a wrapper around the PC calculation from ``scores`` ([scoresa]_, [scoresb]_).
423 Parameters
424 ----------
425 cubes: iris.cube.CubeList
426 A CubeList containing exactly two cubes: a base and an "other" model,
427 this can be an analysis and the model.
428 preserved_coordinates: list[str] | str | None, default is None.
429 The coordinates that you wish to preserve in the calculation of the
430 PC. For example if you want a map of each time you can preserve
431 ["time","grid_latitude", "grid_longitude"] or if you want a time series
432 you can preserve ["time"], if you want to collapse to a single value
433 use `None`. The default is `None`.
435 Returns
436 -------
437 scores_cube: iris.cube.Cube
438 A cube containing the PC between the base and other cube.
439 """
440 base, other = _sort_cubes_for_verification(cubes)
442 # Copy the coordinates of the input cubes.
443 other_xr = xr.DataArray.from_iris(other)
444 base_xr = xr.DataArray.from_iris(base)
445 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
447 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
448 # apply scores, and then transform it back.
449 scores_cube = xr.DataArray.to_iris(
450 scores.continuous.correlation.pearsonr(
451 other_xr,
452 base_xr,
453 preserve_dims=preserve_dims,
454 )
455 )
457 # If time is aggregated out, attach a scalar time coordinate with bounds
458 # so plotting can display the aggregated period in the title.
459 try:
460 if not scores_cube.coords("time"): 460 ↛ 485line 460 didn't jump to line 485 because the condition on line 460 was always true
461 base_time = base.coord("time")
462 time_vals = (
463 base_time.bounds.flatten()
464 if base_time.has_bounds()
465 else base_time.points
466 )
467 t_start = float(time_vals[0])
468 t_end = float(time_vals[-1])
469 t_mid = 0.5 * (t_start + t_end)
471 scores_cube.add_aux_coord(
472 iris.coords.AuxCoord(
473 t_mid,
474 standard_name=base_time.standard_name,
475 long_name=base_time.long_name,
476 var_name=base_time.var_name,
477 units=base_time.units,
478 bounds=np.array([t_start, t_end]),
479 attributes=base_time.attributes.copy(),
480 )
481 )
482 except iris.exceptions.CoordinateNotFoundError:
483 pass
485 scores_cube.rename(f"Pearson_Correlation_of_{base.name()}")
486 return scores_cube
489def scores_crps_for_ensemble(
490 cubes: Cube | CubeList, method: str = "ecdf", control_member: int = 0
491) -> iris.Constraint:
492 r"""Calculate the CRPS for an ensemble.
494 Acts as a wrapper around the crps_for_ensemble from ``scores`` ([scoresa]_, [scoresb]_).
496 Lower CRPS values are better (implies experiment distribution is closer to control distribution/observations),
497 larger values are worse (implies distributions are dissimilar).
498 It is applicable across time and spatial scales as the focus is on the distribution of the values.
499 Default method is ecdf. ecdf is exact value from the empirical distributions,
500 whereas fair produces an approximated value based on a random sample of the underlying distribution.
502 See [CRPS]_ for further information.
504 Parameters
505 ----------
506 cubes: iris.cube.Cube
507 A Cube containing ensembles data
509 Returns
510 -------
511 crps: iris.cube.Cube
512 A cube containing the crps between the ensemble members and the control
513 """
514 if control_member != 0:
515 logger.warning("control member is usual 0")
517 if control_member not in cubes.coords("realization")[0].points:
518 new_control_member = cubes.coords("realization")[0].points[0]
519 logger.warning(
520 f"control member value {control_member} out of bounds, defaulting to control member={new_control_member}"
521 )
522 control_member = new_control_member
524 if cubes.coord("time").shape[0] == 1:
525 raise ValueError("Cube has only one time point.")
527 if cubes.coord("realization").shape[0] < 3:
528 raise ValueError("Cube should have one control member and at least two members")
530 ctrl = cubes.extract(generate_realization_constraint([control_member]))
531 ens_mem = cubes.extract(
532 generate_remove_single_ensemble_member_constraint(control_member)
533 )
535 # Realising the data in advance provides a large speedup
536 _ = ctrl.data
537 _ = ens_mem.data
538 del _
540 ctrl = xr.DataArray.from_iris(ctrl)
541 ens_mem = xr.DataArray.from_iris(ens_mem)
543 crps = xr.DataArray.to_iris(
544 scores.probability.crps_for_ensemble(
545 ens_mem,
546 ctrl,
547 ensemble_member_dim="realization",
548 method=method,
549 preserve_dims="time",
550 )
551 )
553 crps.rename(f"CRPS_of_{cubes[0].name()}")
554 _realization_callback(crps)
555 return crps