Coverage for src/CSET/operators/scoreswrappers.py: 90%
260 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-20 15:26 +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_rmse(
284 cubes: CubeList,
285 preserved_coordinates: list[str] | str | None = None,
286 obs_model_comparison: bool = False,
287):
288 r"""Calculate the Root Mean Square Error (RMSE) using scores.
290 Acts as a wrapper around the RMSE calculation from ``scores`` ([scoresa]_, [scoresb]_).
291 It is calculated as
293 .. math:: RMSE = \sqrt{\frac{1}{N} \Sigma(forecast - observations)^2}
295 Parameters
296 ----------
297 cubes: iris.cube.CubeList
298 A CubeList containing exactly two cubes: a base and an "other" model,
299 this can be an analysis and the model.
300 preserved_coordinates: list[str] | str | None, default is None.
301 The coordinates (or xarray dimension names) that you wish to preserve in the calculaiton of the
302 RMSE. For example if you want a map of each time you can preserve
303 ["time","grid_latitude", "grid_longitude"] or if you want a time series
304 you can preserve ["time"], if you want to collapse to a single value
305 use `None`. The default is `None`.
306 obs_model_comparison: bool, default False
307 Set true if doing model-obs comparison.
309 Returns
310 -------
311 scores_cubelist: iris.cube.CubeList
312 A cubelist containing the RMSE between the base and other cube.
313 """
314 scores_cubelist = CubeList()
315 if obs_model_comparison:
316 for cb in cubes:
317 if "observed" in cb.long_name:
318 base = cb
319 else:
320 others = [cb]
321 else:
322 base, others = _sort_cube_into_base_and_other(cubes)
324 for other in others:
325 base, other = _process_cubes_for_verification(base, other)
327 # Copy the coordinates of the input cubes.
328 other_xr = xr.DataArray.from_iris(other)
329 base_xr = xr.DataArray.from_iris(base)
330 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
332 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
333 # apply scores, and then transform it back.
334 scores_cube = xr.DataArray.to_iris(
335 scores.continuous.rmse(
336 other_xr,
337 base_xr,
338 preserve_dims=preserve_dims,
339 )
340 )
342 # If time is aggregated out, attach a scalar time coordinate with bounds
343 # so plotting can display the aggregated period in the title.
344 try:
345 if not scores_cube.coords("time"):
346 base_time = base.coord("time")
347 time_vals = (
348 base_time.bounds.flatten()
349 if base_time.has_bounds()
350 else base_time.points
351 )
352 t_start = float(time_vals[0])
353 t_end = float(time_vals[-1])
354 t_mid = 0.5 * (t_start + t_end)
356 scores_cube.add_aux_coord(
357 iris.coords.AuxCoord(
358 t_mid,
359 standard_name=base_time.standard_name,
360 long_name=base_time.long_name,
361 var_name=base_time.var_name,
362 units=base_time.units,
363 bounds=np.array([t_start, t_end]),
364 attributes=base_time.attributes.copy(),
365 )
366 )
367 except iris.exceptions.CoordinateNotFoundError:
368 pass
370 scores_cube.rename(f"RMSE_of_{base.name()}")
371 scores_cubelist.append(scores_cube)
373 model_name = other.attributes["model_name"]
374 scores_cube.attributes["model_name"] = model_name
376 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
379def scores_mae(cubes: CubeList, preserved_coordinates: list[str] | str | None = None):
380 r"""Calculate the Mean Absolute Error (MAE) using scores.
382 Acts as a wrapper around the MAE calculation from ``scores`` ([scoresa]_, [scoresb]_).
384 Parameters
385 ----------
386 cubes: iris.cube.CubeList
387 A CubeList containing exactly two cubes: a base and an "other" model,
388 this can be an analysis and the model.
389 preserved_coordinates: list[str] | str | None, default is None.
390 The coordinates that you wish to preserve in the calculaiton of the
391 MAE. For example if you want a map of each time you can preserve
392 ["time","grid_latitude", "grid_longitude"] or if you want a time series
393 you can preserve ["time"], if you want to collapse to a single value
394 use `None`. The default is `None`.
396 Returns
397 -------
398 scores_cubelist: iris.cube.CubeList
399 A cubelist containing the MAE between the base and other cube(s).
400 """
401 base, others = _sort_cube_into_base_and_other(cubes)
402 scores_cubelist = CubeList()
403 for other in others:
404 base, other = _process_cubes_for_verification(base, other)
406 # Copy the coordinates of the input cubes.
407 other_xr = xr.DataArray.from_iris(other)
408 base_xr = xr.DataArray.from_iris(base)
409 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
411 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
412 # apply scores, and then transform it back.
413 scores_cube = xr.DataArray.to_iris(
414 scores.continuous.mae(
415 other_xr,
416 base_xr,
417 preserve_dims=preserve_dims,
418 )
419 )
421 # If time is aggregated out, attach a scalar time coordinate with bounds
422 # so plotting can display the aggregated period in the title.
423 try:
424 if not scores_cube.coords("time"): 424 ↛ 449line 424 didn't jump to line 449 because the condition on line 424 was always true
425 base_time = base.coord("time")
426 time_vals = (
427 base_time.bounds.flatten()
428 if base_time.has_bounds()
429 else base_time.points
430 )
431 t_start = float(time_vals[0])
432 t_end = float(time_vals[-1])
433 t_mid = 0.5 * (t_start + t_end)
435 scores_cube.add_aux_coord(
436 iris.coords.AuxCoord(
437 t_mid,
438 standard_name=base_time.standard_name,
439 long_name=base_time.long_name,
440 var_name=base_time.var_name,
441 units=base_time.units,
442 bounds=np.array([t_start, t_end]),
443 attributes=base_time.attributes.copy(),
444 )
445 )
446 except iris.exceptions.CoordinateNotFoundError:
447 pass
449 scores_cube.rename(f"MAE_of_{base.name()}")
450 scores_cubelist.append(scores_cube)
451 model_name = other.attributes["model_name"]
452 scores_cube.attributes["model_name"] = model_name
454 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
457def scores_additive_bias(
458 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
459):
460 r"""Calculate the Additive Bias (Mean Error) using scores.
462 Acts as a wrapper around the ME calculation from ``scores`` ([scoresa]_, [scoresb]_).
464 Parameters
465 ----------
466 cubes: iris.cube.CubeList
467 A CubeList containing exactly two cubes: a base and an "other" model,
468 this can be an analysis and the model.
469 preserved_coordinates: list[str] | str | None, default is None.
470 The coordinates that you wish to preserve in the calculaiton of the
471 ME. For example if you want a map of each time you can preserve
472 ["time","grid_latitude", "grid_longitude"] or if you want a time series
473 you can preserve ["time"], if you want to collapse to a single value
474 use `None`. The default is `None`.
476 Returns
477 -------
478 scores_cubelist: iris.cube.CubeList
479 A cubelist containing the ME between the base and other cube(s).
480 """
481 base, others = _sort_cube_into_base_and_other(cubes)
482 scores_cubelist = CubeList()
483 for other in others:
484 base, other = _process_cubes_for_verification(base, other)
486 # Copy the coordinates of the input cubes.
487 other_xr = xr.DataArray.from_iris(other)
488 base_xr = xr.DataArray.from_iris(base)
489 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
491 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
492 # apply scores, and then transform it back.
493 scores_cube = xr.DataArray.to_iris(
494 scores.continuous.additive_bias(
495 other_xr,
496 base_xr,
497 preserve_dims=preserve_dims,
498 )
499 )
501 # If time is aggregated out, attach a scalar time coordinate with bounds
502 # so plotting can display the aggregated period in the title.
503 try:
504 if not scores_cube.coords("time"): 504 ↛ 528line 504 didn't jump to line 528 because the condition on line 504 was always true
505 base_time = base.coord("time")
506 time_vals = (
507 base_time.bounds.flatten()
508 if base_time.has_bounds()
509 else base_time.points
510 )
511 t_start = float(time_vals[0])
512 t_end = float(time_vals[-1])
513 t_mid = 0.5 * (t_start + t_end)
515 scores_cube.add_aux_coord(
516 iris.coords.AuxCoord(
517 t_mid,
518 standard_name=base_time.standard_name,
519 long_name=base_time.long_name,
520 var_name=base_time.var_name,
521 units=base_time.units,
522 bounds=np.array([t_start, t_end]),
523 attributes=base_time.attributes.copy(),
524 )
525 )
526 except iris.exceptions.CoordinateNotFoundError:
527 pass
528 scores_cube.rename(f"Additive_Bias_of_{base.name()}")
529 scores_cubelist.append(scores_cube)
530 model_name = other.attributes["model_name"]
531 scores_cube.attributes["model_name"] = model_name
533 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
536def scores_correlation_pearsonr(
537 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
538):
539 r"""Calculate the Pearson's Correlation (PC) coefficient using scores.
541 Acts as a wrapper around the PC calculation from ``scores`` ([scoresa]_, [scoresb]_).
543 Parameters
544 ----------
545 cubes: iris.cube.CubeList
546 A CubeList containing exactly two cubes: a base and an "other" model,
547 this can be an analysis and the model.
548 preserved_coordinates: list[str] | str | None, default is None.
549 The coordinates that you wish to preserve in the calculation of the
550 PC. For example if you want a map of each time you can preserve
551 ["time","grid_latitude", "grid_longitude"] or if you want a time series
552 you can preserve ["time"], if you want to collapse to a single value
553 use `None`. The default is `None`.
555 Returns
556 -------
557 scores_cubelist: iris.cube.CubeList
558 A cubelist containing the PC between the base and other cube(s).
559 """
560 base, others = _sort_cube_into_base_and_other(cubes)
561 scores_cubelist = CubeList()
562 for other in others:
563 base, other = _process_cubes_for_verification(base, other)
565 # Copy the coordinates of the input cubes.
566 other_xr = xr.DataArray.from_iris(other)
567 base_xr = xr.DataArray.from_iris(base)
568 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
570 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
571 # apply scores, and then transform it back.
572 scores_cube = xr.DataArray.to_iris(
573 scores.continuous.correlation.pearsonr(
574 other_xr,
575 base_xr,
576 preserve_dims=preserve_dims,
577 )
578 )
580 # If time is aggregated out, attach a scalar time coordinate with bounds
581 # so plotting can display the aggregated period in the title.
582 try:
583 if not scores_cube.coords("time"): 583 ↛ 608line 583 didn't jump to line 608 because the condition on line 583 was always true
584 base_time = base.coord("time")
585 time_vals = (
586 base_time.bounds.flatten()
587 if base_time.has_bounds()
588 else base_time.points
589 )
590 t_start = float(time_vals[0])
591 t_end = float(time_vals[-1])
592 t_mid = 0.5 * (t_start + t_end)
594 scores_cube.add_aux_coord(
595 iris.coords.AuxCoord(
596 t_mid,
597 standard_name=base_time.standard_name,
598 long_name=base_time.long_name,
599 var_name=base_time.var_name,
600 units=base_time.units,
601 bounds=np.array([t_start, t_end]),
602 attributes=base_time.attributes.copy(),
603 )
604 )
605 except iris.exceptions.CoordinateNotFoundError:
606 pass
608 scores_cube.rename(f"Pearson_Correlation_of_{base.name()}")
609 scores_cubelist.append(scores_cube)
610 model_name = other.attributes["model_name"]
611 scores_cube.attributes["model_name"] = model_name
612 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
615def scores_crps_for_ensemble(
616 cubes: Cube | CubeList, method: str = "ecdf", control_member: int = 0
617) -> iris.Constraint:
618 r"""Calculate the CRPS for an ensemble.
620 Acts as a wrapper around the crps_for_ensemble from ``scores`` ([scoresa]_, [scoresb]_).
622 Lower CRPS values are better (implies experiment distribution is closer to control distribution/observations),
623 larger values are worse (implies distributions are dissimilar).
624 It is applicable across time and spatial scales as the focus is on the distribution of the values.
625 Default method is ecdf. ecdf is exact value from the empirical distributions,
626 whereas fair produces an approximated value based on a random sample of the underlying distribution.
628 See [CRPS]_ for further information.
630 Parameters
631 ----------
632 cubes: iris.cube.Cube
633 A Cube containing ensembles data
635 Returns
636 -------
637 crps: iris.cube.Cube
638 A cube containing the crps between the ensemble members and the control
639 """
640 if control_member != 0:
641 logger.warning("control member is usual 0")
643 if control_member not in cubes.coords("realization")[0].points:
644 new_control_member = cubes.coords("realization")[0].points[0]
645 logger.warning(
646 f"control member value {control_member} out of bounds, defaulting to control member={new_control_member}"
647 )
648 control_member = new_control_member
650 if cubes.coord("time").shape[0] == 1:
651 raise ValueError("Cube has only one time point.")
653 if cubes.coord("realization").shape[0] < 3:
654 raise ValueError("Cube should have one control member and at least two members")
656 ctrl = cubes.extract(generate_realization_constraint([control_member]))
657 ens_mem = cubes.extract(
658 generate_remove_single_ensemble_member_constraint(control_member)
659 )
661 # Realising the data in advance provides a large speedup
662 _ = ctrl.data
663 _ = ens_mem.data
664 del _
666 ctrl = xr.DataArray.from_iris(ctrl)
667 ens_mem = xr.DataArray.from_iris(ens_mem)
669 crps = xr.DataArray.to_iris(
670 scores.probability.crps_for_ensemble(
671 ens_mem,
672 ctrl,
673 ensemble_member_dim="realization",
674 method=method,
675 preserve_dims="time",
676 )
677 )
679 crps.rename(f"CRPS_of_{cubes[0].name()}")
680 _realization_callback(crps)
681 return crps
684def scores_pod_model_obs(
685 cubes: CubeList,
686 preserved_coordinates: list[str] | str | None,
687 threshold: str,
688 op_func: str,
689):
690 r"""
691 Compute the Probability of Detection (POD) score using Scores ([scoresa]_ [scoresb]_).
693 Parameters
694 ----------
695 cubes: iris.cube.CubeList
696 An iris cubelist containing model(s) and an observation cube.
697 preserved_coordinates: list | str | None
698 An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
699 then preserving coordinate 'time' will produce a probability of detection score for each timeslice (shape time). If None,
700 then it will return a single value score for all times/point locations.
701 threshold: str
702 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).
703 op_func: str
704 A string either containing 'lt' for less than or 'gt for greater than, to determine how the threshold is applied to the data
705 to generate the mask.
707 Returns
708 -------
709 cube: iris.cube
710 An iris cube, containing the probability of detection score for further plotting.
712 Notes
713 -----
714 The probability of detection calculates the proportion of observed events that meet a threshold that were correctly forecast by the model.
715 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
716 295k, that would be a positive hit. It does not take into account how far above/below a threshold a model forecasts.
718 It is calculated as .. math:: POD = \frac{true positives}{true positives + false negatives}
720 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.
722 POD produces a range of 0 to 1, where 1 is a perfect score.
723 """
724 # Split out model(s) and obs
725 models = CubeList()
726 for c in cubes:
727 if "observed" in c.long_name:
728 observed = c
729 else:
730 models.append(c)
732 # Setup cubelist to store results
733 scores_results = iris.cube.CubeList()
735 # Setup operators greater than, less than.
736 ops = {
737 "gt": operator.gt,
738 "lt": operator.lt,
739 }
741 try:
742 op = ops[op_func]
743 except KeyError as err:
744 raise ValueError(f"Operator {op_func} not supported.") from err
746 for model in models:
747 # Convert obs cubes to xarray and resolve preserved dimensions.
748 other_xr = xr.DataArray.from_iris(model)
749 base_xr = xr.DataArray.from_iris(observed)
750 preserve_dims = _resolve_preserve_dims(
751 observed, other_xr, preserved_coordinates
752 )
754 # Create event operator object using threshold and operator direction.
755 event_operator = scores.categorical.ThresholdEventOperator(
756 default_event_threshold=float(threshold), default_op_fn=op
757 )
759 # Generate binary fields using the event operator.
760 forecast_binary, observed_binary = event_operator.make_event_tables(
761 other_xr, base_xr
762 )
764 # Create binary contigency manager, as per Scores API, using transform to preserve preserve_dims
765 contingency_manager = scores.categorical.BinaryContingencyManager(
766 forecast_binary, observed_binary
767 ).transform(preserve_dims=preserve_dims)
769 # Get POD from the contigency manager, and convert back to an iris cube.
770 scores_cube = xr.DataArray.to_iris(
771 contingency_manager.probability_of_detection()
772 )
774 # Rename cube so it plots correctly alongside correcting cube units.
775 scores_cube.rename(
776 f"Probability_Of_Detection_{op_func}_{threshold}_{observed.name()}"
777 )
778 scores_cube.units = "1"
779 scores_cube.attributes["model_name"] = model.attributes["model_name"]
781 scores_results.append(scores_cube)
783 return scores_results
786def scores_ets_model_obs(
787 cubes: CubeList,
788 preserved_coordinates: list[str] | str | None,
789 threshold: str,
790 op_func: str,
791):
792 r"""
793 Compute the Equitable Threat Score (ETS) score using Scores ([scoresa]_ [scoresb]_).
795 Parameters
796 ----------
797 cubes: iris.cube.CubeList
798 An iris cubelist containing model(s) and an observation cube.
799 preserved_coordinates: list | str | None
800 An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
801 then preserving coordinate 'time' will produce the equitable threat score for each timeslice (shape time). If None,
802 then it will return a single value score for all times/point locations.
803 threshold: str
804 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).
805 op_func: str
806 A string either containing 'lt' for less than or 'gt for greater than, to determine how the threshold is applied to the data
807 to generate the mask.
809 Returns
810 -------
811 cube: iris.cube
812 An iris cube, containing the probability of detection score for further plotting.
814 Notes
815 -----
816 The Equitable Threat Score (ETS) evaluates the accuracy of forecasts for events that meet a specified threshold,
817 hile accounting for correct forecasts that could occur purely by chance. Unlike the Probability of Detection (POD),
818 ETS considers hits, misses, and false alarms, providing a more balanced assessment of forecast skill.
820 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
821 would be counted as a hit. ETS adjusts the total number of hits by removing the number of hits expected due to random chance.
823 It is calculated as:
825 .. math::
827 ETS = \frac{hits - hits_{random}}
828 {hits + misses + false\ alarms - hits_{random}}
830 where
832 hits_{random} = \frac{(hits + misses)(hits + false\ alarms)}{total count}
834 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
835 random chance.
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 ETS from the contigency manager, and convert back to an iris cube.
883 scores_cube = xr.DataArray.to_iris(contingency_manager.equitable_threat_score())
885 # Rename cube so it plots correctly alongside correcting cube units.
886 scores_cube.rename(
887 f"Equitable_Threat_Score_{op_func}_{threshold}_{observed.name()}"
888 )
889 scores_cube.units = "1"
890 scores_cube.attributes["model_name"] = model.attributes["model_name"]
892 scores_results.append(scores_cube)
894 return scores_results