Coverage for src/CSET/operators/scoreswrappers.py: 93%
323 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-21 13:35 +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
18import operator
20import iris
21import iris.exceptions
22import numpy as np
23import scores
24import scores.categorical
25import scores.continuous
26import scores.probability
27import xarray as xr
28from iris.cube import Cube, CubeList
29from iris.util import reverse
31from CSET._common import is_increasing
32from CSET.operators._utils import fully_equalise_attributes, get_cube_yxcoordname
33from CSET.operators.constraints import (
34 generate_realization_constraint,
35 generate_remove_single_ensemble_member_constraint,
36)
37from CSET.operators.misc import _extract_common_time_points
38from CSET.operators.read import _realization_callback
39from CSET.operators.regrid import regrid_onto_cube
41logger = logging.getLogger(__name__)
44def _sort_cube_into_base_and_other(cubes):
45 """Sorts cube into base and other models.
47 Parameters
48 ----------
49 cubes: iris.cube.CubeList
50 A CubeList of multiple cubes. One base cube and other model cubes.
52 Returns
53 -------
54 base: iris.cube.Cube
55 The cube from the "analysis" in the same format as the other model.
56 others: iris.cube.CubeList
57 The cube list of containing the cube(s) from the model in the same format as the base model.
59 """
60 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
61 others: CubeList = cubes.extract(
62 iris.Constraint(
63 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
64 )
65 )
67 return base, others
70def _ensure_increasing_pressure_coordinates(cubes):
71 """Ensure the pressure coordinate is increasing.
73 Parameters
74 ----------
75 cubes: iris.cube.CubeList
76 A CubeList of n cubes
78 Returns
79 -------
80 Cubes: iris.cube.CubeList
81 The original cube list but where each cube is ensured to have an increasing pressure coordinate.
82 """
83 for cube in cubes:
84 try:
85 if len(cube.coord("pressure").points) > 2 and not is_increasing(
86 cube.coord("pressure").points
87 ):
88 reverse(cube, "pressure")
90 except iris.exceptions.CoordinateNotFoundError:
91 pass
94def _process_cubes_for_verification(base: Cube, other: Cube):
95 """Prepare cubes ready for verification in scores.
97 Parameters
98 ----------
99 cubes: iris.cube.CubeList
100 A CubeList of exact 2 cubes, one from each model.
102 Returns
103 -------
104 base: iris.cube.Cube
105 The cube from the "analysis" in the same format as the other model.
106 other: iris.cube.Cube
107 The cube from the model in the same format as the base model.
109 Raises
110 ------
111 ValueError: "cubes should contain exactly 2 cubes."
112 If any other number of cubes are present.
114 Notes
115 -----
116 This operator is used for sorting the data into the correct format. It
117 is likely going to need to be refactored out of CSET and perhaps moved into
118 `CSET._utils` given common code between here and `misc.difference`.
119 """
120 # Set cubes into correct format using code from difference operator
122 # Extract just common time points.
123 other_model_name = other.attributes["model_name"]
125 base, other = _extract_common_time_points(base, other)
127 # Get spatial coord names.
128 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
129 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
131 # Ensure cubes to compare are on common differencing grid.
132 # This is triggered if either
133 # i) latitude and longitude shapes are not the same. Note grid points
134 # are not compared directly as these can differ through rounding
135 # errors.
136 # ii) or variables are known to often sit on different grid staggering
137 # in different models (e.g. cell center vs cell edge), as is the case
138 # for UM and LFRic comparisons.
139 # In future greater choice of regridding method might be applied depending
140 # on variable type. Linear regridding can in general be appropriate for smooth
141 # variables. Care should be taken with interpretation of differences
142 # given this dependency on regridding.
143 if (
144 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
145 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
146 ) or (
147 base.long_name
148 in [
149 "eastward_wind_at_10m",
150 "northward_wind_at_10m",
151 "northward_wind_at_cell_centres",
152 "eastward_wind_at_cell_centres",
153 "zonal_wind_at_pressure_levels",
154 "meridional_wind_at_pressure_levels",
155 "potential_vorticity_at_pressure_levels",
156 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
157 ]
158 ):
159 logger.debug("Linear regridding base cube to other grid to compute differences")
160 base = regrid_onto_cube(base, other, method="Linear")
162 # Figure out if we are comparing between UM and LFRic; flip array if so.
163 base_lat_direction = is_increasing(base.coord(base_lat_name).points)
164 other_lat_direction = is_increasing(other.coord(other_lat_name).points)
165 if base_lat_direction != other_lat_direction: 165 ↛ 167line 165 didn't jump to line 167 because the condition on line 165 was never true
166 # Copy base cube for correct coordinate information.
167 other_tmp = base.copy()
168 # Flip the data and place in the copied cube.
169 other_tmp.data = np.flip(
170 other.data, other.coord(other_lat_name).cube_dims(other)
171 )
172 # Use original name and units from the other cube.
173 other_tmp.rename(other.name())
174 other_tmp.units = other.units
175 # Replace the cube.
176 other = other_tmp
178 # Equalise attributes so we can merge.
179 fully_equalise_attributes(CubeList([base, other]))
181 other.attributes["model_name"] = other_model_name
182 logger.debug("Base: %s\nOther: %s", base, other)
184 return base, other
187def _resolve_preserve_dims(
188 cube: Cube,
189 data_array: xr.DataArray,
190 preserved_coordinates: list[str] | str | None,
191) -> list[str] | None:
192 """Resolve preserve coordinates to xarray dimension names.
194 The ``scores`` package expects preserve dimensions to match xarray
195 dimension names. In Iris data, commonly used coordinates such as ``time``
196 may be auxiliary coordinates attached to a differently named dimension
197 (e.g. ``dim0``). This helper maps coordinate names to their underlying
198 dimension names and helps to convert from iris to xarray coordinate dimension names.
199 """
200 if preserved_coordinates is None:
201 return None
203 coord_names = (
204 [preserved_coordinates]
205 if isinstance(preserved_coordinates, str)
206 else preserved_coordinates
207 )
208 preserve_dims: list[str] = []
210 for coord_name in coord_names:
211 # Already an xarray dimension name.
212 if coord_name in data_array.dims:
213 if coord_name not in preserve_dims: 213 ↛ 215line 213 didn't jump to line 215 because the condition on line 213 was always true
214 preserve_dims.append(coord_name)
215 continue
217 # Otherwise, map coordinate name to dimension index/indices.
218 try:
219 dim_indices = cube.coord_dims(coord_name)
220 except iris.exceptions.CoordinateNotFoundError:
221 # Keep original name so scores raises a clear error for unknown keys.
222 if coord_name not in preserve_dims:
223 preserve_dims.append(coord_name)
224 continue
226 for dim_index in dim_indices:
227 dim_name = data_array.dims[dim_index]
228 if dim_name not in preserve_dims:
229 preserve_dims.append(dim_name)
231 return preserve_dims
234def scores_rmse_model_obs(
235 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
236):
237 r"""Calculate the Root Mean Square Error (RMSE) using scores.
239 Acts as a wrapper around the RMSE calculation from ``scores`` ([scoresa]_, [scoresb]_).
240 It is calculated as
242 .. math:: RMSE = \sqrt{\frac{1}{N} \Sigma(forecast - observations)^2}
244 Parameters
245 ----------
246 cubes: iris.cube.CubeList
247 A CubeList containing an observation cube and at least one model cube.
248 preserved_coordinates: list[str] | str | None, default is None.
249 The coordinates that you wish to preserve in the calculaiton of the
250 RMSE. For example if you want a map of each time you can preserve
251 ["time","grid_latitude", "grid_longitude"] or if you want a time series
252 you can preserve ["time"], if you want to collapse to a single value
253 use `None`. The default is `None`.
255 Returns
256 -------
257 scores_cubelist: iris.cube.CubeList
258 A cubelist containing the RMSE between the models and observation cube(s).
259 """
260 rmse_cubes = CubeList()
261 model_list = CubeList()
263 for cb in cubes:
264 if "observed" in cb.long_name:
265 observed = cb
266 else:
267 model_list.append(cb)
269 for model in model_list:
270 input_cubelist = CubeList()
271 input_cubelist.append(observed)
272 input_cubelist.append(model)
273 rmse = scores_rmse(
274 input_cubelist, preserved_coordinates, obs_model_comparison=True
275 )
276 model_name = model.attributes["model_name"]
277 rmse.attributes["model_name"] = model_name
278 rmse_cubes.append(rmse)
280 return rmse_cubes
283def scores_mae_model_obs(
284 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
285):
286 r"""Calculate the Mean Absolute Error (MAE) using scores.
288 Acts as a wrapper around the MAE calculation from ``scores`` ([scoresa]_, [scoresb]_).
290 Parameters
291 ----------
292 cubes: iris.cube.CubeList
293 A CubeList containing an observation cube and at least one model cube.
294 preserved_coordinates: list[str] | str | None, default is None.
295 The coordinates that you wish to preserve in the calculaiton of the
296 MAE. For example if you want a map of each time you can preserve
297 ["time","latitude", "longitude"] or if you want a time series
298 you can preserve ["time"], if you want to collapse to a single value
299 use `None`. The default is `None`.
301 Returns
302 -------
303 scores_cube: iris.cube.Cube
304 A cube containing the MAE between the models and observation cube.
305 """
306 mae_cubes = CubeList()
307 model_list = CubeList()
309 for cb in cubes:
310 if "observed" in cb.long_name:
311 observed = cb
312 else:
313 model_list.append(cb)
315 for model in model_list:
316 input_cubelist = CubeList()
317 input_cubelist.append(observed)
318 input_cubelist.append(model)
319 mae = scores_mae(
320 input_cubelist, preserved_coordinates, obs_model_comparison=True
321 )
322 model_name = model.attributes["model_name"]
323 mae.attributes["model_name"] = model_name
324 mae_cubes.append(mae)
326 return mae_cubes
329def scores_additive_bias_model_obs(
330 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
331):
332 r"""Calculate the Additive Bias (Mean Error) using scores.
334 Acts as a wrapper around the ME calculation from ``scores`` ([scoresa]_, [scoresb]_).
336 Parameters
337 ----------
338 cubes: iris.cube.CubeList
339 A CubeList containing an observation cube and at least one model cube.
340 preserved_coordinates: list[str] | str | None, default is None.
341 The coordinates that you wish to preserve in the calculaiton of the
342 ME. For example if you want a map of each time you can preserve
343 ["time","latitude", "longitude"] or if you want a time series
344 you can preserve ["time"], if you want to collapse to a single value
345 use `None`. The default is `None`.
347 Returns
348 -------
349 scores_cube: iris.cube.CubeList
350 A cube list containing the ME between the models and observation cube.
351 """
352 additive_bias_cubes = CubeList()
353 model_list = CubeList()
355 for cb in cubes:
356 if "observed" in cb.long_name:
357 observed = cb
358 else:
359 model_list.append(cb)
361 for model in model_list:
362 input_cubelist = CubeList()
363 input_cubelist.append(observed)
364 input_cubelist.append(model)
365 additive_bias = scores_additive_bias(
366 input_cubelist, preserved_coordinates, obs_model_comparison=True
367 )
368 model_name = model.attributes["model_name"]
369 additive_bias.attributes["model_name"] = model_name
370 additive_bias_cubes.append(additive_bias)
372 return additive_bias_cubes
375def scores_correlation_pearsonr_model_obs(
376 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
377):
378 r"""Calculate the Pearson's Correlation (PC) coefficient using scores.
380 Acts as a wrapper around the PC calculation from ``scores`` ([scoresa]_, [scoresb]_).
382 Parameters
383 ----------
384 cubes: iris.cube.CubeList
385 A CubeList containing exactly two cubes: a base and an "other" model,
386 this can be an analysis and the model.
387 preserved_coordinates: list[str] | str | None, default is None.
388 The coordinates that you wish to preserve in the calculation of the
389 PC. For example if you want a map of each time you can preserve
390 ["time","latitude", "longitude"] or if you want a time series
391 you can preserve ["time"], if you want to collapse to a single value
392 use `None`. The default is `None`.
394 Returns
395 -------
396 scores_cube: iris.cube.CubeList
397 A cube list containing the PC between the models and observation cube.
398 """
399 pearsonr_cubes = CubeList()
400 model_list = CubeList()
402 for cb in cubes:
403 if "observed" in cb.long_name:
404 observed = cb
405 else:
406 model_list.append(cb)
408 for model in model_list:
409 input_cubelist = CubeList()
410 input_cubelist.append(observed)
411 input_cubelist.append(model)
412 pearsonr = scores_correlation_pearsonr(
413 input_cubelist, preserved_coordinates, obs_model_comparison=True
414 )
415 model_name = model.attributes["model_name"]
416 pearsonr.attributes["model_name"] = model_name
417 pearsonr_cubes.append(pearsonr)
419 return pearsonr_cubes
422def scores_rmse(
423 cubes: CubeList,
424 preserved_coordinates: list[str] | str | None = None,
425 obs_model_comparison: bool = False,
426):
427 r"""Calculate the Root Mean Square Error (RMSE) using scores.
429 Acts as a wrapper around the RMSE calculation from ``scores`` ([scoresa]_, [scoresb]_).
430 It is calculated as
432 .. math:: RMSE = \sqrt{\frac{1}{N} \Sigma(forecast - observations)^2}
434 Parameters
435 ----------
436 cubes: iris.cube.CubeList
437 A CubeList containing exactly two cubes: a base and an "other" model,
438 this can be an analysis and the model.
439 preserved_coordinates: list[str] | str | None, default is None.
440 The coordinates (or xarray dimension names) that you wish to preserve in the calculaiton of the
441 RMSE. For example if you want a map of each time you can preserve
442 ["time","grid_latitude", "grid_longitude"] or if you want a time series
443 you can preserve ["time"], if you want to collapse to a single value
444 use `None`. The default is `None`.
445 obs_model_comparison: bool, default False
446 Set true if doing model-obs comparison.
448 Returns
449 -------
450 scores_cubelist: iris.cube.CubeList
451 A cubelist containing the RMSE between the base and other cube.
452 """
453 scores_cubelist = CubeList()
454 if obs_model_comparison:
455 for cb in cubes:
456 if "observed" in cb.long_name:
457 base = cb
458 else:
459 others = [cb]
460 else:
461 base, others = _sort_cube_into_base_and_other(cubes)
463 for other in others:
464 base, other = _process_cubes_for_verification(base, other)
466 # Copy the coordinates of the input cubes.
467 other_xr = xr.DataArray.from_iris(other)
468 base_xr = xr.DataArray.from_iris(base)
469 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
471 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
472 # apply scores, and then transform it back.
473 scores_cube = xr.DataArray.to_iris(
474 scores.continuous.rmse(
475 other_xr,
476 base_xr,
477 preserve_dims=preserve_dims,
478 )
479 )
481 # If time is aggregated out, attach a scalar time coordinate with bounds
482 # so plotting can display the aggregated period in the title.
483 try:
484 if not scores_cube.coords("time"):
485 base_time = base.coord("time")
486 time_vals = (
487 base_time.bounds.flatten()
488 if base_time.has_bounds()
489 else base_time.points
490 )
491 t_start = float(time_vals[0])
492 t_end = float(time_vals[-1])
493 t_mid = 0.5 * (t_start + t_end)
495 scores_cube.add_aux_coord(
496 iris.coords.AuxCoord(
497 t_mid,
498 standard_name=base_time.standard_name,
499 long_name=base_time.long_name,
500 var_name=base_time.var_name,
501 units=base_time.units,
502 bounds=np.array([t_start, t_end]),
503 attributes=base_time.attributes.copy(),
504 )
505 )
506 except iris.exceptions.CoordinateNotFoundError:
507 pass
509 scores_cube.rename(f"RMSE_of_{base.name()}")
510 scores_cubelist.append(scores_cube)
512 model_name = other.attributes["model_name"]
513 scores_cube.attributes["model_name"] = model_name
515 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
518def scores_mae(
519 cubes: CubeList,
520 preserved_coordinates: list[str] | str | None = None,
521 obs_model_comparison: bool = False,
522):
523 r"""Calculate the Mean Absolute Error (MAE) using scores.
525 Acts as a wrapper around the MAE calculation from ``scores`` ([scoresa]_, [scoresb]_).
527 Parameters
528 ----------
529 cubes: iris.cube.CubeList
530 A CubeList containing exactly two cubes: a base and an "other" model,
531 this can be an analysis and the model.
532 preserved_coordinates: list[str] | str | None, default is None.
533 The coordinates that you wish to preserve in the calculaiton of the
534 MAE. For example if you want a map of each time you can preserve
535 ["time","grid_latitude", "grid_longitude"] or if you want a time series
536 you can preserve ["time"], if you want to collapse to a single value
537 use `None`. The default is `None`.
539 Returns
540 -------
541 scores_cubelist: iris.cube.CubeList
542 A cubelist containing the MAE between the base and other cube(s).
543 """
544 scores_cubelist = CubeList()
545 if obs_model_comparison:
546 for cb in cubes:
547 if "observed" in cb.long_name:
548 base = cb
549 else:
550 others = [cb]
551 else:
552 base, others = _sort_cube_into_base_and_other(cubes)
554 for other in others:
555 base, other = _process_cubes_for_verification(base, other)
557 # Copy the coordinates of the input cubes.
558 other_xr = xr.DataArray.from_iris(other)
559 base_xr = xr.DataArray.from_iris(base)
560 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
562 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
563 # apply scores, and then transform it back.
564 scores_cube = xr.DataArray.to_iris(
565 scores.continuous.mae(
566 other_xr,
567 base_xr,
568 preserve_dims=preserve_dims,
569 )
570 )
572 # If time is aggregated out, attach a scalar time coordinate with bounds
573 # so plotting can display the aggregated period in the title.
574 try:
575 if not scores_cube.coords("time"):
576 base_time = base.coord("time")
577 time_vals = (
578 base_time.bounds.flatten()
579 if base_time.has_bounds()
580 else base_time.points
581 )
582 t_start = float(time_vals[0])
583 t_end = float(time_vals[-1])
584 t_mid = 0.5 * (t_start + t_end)
586 scores_cube.add_aux_coord(
587 iris.coords.AuxCoord(
588 t_mid,
589 standard_name=base_time.standard_name,
590 long_name=base_time.long_name,
591 var_name=base_time.var_name,
592 units=base_time.units,
593 bounds=np.array([t_start, t_end]),
594 attributes=base_time.attributes.copy(),
595 )
596 )
597 except iris.exceptions.CoordinateNotFoundError:
598 pass
600 scores_cube.rename(f"MAE_of_{base.name()}")
601 scores_cubelist.append(scores_cube)
602 model_name = other.attributes["model_name"]
603 scores_cube.attributes["model_name"] = model_name
605 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
608def scores_additive_bias(
609 cubes: CubeList,
610 preserved_coordinates: list[str] | str | None = None,
611 obs_model_comparison: bool = False,
612):
613 r"""Calculate the Additive Bias (Mean Error) using scores.
615 Acts as a wrapper around the ME calculation from ``scores`` ([scoresa]_, [scoresb]_).
617 Parameters
618 ----------
619 cubes: iris.cube.CubeList
620 A CubeList containing exactly two cubes: a base and an "other" model,
621 this can be an analysis and the model.
622 preserved_coordinates: list[str] | str | None, default is None.
623 The coordinates that you wish to preserve in the calculaiton of the
624 ME. For example if you want a map of each time you can preserve
625 ["time","grid_latitude", "grid_longitude"] or if you want a time series
626 you can preserve ["time"], if you want to collapse to a single value
627 use `None`. The default is `None`.
629 Returns
630 -------
631 scores_cubelist: iris.cube.CubeList
632 A cubelist containing the ME between the base and other cube(s).
633 """
634 scores_cubelist = CubeList()
635 if obs_model_comparison:
636 for cb in cubes:
637 if "observed" in cb.long_name:
638 base = cb
639 else:
640 others = [cb]
641 else:
642 base, others = _sort_cube_into_base_and_other(cubes)
644 for other in others:
645 base, other = _process_cubes_for_verification(base, other)
647 # Copy the coordinates of the input cubes.
648 other_xr = xr.DataArray.from_iris(other)
649 base_xr = xr.DataArray.from_iris(base)
650 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
652 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
653 # apply scores, and then transform it back.
654 scores_cube = xr.DataArray.to_iris(
655 scores.continuous.additive_bias(
656 other_xr,
657 base_xr,
658 preserve_dims=preserve_dims,
659 )
660 )
662 # If time is aggregated out, attach a scalar time coordinate with bounds
663 # so plotting can display the aggregated period in the title.
664 try:
665 if not scores_cube.coords("time"):
666 base_time = base.coord("time")
667 time_vals = (
668 base_time.bounds.flatten()
669 if base_time.has_bounds()
670 else base_time.points
671 )
672 t_start = float(time_vals[0])
673 t_end = float(time_vals[-1])
674 t_mid = 0.5 * (t_start + t_end)
676 scores_cube.add_aux_coord(
677 iris.coords.AuxCoord(
678 t_mid,
679 standard_name=base_time.standard_name,
680 long_name=base_time.long_name,
681 var_name=base_time.var_name,
682 units=base_time.units,
683 bounds=np.array([t_start, t_end]),
684 attributes=base_time.attributes.copy(),
685 )
686 )
687 except iris.exceptions.CoordinateNotFoundError:
688 pass
689 scores_cube.rename(f"Additive_Bias_of_{base.name()}")
690 scores_cubelist.append(scores_cube)
691 model_name = other.attributes["model_name"]
692 scores_cube.attributes["model_name"] = model_name
694 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
697def scores_correlation_pearsonr(
698 cubes: CubeList,
699 preserved_coordinates: list[str] | str | None = None,
700 obs_model_comparison: bool = False,
701):
702 r"""Calculate the Pearson's Correlation (PC) coefficient using scores.
704 Acts as a wrapper around the PC calculation from ``scores`` ([scoresa]_, [scoresb]_).
706 Parameters
707 ----------
708 cubes: iris.cube.CubeList
709 A CubeList containing exactly two cubes: a base and an "other" model,
710 this can be an analysis and the model.
711 preserved_coordinates: list[str] | str | None, default is None.
712 The coordinates that you wish to preserve in the calculation of the
713 PC. For example if you want a map of each time you can preserve
714 ["time","grid_latitude", "grid_longitude"] or if you want a time series
715 you can preserve ["time"], if you want to collapse to a single value
716 use `None`. The default is `None`.
718 Returns
719 -------
720 scores_cubelist: iris.cube.CubeList
721 A cubelist containing the PC between the base and other cube(s).
722 """
723 scores_cubelist = CubeList()
724 if obs_model_comparison:
725 for cb in cubes:
726 if "observed" in cb.long_name:
727 base = cb
728 else:
729 others = [cb]
730 else:
731 base, others = _sort_cube_into_base_and_other(cubes)
733 for other in others:
734 base, other = _process_cubes_for_verification(base, other)
736 # Copy the coordinates of the input cubes.
737 other_xr = xr.DataArray.from_iris(other)
738 base_xr = xr.DataArray.from_iris(base)
739 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
741 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
742 # apply scores, and then transform it back.
743 scores_cube = xr.DataArray.to_iris(
744 scores.continuous.correlation.pearsonr(
745 other_xr,
746 base_xr,
747 preserve_dims=preserve_dims,
748 )
749 )
751 # If time is aggregated out, attach a scalar time coordinate with bounds
752 # so plotting can display the aggregated period in the title.
753 try:
754 if not scores_cube.coords("time"):
755 base_time = base.coord("time")
756 time_vals = (
757 base_time.bounds.flatten()
758 if base_time.has_bounds()
759 else base_time.points
760 )
761 t_start = float(time_vals[0])
762 t_end = float(time_vals[-1])
763 t_mid = 0.5 * (t_start + t_end)
765 scores_cube.add_aux_coord(
766 iris.coords.AuxCoord(
767 t_mid,
768 standard_name=base_time.standard_name,
769 long_name=base_time.long_name,
770 var_name=base_time.var_name,
771 units=base_time.units,
772 bounds=np.array([t_start, t_end]),
773 attributes=base_time.attributes.copy(),
774 )
775 )
776 except iris.exceptions.CoordinateNotFoundError:
777 pass
779 scores_cube.rename(f"Pearson_Correlation_of_{base.name()}")
780 scores_cubelist.append(scores_cube)
781 model_name = other.attributes["model_name"]
782 scores_cube.attributes["model_name"] = model_name
783 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
786def scores_crps_for_ensemble(
787 cubes: Cube | CubeList, method: str = "ecdf", control_member: int = 0
788) -> iris.Constraint:
789 r"""Calculate the CRPS for an ensemble.
791 Acts as a wrapper around the crps_for_ensemble from ``scores`` ([scoresa]_, [scoresb]_).
793 Lower CRPS values are better (implies experiment distribution is closer to control distribution/observations),
794 larger values are worse (implies distributions are dissimilar).
795 It is applicable across time and spatial scales as the focus is on the distribution of the values.
796 Default method is ecdf. ecdf is exact value from the empirical distributions,
797 whereas fair produces an approximated value based on a random sample of the underlying distribution.
799 See [CRPS]_ for further information.
801 Parameters
802 ----------
803 cubes: iris.cube.Cube
804 A Cube containing ensembles data
806 Returns
807 -------
808 crps: iris.cube.Cube
809 A cube containing the crps between the ensemble members and the control
810 """
811 if control_member != 0:
812 logger.warning("control member is usual 0")
814 if control_member not in cubes.coords("realization")[0].points:
815 new_control_member = cubes.coords("realization")[0].points[0]
816 logger.warning(
817 f"control member value {control_member} out of bounds, defaulting to control member={new_control_member}"
818 )
819 control_member = new_control_member
821 if cubes.coord("time").shape[0] == 1:
822 raise ValueError("Cube has only one time point.")
824 if cubes.coord("realization").shape[0] < 3:
825 raise ValueError("Cube should have one control member and at least two members")
827 ctrl = cubes.extract(generate_realization_constraint([control_member]))
828 ens_mem = cubes.extract(
829 generate_remove_single_ensemble_member_constraint(control_member)
830 )
832 # Realising the data in advance provides a large speedup
833 _ = ctrl.data
834 _ = ens_mem.data
835 del _
837 ctrl = xr.DataArray.from_iris(ctrl)
838 ens_mem = xr.DataArray.from_iris(ens_mem)
840 crps = xr.DataArray.to_iris(
841 scores.probability.crps_for_ensemble(
842 ens_mem,
843 ctrl,
844 ensemble_member_dim="realization",
845 method=method,
846 preserve_dims="time",
847 )
848 )
850 crps.rename(f"CRPS_of_{cubes[0].name()}")
851 _realization_callback(crps)
852 return crps
855def scores_pod_model_obs(
856 cubes: CubeList,
857 preserved_coordinates: list[str] | str | None,
858 threshold: str,
859 op_func: str,
860):
861 r"""
862 Compute the Probability of Detection (POD) score using Scores ([scoresa]_ [scoresb]_).
864 Parameters
865 ----------
866 cubes: iris.cube.CubeList
867 An iris cubelist containing model(s) and an observation cube.
868 preserved_coordinates: list | str | None
869 An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
870 then preserving coordinate 'time' will produce a probability of detection score for each timeslice (shape time). If None,
871 then it will return a single value score for all times/point locations.
872 threshold: str
873 A str containing the threshold to use to generate the binary masks, which subsequently gets turned to a float (but passed as str around the recipe templating).
874 op_func: str
875 A string either containing 'lt' for less than or 'gt for greater than, to determine how the threshold is applied to the data
876 to generate the mask.
878 Returns
879 -------
880 cube: iris.cube
881 An iris cube, containing the probability of detection score for further plotting.
883 Notes
884 -----
885 The probability of detection calculates the proportion of observed events that meet a threshold that were correctly forecast by the model.
886 For example, if threshold is 290K and op_func is gt (greater than), and at some station a temperature was recorded as 292K and the model produced
887 295k, that would be a positive hit. It does not take into account how far above/below a threshold a model forecasts.
889 It is calculated as .. math:: POD = \frac{true positives}{true positives + false negatives}
891 It is equivalent to the hit rate. Note if there are no events that meet the threshold in model and observations, a POD of zero is returned.
893 POD produces a range of 0 to 1, where 1 is a perfect score.
894 """
895 # Split out model(s) and obs
896 models = CubeList()
897 for c in cubes:
898 if "observed" in c.long_name:
899 observed = c
900 else:
901 models.append(c)
903 # Setup cubelist to store results
904 scores_results = iris.cube.CubeList()
906 # Setup operators greater than, less than.
907 ops = {
908 "gt": operator.gt,
909 "lt": operator.lt,
910 }
912 try:
913 op = ops[op_func]
914 except KeyError as err:
915 raise ValueError(f"Operator {op_func} not supported.") from err
917 for model in models:
918 # Convert obs cubes to xarray and resolve preserved dimensions.
919 other_xr = xr.DataArray.from_iris(model)
920 base_xr = xr.DataArray.from_iris(observed)
921 preserve_dims = _resolve_preserve_dims(
922 observed, other_xr, preserved_coordinates
923 )
925 # Create event operator object using threshold and operator direction.
926 event_operator = scores.categorical.ThresholdEventOperator(
927 default_event_threshold=float(threshold), default_op_fn=op
928 )
930 # Generate binary fields using the event operator.
931 forecast_binary, observed_binary = event_operator.make_event_tables(
932 other_xr, base_xr
933 )
935 # Create binary contigency manager, as per Scores API, using transform to preserve preserve_dims
936 contingency_manager = scores.categorical.BinaryContingencyManager(
937 forecast_binary, observed_binary
938 ).transform(preserve_dims=preserve_dims)
940 # Get POD from the contigency manager, and convert back to an iris cube.
941 scores_cube = xr.DataArray.to_iris(
942 contingency_manager.probability_of_detection()
943 )
945 # Rename cube so it plots correctly alongside correcting cube units.
946 scores_cube.rename(
947 f"Probability_Of_Detection_{op_func}_{threshold}_{observed.name()}"
948 )
949 scores_cube.units = "1"
950 scores_cube.attributes["model_name"] = model.attributes["model_name"]
952 scores_results.append(scores_cube)
954 return scores_results
957def scores_ets_model_obs(
958 cubes: CubeList,
959 preserved_coordinates: list[str] | str | None,
960 threshold: str,
961 op_func: str,
962):
963 r"""
964 Compute the Equitable Threat Score (ETS) score using Scores ([scoresa]_ [scoresb]_).
966 Parameters
967 ----------
968 cubes: iris.cube.CubeList
969 An iris cubelist containing model(s) and an observation cube.
970 preserved_coordinates: list | str | None
971 An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
972 then preserving coordinate 'time' will produce the equitable threat score for each timeslice (shape time). If None,
973 then it will return a single value score for all times/point locations.
974 threshold: str
975 A str containing the threshold to use to generate the binary masks, which subsequently gets turned to a float (but passed as str around the recipe templating).
976 op_func: str
977 A string either containing 'lt' for less than or 'gt for greater than, to determine how the threshold is applied to the data
978 to generate the mask.
980 Returns
981 -------
982 cube: iris.cube
983 An iris cube, containing the probability of detection score for further plotting.
985 Notes
986 -----
987 The Equitable Threat Score (ETS) evaluates the accuracy of forecasts for events that meet a specified threshold,
988 hile accounting for correct forecasts that could occur purely by chance. Unlike the Probability of Detection (POD),
989 ETS considers hits, misses, and false alarms, providing a more balanced assessment of forecast skill.
991 For example, if the threshold is 290 K and op_func is gt (greater than), an observation of 292 K and a forecast of 295 K
992 would be counted as a hit. ETS adjusts the total number of hits by removing the number of hits expected due to random chance.
994 It is calculated as:
996 .. math::
998 ETS = \frac{hits - hits_{random}}
999 {hits + misses + false\ alarms - hits_{random}}
1001 where
1003 hits_{random} = \frac{(hits + misses)(hits + false\ alarms)}{total count}
1005 ETS ranges from -1/3 to 1, where 1 indicates a perfect forecast, 0 indicates no skill beyond random chance, and negative values indicate worse than
1006 random chance.
1007 """
1008 # Split out model(s) and obs
1009 models = CubeList()
1010 for c in cubes:
1011 if "observed" in c.long_name:
1012 observed = c
1013 else:
1014 models.append(c)
1016 # Setup cubelist to store results
1017 scores_results = iris.cube.CubeList()
1019 # Setup operators greater than, less than.
1020 ops = {
1021 "gt": operator.gt,
1022 "lt": operator.lt,
1023 }
1025 try:
1026 op = ops[op_func]
1027 except KeyError as err:
1028 raise ValueError(f"Operator {op_func} not supported.") from err
1030 for model in models:
1031 # Convert obs cubes to xarray and resolve preserved dimensions.
1032 other_xr = xr.DataArray.from_iris(model)
1033 base_xr = xr.DataArray.from_iris(observed)
1034 preserve_dims = _resolve_preserve_dims(
1035 observed, other_xr, preserved_coordinates
1036 )
1038 # Create event operator object using threshold and operator direction.
1039 event_operator = scores.categorical.ThresholdEventOperator(
1040 default_event_threshold=float(threshold), default_op_fn=op
1041 )
1043 # Generate binary fields using the event operator.
1044 forecast_binary, observed_binary = event_operator.make_event_tables(
1045 other_xr, base_xr
1046 )
1048 # Create binary contigency manager, as per Scores API, using transform to preserve preserve_dims
1049 contingency_manager = scores.categorical.BinaryContingencyManager(
1050 forecast_binary, observed_binary
1051 ).transform(preserve_dims=preserve_dims)
1053 # Get ETS from the contigency manager, and convert back to an iris cube.
1054 scores_cube = xr.DataArray.to_iris(contingency_manager.equitable_threat_score())
1056 # Rename cube so it plots correctly alongside correcting cube units.
1057 scores_cube.rename(
1058 f"Equitable_Threat_Score_{op_func}_{threshold}_{observed.name()}"
1059 )
1060 scores_cube.units = "1"
1061 scores_cube.attributes["model_name"] = model.attributes["model_name"]
1063 scores_results.append(scores_cube)
1065 return scores_results