Coverage for src/CSET/operators/scoreswrappers.py: 92%
302 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-21 13:10 +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_additive_bias_model_obs(
284 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
285):
286 r"""Calculate the Additive Bias (Mean Error) using scores.
288 Acts as a wrapper around the ME 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 ME. 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.CubeList
304 A cube list containing the ME between the models and observation cube.
305 """
306 additive_bias_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 additive_bias = scores_additive_bias(
320 input_cubelist, preserved_coordinates, obs_model_comparison=True
321 )
322 model_name = model.attributes["model_name"]
323 additive_bias.attributes["model_name"] = model_name
324 additive_bias_cubes.append(additive_bias)
326 return additive_bias_cubes
329def scores_correlation_pearsonr_model_obs(
330 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
331):
332 r"""Calculate the Pearson's Correlation (PC) coefficient using scores.
334 Acts as a wrapper around the PC calculation from ``scores`` ([scoresa]_, [scoresb]_).
336 Parameters
337 ----------
338 cubes: iris.cube.CubeList
339 A CubeList containing exactly two cubes: a base and an "other" model,
340 this can be an analysis and the model.
341 preserved_coordinates: list[str] | str | None, default is None.
342 The coordinates that you wish to preserve in the calculation of the
343 PC. For example if you want a map of each time you can preserve
344 ["time","latitude", "longitude"] or if you want a time series
345 you can preserve ["time"], if you want to collapse to a single value
346 use `None`. The default is `None`.
348 Returns
349 -------
350 scores_cube: iris.cube.CubeList
351 A cube list containing the PC between the models and observation cube.
352 """
353 pearsonr_cubes = CubeList()
354 model_list = CubeList()
356 for cb in cubes:
357 if "observed" in cb.long_name:
358 observed = cb
359 else:
360 model_list.append(cb)
362 for model in model_list:
363 input_cubelist = CubeList()
364 input_cubelist.append(observed)
365 input_cubelist.append(model)
366 pearsonr = scores_correlation_pearsonr(
367 input_cubelist, preserved_coordinates, obs_model_comparison=True
368 )
369 model_name = model.attributes["model_name"]
370 pearsonr.attributes["model_name"] = model_name
371 pearsonr_cubes.append(pearsonr)
373 return pearsonr_cubes
376def scores_rmse(
377 cubes: CubeList,
378 preserved_coordinates: list[str] | str | None = None,
379 obs_model_comparison: bool = False,
380):
381 r"""Calculate the Root Mean Square Error (RMSE) using scores.
383 Acts as a wrapper around the RMSE calculation from ``scores`` ([scoresa]_, [scoresb]_).
384 It is calculated as
386 .. math:: RMSE = \sqrt{\frac{1}{N} \Sigma(forecast - observations)^2}
388 Parameters
389 ----------
390 cubes: iris.cube.CubeList
391 A CubeList containing exactly two cubes: a base and an "other" model,
392 this can be an analysis and the model.
393 preserved_coordinates: list[str] | str | None, default is None.
394 The coordinates (or xarray dimension names) that you wish to preserve in the calculaiton of the
395 RMSE. For example if you want a map of each time you can preserve
396 ["time","grid_latitude", "grid_longitude"] or if you want a time series
397 you can preserve ["time"], if you want to collapse to a single value
398 use `None`. The default is `None`.
399 obs_model_comparison: bool, default False
400 Set true if doing model-obs comparison.
402 Returns
403 -------
404 scores_cubelist: iris.cube.CubeList
405 A cubelist containing the RMSE between the base and other cube.
406 """
407 scores_cubelist = CubeList()
408 if obs_model_comparison:
409 for cb in cubes:
410 if "observed" in cb.long_name:
411 base = cb
412 else:
413 others = [cb]
414 else:
415 base, others = _sort_cube_into_base_and_other(cubes)
417 for other in others:
418 base, other = _process_cubes_for_verification(base, other)
420 # Copy the coordinates of the input cubes.
421 other_xr = xr.DataArray.from_iris(other)
422 base_xr = xr.DataArray.from_iris(base)
423 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
425 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
426 # apply scores, and then transform it back.
427 scores_cube = xr.DataArray.to_iris(
428 scores.continuous.rmse(
429 other_xr,
430 base_xr,
431 preserve_dims=preserve_dims,
432 )
433 )
435 # If time is aggregated out, attach a scalar time coordinate with bounds
436 # so plotting can display the aggregated period in the title.
437 try:
438 if not scores_cube.coords("time"):
439 base_time = base.coord("time")
440 time_vals = (
441 base_time.bounds.flatten()
442 if base_time.has_bounds()
443 else base_time.points
444 )
445 t_start = float(time_vals[0])
446 t_end = float(time_vals[-1])
447 t_mid = 0.5 * (t_start + t_end)
449 scores_cube.add_aux_coord(
450 iris.coords.AuxCoord(
451 t_mid,
452 standard_name=base_time.standard_name,
453 long_name=base_time.long_name,
454 var_name=base_time.var_name,
455 units=base_time.units,
456 bounds=np.array([t_start, t_end]),
457 attributes=base_time.attributes.copy(),
458 )
459 )
460 except iris.exceptions.CoordinateNotFoundError:
461 pass
463 scores_cube.rename(f"RMSE_of_{base.name()}")
464 scores_cubelist.append(scores_cube)
466 model_name = other.attributes["model_name"]
467 scores_cube.attributes["model_name"] = model_name
469 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
472def scores_mae(cubes: CubeList, preserved_coordinates: list[str] | str | None = None):
473 r"""Calculate the Mean Absolute Error (MAE) using scores.
475 Acts as a wrapper around the MAE calculation from ``scores`` ([scoresa]_, [scoresb]_).
477 Parameters
478 ----------
479 cubes: iris.cube.CubeList
480 A CubeList containing exactly two cubes: a base and an "other" model,
481 this can be an analysis and the model.
482 preserved_coordinates: list[str] | str | None, default is None.
483 The coordinates that you wish to preserve in the calculaiton of the
484 MAE. For example if you want a map of each time you can preserve
485 ["time","grid_latitude", "grid_longitude"] or if you want a time series
486 you can preserve ["time"], if you want to collapse to a single value
487 use `None`. The default is `None`.
489 Returns
490 -------
491 scores_cubelist: iris.cube.CubeList
492 A cubelist containing the MAE between the base and other cube(s).
493 """
494 base, others = _sort_cube_into_base_and_other(cubes)
495 scores_cubelist = CubeList()
496 for other in others:
497 base, other = _process_cubes_for_verification(base, other)
499 # Copy the coordinates of the input cubes.
500 other_xr = xr.DataArray.from_iris(other)
501 base_xr = xr.DataArray.from_iris(base)
502 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
504 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
505 # apply scores, and then transform it back.
506 scores_cube = xr.DataArray.to_iris(
507 scores.continuous.mae(
508 other_xr,
509 base_xr,
510 preserve_dims=preserve_dims,
511 )
512 )
514 # If time is aggregated out, attach a scalar time coordinate with bounds
515 # so plotting can display the aggregated period in the title.
516 try:
517 if not scores_cube.coords("time"): 517 ↛ 542line 517 didn't jump to line 542 because the condition on line 517 was always true
518 base_time = base.coord("time")
519 time_vals = (
520 base_time.bounds.flatten()
521 if base_time.has_bounds()
522 else base_time.points
523 )
524 t_start = float(time_vals[0])
525 t_end = float(time_vals[-1])
526 t_mid = 0.5 * (t_start + t_end)
528 scores_cube.add_aux_coord(
529 iris.coords.AuxCoord(
530 t_mid,
531 standard_name=base_time.standard_name,
532 long_name=base_time.long_name,
533 var_name=base_time.var_name,
534 units=base_time.units,
535 bounds=np.array([t_start, t_end]),
536 attributes=base_time.attributes.copy(),
537 )
538 )
539 except iris.exceptions.CoordinateNotFoundError:
540 pass
542 scores_cube.rename(f"MAE_of_{base.name()}")
543 scores_cubelist.append(scores_cube)
544 model_name = other.attributes["model_name"]
545 scores_cube.attributes["model_name"] = model_name
547 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
550def scores_additive_bias(
551 cubes: CubeList,
552 preserved_coordinates: list[str] | str | None = None,
553 obs_model_comparison: bool = False,
554):
555 r"""Calculate the Additive Bias (Mean Error) using scores.
557 Acts as a wrapper around the ME calculation from ``scores`` ([scoresa]_, [scoresb]_).
559 Parameters
560 ----------
561 cubes: iris.cube.CubeList
562 A CubeList containing exactly two cubes: a base and an "other" model,
563 this can be an analysis and the model.
564 preserved_coordinates: list[str] | str | None, default is None.
565 The coordinates that you wish to preserve in the calculaiton of the
566 ME. For example if you want a map of each time you can preserve
567 ["time","grid_latitude", "grid_longitude"] or if you want a time series
568 you can preserve ["time"], if you want to collapse to a single value
569 use `None`. The default is `None`.
571 Returns
572 -------
573 scores_cubelist: iris.cube.CubeList
574 A cubelist containing the ME between the base and other cube(s).
575 """
576 scores_cubelist = CubeList()
577 if obs_model_comparison:
578 for cb in cubes:
579 if "observed" in cb.long_name:
580 base = cb
581 else:
582 others = [cb]
583 else:
584 base, others = _sort_cube_into_base_and_other(cubes)
586 for other in others:
587 base, other = _process_cubes_for_verification(base, other)
589 # Copy the coordinates of the input cubes.
590 other_xr = xr.DataArray.from_iris(other)
591 base_xr = xr.DataArray.from_iris(base)
592 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
594 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
595 # apply scores, and then transform it back.
596 scores_cube = xr.DataArray.to_iris(
597 scores.continuous.additive_bias(
598 other_xr,
599 base_xr,
600 preserve_dims=preserve_dims,
601 )
602 )
604 # If time is aggregated out, attach a scalar time coordinate with bounds
605 # so plotting can display the aggregated period in the title.
606 try:
607 if not scores_cube.coords("time"):
608 base_time = base.coord("time")
609 time_vals = (
610 base_time.bounds.flatten()
611 if base_time.has_bounds()
612 else base_time.points
613 )
614 t_start = float(time_vals[0])
615 t_end = float(time_vals[-1])
616 t_mid = 0.5 * (t_start + t_end)
618 scores_cube.add_aux_coord(
619 iris.coords.AuxCoord(
620 t_mid,
621 standard_name=base_time.standard_name,
622 long_name=base_time.long_name,
623 var_name=base_time.var_name,
624 units=base_time.units,
625 bounds=np.array([t_start, t_end]),
626 attributes=base_time.attributes.copy(),
627 )
628 )
629 except iris.exceptions.CoordinateNotFoundError:
630 pass
631 scores_cube.rename(f"Additive_Bias_of_{base.name()}")
632 scores_cubelist.append(scores_cube)
633 model_name = other.attributes["model_name"]
634 scores_cube.attributes["model_name"] = model_name
636 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
639def scores_correlation_pearsonr(
640 cubes: CubeList,
641 preserved_coordinates: list[str] | str | None = None,
642 obs_model_comparison: bool = False,
643):
644 r"""Calculate the Pearson's Correlation (PC) coefficient using scores.
646 Acts as a wrapper around the PC calculation from ``scores`` ([scoresa]_, [scoresb]_).
648 Parameters
649 ----------
650 cubes: iris.cube.CubeList
651 A CubeList containing exactly two cubes: a base and an "other" model,
652 this can be an analysis and the model.
653 preserved_coordinates: list[str] | str | None, default is None.
654 The coordinates that you wish to preserve in the calculation of the
655 PC. For example if you want a map of each time you can preserve
656 ["time","grid_latitude", "grid_longitude"] or if you want a time series
657 you can preserve ["time"], if you want to collapse to a single value
658 use `None`. The default is `None`.
660 Returns
661 -------
662 scores_cubelist: iris.cube.CubeList
663 A cubelist containing the PC between the base and other cube(s).
664 """
665 scores_cubelist = CubeList()
666 if obs_model_comparison:
667 for cb in cubes:
668 if "observed" in cb.long_name:
669 base = cb
670 else:
671 others = [cb]
672 else:
673 base, others = _sort_cube_into_base_and_other(cubes)
675 for other in others:
676 base, other = _process_cubes_for_verification(base, other)
678 # Copy the coordinates of the input cubes.
679 other_xr = xr.DataArray.from_iris(other)
680 base_xr = xr.DataArray.from_iris(base)
681 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
683 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
684 # apply scores, and then transform it back.
685 scores_cube = xr.DataArray.to_iris(
686 scores.continuous.correlation.pearsonr(
687 other_xr,
688 base_xr,
689 preserve_dims=preserve_dims,
690 )
691 )
693 # If time is aggregated out, attach a scalar time coordinate with bounds
694 # so plotting can display the aggregated period in the title.
695 try:
696 if not scores_cube.coords("time"):
697 base_time = base.coord("time")
698 time_vals = (
699 base_time.bounds.flatten()
700 if base_time.has_bounds()
701 else base_time.points
702 )
703 t_start = float(time_vals[0])
704 t_end = float(time_vals[-1])
705 t_mid = 0.5 * (t_start + t_end)
707 scores_cube.add_aux_coord(
708 iris.coords.AuxCoord(
709 t_mid,
710 standard_name=base_time.standard_name,
711 long_name=base_time.long_name,
712 var_name=base_time.var_name,
713 units=base_time.units,
714 bounds=np.array([t_start, t_end]),
715 attributes=base_time.attributes.copy(),
716 )
717 )
718 except iris.exceptions.CoordinateNotFoundError:
719 pass
721 scores_cube.rename(f"Pearson_Correlation_of_{base.name()}")
722 scores_cubelist.append(scores_cube)
723 model_name = other.attributes["model_name"]
724 scores_cube.attributes["model_name"] = model_name
725 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
728def scores_crps_for_ensemble(
729 cubes: Cube | CubeList, method: str = "ecdf", control_member: int = 0
730) -> iris.Constraint:
731 r"""Calculate the CRPS for an ensemble.
733 Acts as a wrapper around the crps_for_ensemble from ``scores`` ([scoresa]_, [scoresb]_).
735 Lower CRPS values are better (implies experiment distribution is closer to control distribution/observations),
736 larger values are worse (implies distributions are dissimilar).
737 It is applicable across time and spatial scales as the focus is on the distribution of the values.
738 Default method is ecdf. ecdf is exact value from the empirical distributions,
739 whereas fair produces an approximated value based on a random sample of the underlying distribution.
741 See [CRPS]_ for further information.
743 Parameters
744 ----------
745 cubes: iris.cube.Cube
746 A Cube containing ensembles data
748 Returns
749 -------
750 crps: iris.cube.Cube
751 A cube containing the crps between the ensemble members and the control
752 """
753 if control_member != 0:
754 logger.warning("control member is usual 0")
756 if control_member not in cubes.coords("realization")[0].points:
757 new_control_member = cubes.coords("realization")[0].points[0]
758 logger.warning(
759 f"control member value {control_member} out of bounds, defaulting to control member={new_control_member}"
760 )
761 control_member = new_control_member
763 if cubes.coord("time").shape[0] == 1:
764 raise ValueError("Cube has only one time point.")
766 if cubes.coord("realization").shape[0] < 3:
767 raise ValueError("Cube should have one control member and at least two members")
769 ctrl = cubes.extract(generate_realization_constraint([control_member]))
770 ens_mem = cubes.extract(
771 generate_remove_single_ensemble_member_constraint(control_member)
772 )
774 # Realising the data in advance provides a large speedup
775 _ = ctrl.data
776 _ = ens_mem.data
777 del _
779 ctrl = xr.DataArray.from_iris(ctrl)
780 ens_mem = xr.DataArray.from_iris(ens_mem)
782 crps = xr.DataArray.to_iris(
783 scores.probability.crps_for_ensemble(
784 ens_mem,
785 ctrl,
786 ensemble_member_dim="realization",
787 method=method,
788 preserve_dims="time",
789 )
790 )
792 crps.rename(f"CRPS_of_{cubes[0].name()}")
793 _realization_callback(crps)
794 return crps
797def scores_pod_model_obs(
798 cubes: CubeList,
799 preserved_coordinates: list[str] | str | None,
800 threshold: str,
801 op_func: str,
802):
803 r"""
804 Compute the Probability of Detection (POD) score using Scores ([scoresa]_ [scoresb]_).
806 Parameters
807 ----------
808 cubes: iris.cube.CubeList
809 An iris cubelist containing model(s) and an observation cube.
810 preserved_coordinates: list | str | None
811 An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
812 then preserving coordinate 'time' will produce a probability of detection score for each timeslice (shape time). If None,
813 then it will return a single value score for all times/point locations.
814 threshold: str
815 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).
816 op_func: str
817 A string either containing 'lt' for less than or 'gt for greater than, to determine how the threshold is applied to the data
818 to generate the mask.
820 Returns
821 -------
822 cube: iris.cube
823 An iris cube, containing the probability of detection score for further plotting.
825 Notes
826 -----
827 The probability of detection calculates the proportion of observed events that meet a threshold that were correctly forecast by the model.
828 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
829 295k, that would be a positive hit. It does not take into account how far above/below a threshold a model forecasts.
831 It is calculated as .. math:: POD = \frac{true positives}{true positives + false negatives}
833 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.
835 POD produces a range of 0 to 1, where 1 is a perfect score.
836 """
837 # Split out model(s) and obs
838 models = CubeList()
839 for c in cubes:
840 if "observed" in c.long_name:
841 observed = c
842 else:
843 models.append(c)
845 # Setup cubelist to store results
846 scores_results = iris.cube.CubeList()
848 # Setup operators greater than, less than.
849 ops = {
850 "gt": operator.gt,
851 "lt": operator.lt,
852 }
854 try:
855 op = ops[op_func]
856 except KeyError as err:
857 raise ValueError(f"Operator {op_func} not supported.") from err
859 for model in models:
860 # Convert obs cubes to xarray and resolve preserved dimensions.
861 other_xr = xr.DataArray.from_iris(model)
862 base_xr = xr.DataArray.from_iris(observed)
863 preserve_dims = _resolve_preserve_dims(
864 observed, other_xr, preserved_coordinates
865 )
867 # Create event operator object using threshold and operator direction.
868 event_operator = scores.categorical.ThresholdEventOperator(
869 default_event_threshold=float(threshold), default_op_fn=op
870 )
872 # Generate binary fields using the event operator.
873 forecast_binary, observed_binary = event_operator.make_event_tables(
874 other_xr, base_xr
875 )
877 # Create binary contigency manager, as per Scores API, using transform to preserve preserve_dims
878 contingency_manager = scores.categorical.BinaryContingencyManager(
879 forecast_binary, observed_binary
880 ).transform(preserve_dims=preserve_dims)
882 # Get POD from the contigency manager, and convert back to an iris cube.
883 scores_cube = xr.DataArray.to_iris(
884 contingency_manager.probability_of_detection()
885 )
887 # Rename cube so it plots correctly alongside correcting cube units.
888 scores_cube.rename(
889 f"Probability_Of_Detection_{op_func}_{threshold}_{observed.name()}"
890 )
891 scores_cube.units = "1"
892 scores_cube.attributes["model_name"] = model.attributes["model_name"]
894 scores_results.append(scores_cube)
896 return scores_results
899def scores_ets_model_obs(
900 cubes: CubeList,
901 preserved_coordinates: list[str] | str | None,
902 threshold: str,
903 op_func: str,
904):
905 r"""
906 Compute the Equitable Threat Score (ETS) score using Scores ([scoresa]_ [scoresb]_).
908 Parameters
909 ----------
910 cubes: iris.cube.CubeList
911 An iris cubelist containing model(s) and an observation cube.
912 preserved_coordinates: list | str | None
913 An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
914 then preserving coordinate 'time' will produce the equitable threat score for each timeslice (shape time). If None,
915 then it will return a single value score for all times/point locations.
916 threshold: str
917 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).
918 op_func: str
919 A string either containing 'lt' for less than or 'gt for greater than, to determine how the threshold is applied to the data
920 to generate the mask.
922 Returns
923 -------
924 cube: iris.cube
925 An iris cube, containing the probability of detection score for further plotting.
927 Notes
928 -----
929 The Equitable Threat Score (ETS) evaluates the accuracy of forecasts for events that meet a specified threshold,
930 hile accounting for correct forecasts that could occur purely by chance. Unlike the Probability of Detection (POD),
931 ETS considers hits, misses, and false alarms, providing a more balanced assessment of forecast skill.
933 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
934 would be counted as a hit. ETS adjusts the total number of hits by removing the number of hits expected due to random chance.
936 It is calculated as:
938 .. math::
940 ETS = \frac{hits - hits_{random}}
941 {hits + misses + false\ alarms - hits_{random}}
943 where
945 hits_{random} = \frac{(hits + misses)(hits + false\ alarms)}{total count}
947 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
948 random chance.
949 """
950 # Split out model(s) and obs
951 models = CubeList()
952 for c in cubes:
953 if "observed" in c.long_name:
954 observed = c
955 else:
956 models.append(c)
958 # Setup cubelist to store results
959 scores_results = iris.cube.CubeList()
961 # Setup operators greater than, less than.
962 ops = {
963 "gt": operator.gt,
964 "lt": operator.lt,
965 }
967 try:
968 op = ops[op_func]
969 except KeyError as err:
970 raise ValueError(f"Operator {op_func} not supported.") from err
972 for model in models:
973 # Convert obs cubes to xarray and resolve preserved dimensions.
974 other_xr = xr.DataArray.from_iris(model)
975 base_xr = xr.DataArray.from_iris(observed)
976 preserve_dims = _resolve_preserve_dims(
977 observed, other_xr, preserved_coordinates
978 )
980 # Create event operator object using threshold and operator direction.
981 event_operator = scores.categorical.ThresholdEventOperator(
982 default_event_threshold=float(threshold), default_op_fn=op
983 )
985 # Generate binary fields using the event operator.
986 forecast_binary, observed_binary = event_operator.make_event_tables(
987 other_xr, base_xr
988 )
990 # Create binary contigency manager, as per Scores API, using transform to preserve preserve_dims
991 contingency_manager = scores.categorical.BinaryContingencyManager(
992 forecast_binary, observed_binary
993 ).transform(preserve_dims=preserve_dims)
995 # Get ETS from the contigency manager, and convert back to an iris cube.
996 scores_cube = xr.DataArray.to_iris(contingency_manager.equitable_threat_score())
998 # Rename cube so it plots correctly alongside correcting cube units.
999 scores_cube.rename(
1000 f"Equitable_Threat_Score_{op_func}_{threshold}_{observed.name()}"
1001 )
1002 scores_cube.units = "1"
1003 scores_cube.attributes["model_name"] = model.attributes["model_name"]
1005 scores_results.append(scores_cube)
1007 return scores_results