Coverage for src/CSET/operators/scoreswrappers.py: 91%
281 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-21 07:38 +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_correlation_pearsonr_model_obs(
284 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
285):
286 r"""Calculate the Pearson's Correlation (PC) coefficient using scores.
288 Acts as a wrapper around the PC calculation from ``scores`` ([scoresa]_, [scoresb]_).
290 Parameters
291 ----------
292 cubes: iris.cube.CubeList
293 A CubeList containing exactly two cubes: a base and an "other" model,
294 this can be an analysis and the model.
295 preserved_coordinates: list[str] | str | None, default is None.
296 The coordinates that you wish to preserve in the calculation of the
297 PC. For example if you want a map of each time you can preserve
298 ["time","latitude", "longitude"] or if you want a time series
299 you can preserve ["time"], if you want to collapse to a single value
300 use `None`. The default is `None`.
302 Returns
303 -------
304 scores_cube: iris.cube.CubeList
305 A cube list containing the PC between the models and observation cube.
306 """
307 pearsonr_cubes = CubeList()
308 model_list = CubeList()
310 for cb in cubes:
311 if "observed" in cb.long_name:
312 observed = cb
313 else:
314 model_list.append(cb)
316 for model in model_list:
317 input_cubelist = CubeList()
318 input_cubelist.append(observed)
319 input_cubelist.append(model)
320 pearsonr = scores_correlation_pearsonr(
321 input_cubelist, preserved_coordinates, obs_model_comparison=True
322 )
323 model_name = model.attributes["model_name"]
324 pearsonr.attributes["model_name"] = model_name
325 pearsonr_cubes.append(pearsonr)
327 return pearsonr_cubes
330def scores_rmse(
331 cubes: CubeList,
332 preserved_coordinates: list[str] | str | None = None,
333 obs_model_comparison: bool = False,
334):
335 r"""Calculate the Root Mean Square Error (RMSE) using scores.
337 Acts as a wrapper around the RMSE calculation from ``scores`` ([scoresa]_, [scoresb]_).
338 It is calculated as
340 .. math:: RMSE = \sqrt{\frac{1}{N} \Sigma(forecast - observations)^2}
342 Parameters
343 ----------
344 cubes: iris.cube.CubeList
345 A CubeList containing exactly two cubes: a base and an "other" model,
346 this can be an analysis and the model.
347 preserved_coordinates: list[str] | str | None, default is None.
348 The coordinates (or xarray dimension names) that you wish to preserve in the calculaiton of the
349 RMSE. For example if you want a map of each time you can preserve
350 ["time","grid_latitude", "grid_longitude"] or if you want a time series
351 you can preserve ["time"], if you want to collapse to a single value
352 use `None`. The default is `None`.
353 obs_model_comparison: bool, default False
354 Set true if doing model-obs comparison.
356 Returns
357 -------
358 scores_cubelist: iris.cube.CubeList
359 A cubelist containing the RMSE between the base and other cube.
360 """
361 scores_cubelist = CubeList()
362 if obs_model_comparison:
363 for cb in cubes:
364 if "observed" in cb.long_name:
365 base = cb
366 else:
367 others = [cb]
368 else:
369 base, others = _sort_cube_into_base_and_other(cubes)
371 for other in others:
372 base, other = _process_cubes_for_verification(base, other)
374 # Copy the coordinates of the input cubes.
375 other_xr = xr.DataArray.from_iris(other)
376 base_xr = xr.DataArray.from_iris(base)
377 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
379 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
380 # apply scores, and then transform it back.
381 scores_cube = xr.DataArray.to_iris(
382 scores.continuous.rmse(
383 other_xr,
384 base_xr,
385 preserve_dims=preserve_dims,
386 )
387 )
389 # If time is aggregated out, attach a scalar time coordinate with bounds
390 # so plotting can display the aggregated period in the title.
391 try:
392 if not scores_cube.coords("time"):
393 base_time = base.coord("time")
394 time_vals = (
395 base_time.bounds.flatten()
396 if base_time.has_bounds()
397 else base_time.points
398 )
399 t_start = float(time_vals[0])
400 t_end = float(time_vals[-1])
401 t_mid = 0.5 * (t_start + t_end)
403 scores_cube.add_aux_coord(
404 iris.coords.AuxCoord(
405 t_mid,
406 standard_name=base_time.standard_name,
407 long_name=base_time.long_name,
408 var_name=base_time.var_name,
409 units=base_time.units,
410 bounds=np.array([t_start, t_end]),
411 attributes=base_time.attributes.copy(),
412 )
413 )
414 except iris.exceptions.CoordinateNotFoundError:
415 pass
417 scores_cube.rename(f"RMSE_of_{base.name()}")
418 scores_cubelist.append(scores_cube)
420 model_name = other.attributes["model_name"]
421 scores_cube.attributes["model_name"] = model_name
423 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
426def scores_mae(cubes: CubeList, preserved_coordinates: list[str] | str | None = None):
427 r"""Calculate the Mean Absolute Error (MAE) using scores.
429 Acts as a wrapper around the MAE calculation from ``scores`` ([scoresa]_, [scoresb]_).
431 Parameters
432 ----------
433 cubes: iris.cube.CubeList
434 A CubeList containing exactly two cubes: a base and an "other" model,
435 this can be an analysis and the model.
436 preserved_coordinates: list[str] | str | None, default is None.
437 The coordinates that you wish to preserve in the calculaiton of the
438 MAE. For example if you want a map of each time you can preserve
439 ["time","grid_latitude", "grid_longitude"] or if you want a time series
440 you can preserve ["time"], if you want to collapse to a single value
441 use `None`. The default is `None`.
443 Returns
444 -------
445 scores_cubelist: iris.cube.CubeList
446 A cubelist containing the MAE between the base and other cube(s).
447 """
448 base, others = _sort_cube_into_base_and_other(cubes)
449 scores_cubelist = CubeList()
450 for other in others:
451 base, other = _process_cubes_for_verification(base, other)
453 # Copy the coordinates of the input cubes.
454 other_xr = xr.DataArray.from_iris(other)
455 base_xr = xr.DataArray.from_iris(base)
456 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
458 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
459 # apply scores, and then transform it back.
460 scores_cube = xr.DataArray.to_iris(
461 scores.continuous.mae(
462 other_xr,
463 base_xr,
464 preserve_dims=preserve_dims,
465 )
466 )
468 # If time is aggregated out, attach a scalar time coordinate with bounds
469 # so plotting can display the aggregated period in the title.
470 try:
471 if not scores_cube.coords("time"): 471 ↛ 496line 471 didn't jump to line 496 because the condition on line 471 was always true
472 base_time = base.coord("time")
473 time_vals = (
474 base_time.bounds.flatten()
475 if base_time.has_bounds()
476 else base_time.points
477 )
478 t_start = float(time_vals[0])
479 t_end = float(time_vals[-1])
480 t_mid = 0.5 * (t_start + t_end)
482 scores_cube.add_aux_coord(
483 iris.coords.AuxCoord(
484 t_mid,
485 standard_name=base_time.standard_name,
486 long_name=base_time.long_name,
487 var_name=base_time.var_name,
488 units=base_time.units,
489 bounds=np.array([t_start, t_end]),
490 attributes=base_time.attributes.copy(),
491 )
492 )
493 except iris.exceptions.CoordinateNotFoundError:
494 pass
496 scores_cube.rename(f"MAE_of_{base.name()}")
497 scores_cubelist.append(scores_cube)
498 model_name = other.attributes["model_name"]
499 scores_cube.attributes["model_name"] = model_name
501 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
504def scores_additive_bias(
505 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
506):
507 r"""Calculate the Additive Bias (Mean Error) using scores.
509 Acts as a wrapper around the ME calculation from ``scores`` ([scoresa]_, [scoresb]_).
511 Parameters
512 ----------
513 cubes: iris.cube.CubeList
514 A CubeList containing exactly two cubes: a base and an "other" model,
515 this can be an analysis and the model.
516 preserved_coordinates: list[str] | str | None, default is None.
517 The coordinates that you wish to preserve in the calculaiton of the
518 ME. For example if you want a map of each time you can preserve
519 ["time","grid_latitude", "grid_longitude"] or if you want a time series
520 you can preserve ["time"], if you want to collapse to a single value
521 use `None`. The default is `None`.
523 Returns
524 -------
525 scores_cubelist: iris.cube.CubeList
526 A cubelist containing the ME between the base and other cube(s).
527 """
528 base, others = _sort_cube_into_base_and_other(cubes)
529 scores_cubelist = CubeList()
530 for other in others:
531 base, other = _process_cubes_for_verification(base, other)
533 # Copy the coordinates of the input cubes.
534 other_xr = xr.DataArray.from_iris(other)
535 base_xr = xr.DataArray.from_iris(base)
536 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
538 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
539 # apply scores, and then transform it back.
540 scores_cube = xr.DataArray.to_iris(
541 scores.continuous.additive_bias(
542 other_xr,
543 base_xr,
544 preserve_dims=preserve_dims,
545 )
546 )
548 # If time is aggregated out, attach a scalar time coordinate with bounds
549 # so plotting can display the aggregated period in the title.
550 try:
551 if not scores_cube.coords("time"): 551 ↛ 575line 551 didn't jump to line 575 because the condition on line 551 was always true
552 base_time = base.coord("time")
553 time_vals = (
554 base_time.bounds.flatten()
555 if base_time.has_bounds()
556 else base_time.points
557 )
558 t_start = float(time_vals[0])
559 t_end = float(time_vals[-1])
560 t_mid = 0.5 * (t_start + t_end)
562 scores_cube.add_aux_coord(
563 iris.coords.AuxCoord(
564 t_mid,
565 standard_name=base_time.standard_name,
566 long_name=base_time.long_name,
567 var_name=base_time.var_name,
568 units=base_time.units,
569 bounds=np.array([t_start, t_end]),
570 attributes=base_time.attributes.copy(),
571 )
572 )
573 except iris.exceptions.CoordinateNotFoundError:
574 pass
575 scores_cube.rename(f"Additive_Bias_of_{base.name()}")
576 scores_cubelist.append(scores_cube)
577 model_name = other.attributes["model_name"]
578 scores_cube.attributes["model_name"] = model_name
580 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
583def scores_correlation_pearsonr(
584 cubes: CubeList,
585 preserved_coordinates: list[str] | str | None = None,
586 obs_model_comparison: bool = False,
587):
588 r"""Calculate the Pearson's Correlation (PC) coefficient using scores.
590 Acts as a wrapper around the PC calculation from ``scores`` ([scoresa]_, [scoresb]_).
592 Parameters
593 ----------
594 cubes: iris.cube.CubeList
595 A CubeList containing exactly two cubes: a base and an "other" model,
596 this can be an analysis and the model.
597 preserved_coordinates: list[str] | str | None, default is None.
598 The coordinates that you wish to preserve in the calculation of the
599 PC. For example if you want a map of each time you can preserve
600 ["time","grid_latitude", "grid_longitude"] or if you want a time series
601 you can preserve ["time"], if you want to collapse to a single value
602 use `None`. The default is `None`.
604 Returns
605 -------
606 scores_cubelist: iris.cube.CubeList
607 A cubelist containing the PC between the base and other cube(s).
608 """
609 scores_cubelist = CubeList()
610 if obs_model_comparison:
611 for cb in cubes:
612 if "observed" in cb.long_name:
613 base = cb
614 else:
615 others = [cb]
616 else:
617 base, others = _sort_cube_into_base_and_other(cubes)
619 for other in others:
620 base, other = _process_cubes_for_verification(base, other)
622 # Copy the coordinates of the input cubes.
623 other_xr = xr.DataArray.from_iris(other)
624 base_xr = xr.DataArray.from_iris(base)
625 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
627 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
628 # apply scores, and then transform it back.
629 scores_cube = xr.DataArray.to_iris(
630 scores.continuous.correlation.pearsonr(
631 other_xr,
632 base_xr,
633 preserve_dims=preserve_dims,
634 )
635 )
637 # If time is aggregated out, attach a scalar time coordinate with bounds
638 # so plotting can display the aggregated period in the title.
639 try:
640 if not scores_cube.coords("time"):
641 base_time = base.coord("time")
642 time_vals = (
643 base_time.bounds.flatten()
644 if base_time.has_bounds()
645 else base_time.points
646 )
647 t_start = float(time_vals[0])
648 t_end = float(time_vals[-1])
649 t_mid = 0.5 * (t_start + t_end)
651 scores_cube.add_aux_coord(
652 iris.coords.AuxCoord(
653 t_mid,
654 standard_name=base_time.standard_name,
655 long_name=base_time.long_name,
656 var_name=base_time.var_name,
657 units=base_time.units,
658 bounds=np.array([t_start, t_end]),
659 attributes=base_time.attributes.copy(),
660 )
661 )
662 except iris.exceptions.CoordinateNotFoundError:
663 pass
665 scores_cube.rename(f"Pearson_Correlation_of_{base.name()}")
666 scores_cubelist.append(scores_cube)
667 model_name = other.attributes["model_name"]
668 scores_cube.attributes["model_name"] = model_name
669 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
672def scores_crps_for_ensemble(
673 cubes: Cube | CubeList, method: str = "ecdf", control_member: int = 0
674) -> iris.Constraint:
675 r"""Calculate the CRPS for an ensemble.
677 Acts as a wrapper around the crps_for_ensemble from ``scores`` ([scoresa]_, [scoresb]_).
679 Lower CRPS values are better (implies experiment distribution is closer to control distribution/observations),
680 larger values are worse (implies distributions are dissimilar).
681 It is applicable across time and spatial scales as the focus is on the distribution of the values.
682 Default method is ecdf. ecdf is exact value from the empirical distributions,
683 whereas fair produces an approximated value based on a random sample of the underlying distribution.
685 See [CRPS]_ for further information.
687 Parameters
688 ----------
689 cubes: iris.cube.Cube
690 A Cube containing ensembles data
692 Returns
693 -------
694 crps: iris.cube.Cube
695 A cube containing the crps between the ensemble members and the control
696 """
697 if control_member != 0:
698 logger.warning("control member is usual 0")
700 if control_member not in cubes.coords("realization")[0].points:
701 new_control_member = cubes.coords("realization")[0].points[0]
702 logger.warning(
703 f"control member value {control_member} out of bounds, defaulting to control member={new_control_member}"
704 )
705 control_member = new_control_member
707 if cubes.coord("time").shape[0] == 1:
708 raise ValueError("Cube has only one time point.")
710 if cubes.coord("realization").shape[0] < 3:
711 raise ValueError("Cube should have one control member and at least two members")
713 ctrl = cubes.extract(generate_realization_constraint([control_member]))
714 ens_mem = cubes.extract(
715 generate_remove_single_ensemble_member_constraint(control_member)
716 )
718 # Realising the data in advance provides a large speedup
719 _ = ctrl.data
720 _ = ens_mem.data
721 del _
723 ctrl = xr.DataArray.from_iris(ctrl)
724 ens_mem = xr.DataArray.from_iris(ens_mem)
726 crps = xr.DataArray.to_iris(
727 scores.probability.crps_for_ensemble(
728 ens_mem,
729 ctrl,
730 ensemble_member_dim="realization",
731 method=method,
732 preserve_dims="time",
733 )
734 )
736 crps.rename(f"CRPS_of_{cubes[0].name()}")
737 _realization_callback(crps)
738 return crps
741def scores_pod_model_obs(
742 cubes: CubeList,
743 preserved_coordinates: list[str] | str | None,
744 threshold: str,
745 op_func: str,
746):
747 r"""
748 Compute the Probability of Detection (POD) score using Scores ([scoresa]_ [scoresb]_).
750 Parameters
751 ----------
752 cubes: iris.cube.CubeList
753 An iris cubelist containing model(s) and an observation cube.
754 preserved_coordinates: list | str | None
755 An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
756 then preserving coordinate 'time' will produce a probability of detection score for each timeslice (shape time). If None,
757 then it will return a single value score for all times/point locations.
758 threshold: str
759 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).
760 op_func: str
761 A string either containing 'lt' for less than or 'gt for greater than, to determine how the threshold is applied to the data
762 to generate the mask.
764 Returns
765 -------
766 cube: iris.cube
767 An iris cube, containing the probability of detection score for further plotting.
769 Notes
770 -----
771 The probability of detection calculates the proportion of observed events that meet a threshold that were correctly forecast by the model.
772 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
773 295k, that would be a positive hit. It does not take into account how far above/below a threshold a model forecasts.
775 It is calculated as .. math:: POD = \frac{true positives}{true positives + false negatives}
777 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.
779 POD produces a range of 0 to 1, where 1 is a perfect score.
780 """
781 # Split out model(s) and obs
782 models = CubeList()
783 for c in cubes:
784 if "observed" in c.long_name:
785 observed = c
786 else:
787 models.append(c)
789 # Setup cubelist to store results
790 scores_results = iris.cube.CubeList()
792 # Setup operators greater than, less than.
793 ops = {
794 "gt": operator.gt,
795 "lt": operator.lt,
796 }
798 try:
799 op = ops[op_func]
800 except KeyError as err:
801 raise ValueError(f"Operator {op_func} not supported.") from err
803 for model in models:
804 # Convert obs cubes to xarray and resolve preserved dimensions.
805 other_xr = xr.DataArray.from_iris(model)
806 base_xr = xr.DataArray.from_iris(observed)
807 preserve_dims = _resolve_preserve_dims(
808 observed, other_xr, preserved_coordinates
809 )
811 # Create event operator object using threshold and operator direction.
812 event_operator = scores.categorical.ThresholdEventOperator(
813 default_event_threshold=float(threshold), default_op_fn=op
814 )
816 # Generate binary fields using the event operator.
817 forecast_binary, observed_binary = event_operator.make_event_tables(
818 other_xr, base_xr
819 )
821 # Create binary contigency manager, as per Scores API, using transform to preserve preserve_dims
822 contingency_manager = scores.categorical.BinaryContingencyManager(
823 forecast_binary, observed_binary
824 ).transform(preserve_dims=preserve_dims)
826 # Get POD from the contigency manager, and convert back to an iris cube.
827 scores_cube = xr.DataArray.to_iris(
828 contingency_manager.probability_of_detection()
829 )
831 # Rename cube so it plots correctly alongside correcting cube units.
832 scores_cube.rename(
833 f"Probability_Of_Detection_{op_func}_{threshold}_{observed.name()}"
834 )
835 scores_cube.units = "1"
836 scores_cube.attributes["model_name"] = model.attributes["model_name"]
838 scores_results.append(scores_cube)
840 return scores_results
843def scores_ets_model_obs(
844 cubes: CubeList,
845 preserved_coordinates: list[str] | str | None,
846 threshold: str,
847 op_func: str,
848):
849 r"""
850 Compute the Equitable Threat Score (ETS) score using Scores ([scoresa]_ [scoresb]_).
852 Parameters
853 ----------
854 cubes: iris.cube.CubeList
855 An iris cubelist containing model(s) and an observation cube.
856 preserved_coordinates: list | str | None
857 An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
858 then preserving coordinate 'time' will produce the equitable threat score for each timeslice (shape time). If None,
859 then it will return a single value score for all times/point locations.
860 threshold: str
861 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).
862 op_func: str
863 A string either containing 'lt' for less than or 'gt for greater than, to determine how the threshold is applied to the data
864 to generate the mask.
866 Returns
867 -------
868 cube: iris.cube
869 An iris cube, containing the probability of detection score for further plotting.
871 Notes
872 -----
873 The Equitable Threat Score (ETS) evaluates the accuracy of forecasts for events that meet a specified threshold,
874 hile accounting for correct forecasts that could occur purely by chance. Unlike the Probability of Detection (POD),
875 ETS considers hits, misses, and false alarms, providing a more balanced assessment of forecast skill.
877 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
878 would be counted as a hit. ETS adjusts the total number of hits by removing the number of hits expected due to random chance.
880 It is calculated as:
882 .. math::
884 ETS = \frac{hits - hits_{random}}
885 {hits + misses + false\ alarms - hits_{random}}
887 where
889 hits_{random} = \frac{(hits + misses)(hits + false\ alarms)}{total count}
891 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
892 random chance.
893 """
894 # Split out model(s) and obs
895 models = CubeList()
896 for c in cubes:
897 if "observed" in c.long_name:
898 observed = c
899 else:
900 models.append(c)
902 # Setup cubelist to store results
903 scores_results = iris.cube.CubeList()
905 # Setup operators greater than, less than.
906 ops = {
907 "gt": operator.gt,
908 "lt": operator.lt,
909 }
911 try:
912 op = ops[op_func]
913 except KeyError as err:
914 raise ValueError(f"Operator {op_func} not supported.") from err
916 for model in models:
917 # Convert obs cubes to xarray and resolve preserved dimensions.
918 other_xr = xr.DataArray.from_iris(model)
919 base_xr = xr.DataArray.from_iris(observed)
920 preserve_dims = _resolve_preserve_dims(
921 observed, other_xr, preserved_coordinates
922 )
924 # Create event operator object using threshold and operator direction.
925 event_operator = scores.categorical.ThresholdEventOperator(
926 default_event_threshold=float(threshold), default_op_fn=op
927 )
929 # Generate binary fields using the event operator.
930 forecast_binary, observed_binary = event_operator.make_event_tables(
931 other_xr, base_xr
932 )
934 # Create binary contigency manager, as per Scores API, using transform to preserve preserve_dims
935 contingency_manager = scores.categorical.BinaryContingencyManager(
936 forecast_binary, observed_binary
937 ).transform(preserve_dims=preserve_dims)
939 # Get ETS from the contigency manager, and convert back to an iris cube.
940 scores_cube = xr.DataArray.to_iris(contingency_manager.equitable_threat_score())
942 # Rename cube so it plots correctly alongside correcting cube units.
943 scores_cube.rename(
944 f"Equitable_Threat_Score_{op_func}_{threshold}_{observed.name()}"
945 )
946 scores_cube.units = "1"
947 scores_cube.attributes["model_name"] = model.attributes["model_name"]
949 scores_results.append(scores_cube)
951 return scores_results