Coverage for src/CSET/operators/scoreswrappers.py: 90%
181 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-13 14:57 +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:
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_model_obs(
198 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
199):
200 r"""Calculate the Root Mean Square Error (RMSE) using scores.
202 Acts as a wrapper around the RMSE calculation from ``scores`` ([scoresa]_, [scoresb]_).
203 It is calculated as
205 .. math:: RMSE = \sqrt{\frac{1}{N} \Sigma(forecast - observations)^2}
207 Parameters
208 ----------
209 cubes: iris.cube.CubeList
210 A CubeList containing an observation cube and at least one model cube.
211 preserved_coordinates: list[str] | str | None, default is None.
212 The coordinates that you wish to preserve in the calculaiton of the
213 RMSE. For example if you want a map of each time you can preserve
214 ["time","grid_latitude", "grid_longitude"] or if you want a time series
215 you can preserve ["time"], if you want to collapse to a single value
216 use `None`. The default is `None`.
218 Returns
219 -------
220 scores_cube: iris.cube.Cube
221 A cube containing the RMSE between the models and observation cube.
222 """
223 rmse_cubes = CubeList()
224 model_list = CubeList()
226 for cb in cubes:
227 if "observed" in cb.long_name:
228 observed = cb
229 else:
230 model_list.append(cb)
232 for model in model_list:
233 input_cubelist = CubeList()
234 input_cubelist.append(observed)
235 input_cubelist.append(model)
236 rmse = scores_rmse(
237 input_cubelist, preserved_coordinates, obs_model_comparison=True
238 )
239 model_name = model.attributes["model_name"]
240 rmse.attributes["model_name"] = model_name
241 rmse_cubes.append(rmse)
243 return rmse_cubes
246def scores_rmse(
247 cubes: CubeList,
248 preserved_coordinates: list[str] | str | None = None,
249 obs_model_comparison: bool = False,
250):
251 r"""Calculate the Root Mean Square Error (RMSE) using scores.
253 Acts as a wrapper around the RMSE calculation from ``scores`` ([scoresa]_, [scoresb]_).
254 It is calculated as
256 .. math:: RMSE = \sqrt{\frac{1}{N} \Sigma(forecast - observations)^2}
258 Parameters
259 ----------
260 cubes: iris.cube.CubeList
261 A CubeList containing exactly two cubes: a base and an "other" model,
262 this can be an analysis and the model.
263 preserved_coordinates: list[str] | str | None, default is None.
264 The coordinates (or xarray dimension names) that you wish to preserve in the calculaiton of the
265 RMSE. For example if you want a map of each time you can preserve
266 ["time","grid_latitude", "grid_longitude"] or if you want a time series
267 you can preserve ["time"], if you want to collapse to a single value
268 use `None`. The default is `None`.
270 Returns
271 -------
272 scores_cube: iris.cube.Cube
273 A cube containing the RMSE between the base and other cube.
274 """
275 if obs_model_comparison:
276 for cb in cubes:
277 if "observed" in cb.long_name:
278 base = cb
279 else:
280 other = cb
281 else:
282 base, other = _sort_cubes_for_verification(cubes)
284 # Copy the coordinates of the input cubes.
285 other_xr = xr.DataArray.from_iris(other)
286 base_xr = xr.DataArray.from_iris(base)
287 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
289 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
290 # apply scores, and then transform it back.
291 scores_cube = xr.DataArray.to_iris(
292 scores.continuous.rmse(
293 other_xr,
294 base_xr,
295 preserve_dims=preserve_dims,
296 )
297 )
299 # If time is aggregated out, attach a scalar time coordinate with bounds
300 # so plotting can display the aggregated period in the title.
301 try:
302 if not scores_cube.coords("time"):
303 base_time = base.coord("time")
304 time_vals = (
305 base_time.bounds.flatten()
306 if base_time.has_bounds()
307 else base_time.points
308 )
309 t_start = float(time_vals[0])
310 t_end = float(time_vals[-1])
311 t_mid = 0.5 * (t_start + t_end)
313 scores_cube.add_aux_coord(
314 iris.coords.AuxCoord(
315 t_mid,
316 standard_name=base_time.standard_name,
317 long_name=base_time.long_name,
318 var_name=base_time.var_name,
319 units=base_time.units,
320 bounds=np.array([t_start, t_end]),
321 attributes=base_time.attributes.copy(),
322 )
323 )
324 except iris.exceptions.CoordinateNotFoundError:
325 pass
327 scores_cube.rename(f"RMSE_of_{base.name()}")
328 # if preserved_coordinates == ["grid_latitude", "grid_longitude"]:
329 # scores_cube.add_aux_coord(time_coord)
330 return scores_cube
333def scores_mae(cubes: CubeList, preserved_coordinates: list[str] | str | None = None):
334 r"""Calculate the Mean Absolute Error (MAE) using scores.
336 Acts as a wrapper around the MAE calculation from ``scores`` ([scoresa]_, [scoresb]_).
338 Parameters
339 ----------
340 cubes: iris.cube.CubeList
341 A CubeList containing exactly two cubes: a base and an "other" model,
342 this can be an analysis and the model.
343 preserved_coordinates: list[str] | str | None, default is None.
344 The coordinates that you wish to preserve in the calculaiton of the
345 MAE. For example if you want a map of each time you can preserve
346 ["time","grid_latitude", "grid_longitude"] or if you want a time series
347 you can preserve ["time"], if you want to collapse to a single value
348 use `None`. The default is `None`.
350 Returns
351 -------
352 scores_cube: iris.cube.Cube
353 A cube containing the MAE between the base and other cube.
354 """
355 base, other = _sort_cubes_for_verification(cubes)
357 # Copy the coordinates of the input cubes.
358 other_xr = xr.DataArray.from_iris(other)
359 base_xr = xr.DataArray.from_iris(base)
360 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
362 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
363 # apply scores, and then transform it back.
364 scores_cube = xr.DataArray.to_iris(
365 scores.continuous.mae(
366 other_xr,
367 base_xr,
368 preserve_dims=preserve_dims,
369 )
370 )
372 # If time is aggregated out, attach a scalar time coordinate with bounds
373 # so plotting can display the aggregated period in the title.
374 try:
375 if not scores_cube.coords("time"): 375 ↛ 400line 375 didn't jump to line 400 because the condition on line 375 was always true
376 base_time = base.coord("time")
377 time_vals = (
378 base_time.bounds.flatten()
379 if base_time.has_bounds()
380 else base_time.points
381 )
382 t_start = float(time_vals[0])
383 t_end = float(time_vals[-1])
384 t_mid = 0.5 * (t_start + t_end)
386 scores_cube.add_aux_coord(
387 iris.coords.AuxCoord(
388 t_mid,
389 standard_name=base_time.standard_name,
390 long_name=base_time.long_name,
391 var_name=base_time.var_name,
392 units=base_time.units,
393 bounds=np.array([t_start, t_end]),
394 attributes=base_time.attributes.copy(),
395 )
396 )
397 except iris.exceptions.CoordinateNotFoundError:
398 pass
400 scores_cube.rename(f"MAE_of_{base.name()}")
401 return scores_cube
404def scores_additive_bias(
405 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
406):
407 r"""Calculate the Additive Bias (Mean Error) using scores.
409 Acts as a wrapper around the ME calculation from ``scores`` ([scoresa]_, [scoresb]_).
411 Parameters
412 ----------
413 cubes: iris.cube.CubeList
414 A CubeList containing exactly two cubes: a base and an "other" model,
415 this can be an analysis and the model.
416 preserved_coordinates: list[str] | str | None, default is None.
417 The coordinates that you wish to preserve in the calculaiton of the
418 ME. For example if you want a map of each time you can preserve
419 ["time","grid_latitude", "grid_longitude"] or if you want a time series
420 you can preserve ["time"], if you want to collapse to a single value
421 use `None`. The default is `None`.
423 Returns
424 -------
425 scores_cube: iris.cube.Cube
426 A cube containing the ME between the base and other cube.
427 """
428 base, other = _sort_cubes_for_verification(cubes)
430 # Copy the coordinates of the input cubes.
431 other_xr = xr.DataArray.from_iris(other)
432 base_xr = xr.DataArray.from_iris(base)
433 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
435 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
436 # apply scores, and then transform it back.
437 scores_cube = xr.DataArray.to_iris(
438 scores.continuous.additive_bias(
439 other_xr,
440 base_xr,
441 preserve_dims=preserve_dims,
442 )
443 )
445 # If time is aggregated out, attach a scalar time coordinate with bounds
446 # so plotting can display the aggregated period in the title.
447 try:
448 if not scores_cube.coords("time"): 448 ↛ 472line 448 didn't jump to line 472 because the condition on line 448 was always true
449 base_time = base.coord("time")
450 time_vals = (
451 base_time.bounds.flatten()
452 if base_time.has_bounds()
453 else base_time.points
454 )
455 t_start = float(time_vals[0])
456 t_end = float(time_vals[-1])
457 t_mid = 0.5 * (t_start + t_end)
459 scores_cube.add_aux_coord(
460 iris.coords.AuxCoord(
461 t_mid,
462 standard_name=base_time.standard_name,
463 long_name=base_time.long_name,
464 var_name=base_time.var_name,
465 units=base_time.units,
466 bounds=np.array([t_start, t_end]),
467 attributes=base_time.attributes.copy(),
468 )
469 )
470 except iris.exceptions.CoordinateNotFoundError:
471 pass
472 scores_cube.rename(f"Additive_Bias_of_{base.name()}")
473 return scores_cube
476def scores_correlation_pearsonr(
477 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
478):
479 r"""Calculate the Pearson's Correlation (PC) coefficient using scores.
481 Acts as a wrapper around the PC calculation from ``scores`` ([scoresa]_, [scoresb]_).
483 Parameters
484 ----------
485 cubes: iris.cube.CubeList
486 A CubeList containing exactly two cubes: a base and an "other" model,
487 this can be an analysis and the model.
488 preserved_coordinates: list[str] | str | None, default is None.
489 The coordinates that you wish to preserve in the calculation of the
490 PC. For example if you want a map of each time you can preserve
491 ["time","grid_latitude", "grid_longitude"] or if you want a time series
492 you can preserve ["time"], if you want to collapse to a single value
493 use `None`. The default is `None`.
495 Returns
496 -------
497 scores_cube: iris.cube.Cube
498 A cube containing the PC between the base and other cube.
499 """
500 base, other = _sort_cubes_for_verification(cubes)
502 # Copy the coordinates of the input cubes.
503 other_xr = xr.DataArray.from_iris(other)
504 base_xr = xr.DataArray.from_iris(base)
505 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
507 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
508 # apply scores, and then transform it back.
509 scores_cube = xr.DataArray.to_iris(
510 scores.continuous.correlation.pearsonr(
511 other_xr,
512 base_xr,
513 preserve_dims=preserve_dims,
514 )
515 )
517 # If time is aggregated out, attach a scalar time coordinate with bounds
518 # so plotting can display the aggregated period in the title.
519 try:
520 if not scores_cube.coords("time"): 520 ↛ 545line 520 didn't jump to line 545 because the condition on line 520 was always true
521 base_time = base.coord("time")
522 time_vals = (
523 base_time.bounds.flatten()
524 if base_time.has_bounds()
525 else base_time.points
526 )
527 t_start = float(time_vals[0])
528 t_end = float(time_vals[-1])
529 t_mid = 0.5 * (t_start + t_end)
531 scores_cube.add_aux_coord(
532 iris.coords.AuxCoord(
533 t_mid,
534 standard_name=base_time.standard_name,
535 long_name=base_time.long_name,
536 var_name=base_time.var_name,
537 units=base_time.units,
538 bounds=np.array([t_start, t_end]),
539 attributes=base_time.attributes.copy(),
540 )
541 )
542 except iris.exceptions.CoordinateNotFoundError:
543 pass
545 scores_cube.rename(f"Pearson_Correlation_of_{base.name()}")
546 return scores_cube
549def scores_crps_for_ensemble(
550 cubes: Cube | CubeList, method: str = "ecdf", control_member: int = 0
551) -> iris.Constraint:
552 r"""Calculate the CRPS for an ensemble.
554 Acts as a wrapper around the crps_for_ensemble from ``scores`` ([scoresa]_, [scoresb]_).
556 Lower CRPS values are better (implies experiment distribution is closer to control distribution/observations),
557 larger values are worse (implies distributions are dissimilar).
558 It is applicable across time and spatial scales as the focus is on the distribution of the values.
559 Default method is ecdf. ecdf is exact value from the empirical distributions,
560 whereas fair produces an approximated value based on a random sample of the underlying distribution.
562 See [CRPS]_ for further information.
564 Parameters
565 ----------
566 cubes: iris.cube.Cube
567 A Cube containing ensembles data
569 Returns
570 -------
571 crps: iris.cube.Cube
572 A cube containing the crps between the ensemble members and the control
573 """
574 if control_member != 0:
575 logger.warning("control member is usual 0")
577 if control_member not in cubes.coords("realization")[0].points:
578 new_control_member = cubes.coords("realization")[0].points[0]
579 logger.warning(
580 f"control member value {control_member} out of bounds, defaulting to control member={new_control_member}"
581 )
582 control_member = new_control_member
584 if cubes.coord("time").shape[0] == 1:
585 raise ValueError("Cube has only one time point.")
587 if cubes.coord("realization").shape[0] < 3:
588 raise ValueError("Cube should have one control member and at least two members")
590 ctrl = cubes.extract(generate_realization_constraint([control_member]))
591 ens_mem = cubes.extract(
592 generate_remove_single_ensemble_member_constraint(control_member)
593 )
595 # Realising the data in advance provides a large speedup
596 _ = ctrl.data
597 _ = ens_mem.data
598 del _
600 ctrl = xr.DataArray.from_iris(ctrl)
601 ens_mem = xr.DataArray.from_iris(ens_mem)
603 crps = xr.DataArray.to_iris(
604 scores.probability.crps_for_ensemble(
605 ens_mem,
606 ctrl,
607 ensemble_member_dim="realization",
608 method=method,
609 preserve_dims="time",
610 )
611 )
613 crps.rename(f"CRPS_of_{cubes[0].name()}")
614 _realization_callback(crps)
615 return crps