Coverage for src/CSET/operators/scoreswrappers.py: 92%
229 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-17 11:00 +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_cubes_for_verification(cubes: CubeList):
45 """Prepare cubes ready for verification in scores.
47 Parameters
48 ----------
49 cubes: iris.cube.CubeList
50 A CubeList of exact 2 cubes, one from each model.
52 Returns
53 -------
54 base: iris.cube.Cube
55 The cube from the "analysis" in the same format as the other model.
56 other: iris.cube.Cube
57 The cube from the model in the same format as the base model.
59 Raises
60 ------
61 ValueError: "cubes should contain exactly 2 cubes."
62 If any other number of cubes are present.
64 Notes
65 -----
66 This operator is used for sorting the data into the correct format. It
67 is likely going to need to be refactored out of CSET and perhaps moved into
68 `CSET._utils` given common code between here and `misc.difference`.
69 """
70 # Set cubes into correct format using code from difference operator
71 if len(cubes) != 2:
72 raise ValueError("cubes should contain exactly 2 cubes.")
73 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
74 other: Cube = cubes.extract_cube(
75 iris.Constraint(
76 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
77 )
78 )
80 # If cubes contain a pressure coordinate, ensure it is increasing.
81 for cube in cubes:
82 try:
83 if len(cube.coord("pressure").points) > 2 and not is_increasing( 83 ↛ 86line 83 didn't jump to line 86 because the condition on line 83 was never true
84 cube.coord("pressure").points
85 ):
86 reverse(cube, "pressure")
88 except iris.exceptions.CoordinateNotFoundError:
89 pass
91 # Extract just common time points.
92 base, other = _extract_common_time_points(base, other)
94 # Get spatial coord names.
95 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
96 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
98 # Ensure cubes to compare are on common differencing grid.
99 # This is triggered if either
100 # i) latitude and longitude shapes are not the same. Note grid points
101 # are not compared directly as these can differ through rounding
102 # errors.
103 # ii) or variables are known to often sit on different grid staggering
104 # in different models (e.g. cell center vs cell edge), as is the case
105 # for UM and LFRic comparisons.
106 # In future greater choice of regridding method might be applied depending
107 # on variable type. Linear regridding can in general be appropriate for smooth
108 # variables. Care should be taken with interpretation of differences
109 # given this dependency on regridding.
110 if (
111 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
112 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
113 ) or (
114 base.long_name
115 in [
116 "eastward_wind_at_10m",
117 "northward_wind_at_10m",
118 "northward_wind_at_cell_centres",
119 "eastward_wind_at_cell_centres",
120 "zonal_wind_at_pressure_levels",
121 "meridional_wind_at_pressure_levels",
122 "potential_vorticity_at_pressure_levels",
123 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
124 ]
125 ):
126 logger.debug("Linear regridding base cube to other grid to compute differences")
127 base = regrid_onto_cube(base, other, method="Linear")
129 # Figure out if we are comparing between UM and LFRic; flip array if so.
130 base_lat_direction = is_increasing(base.coord(base_lat_name).points)
131 other_lat_direction = is_increasing(other.coord(other_lat_name).points)
132 if base_lat_direction != other_lat_direction: 132 ↛ 134line 132 didn't jump to line 134 because the condition on line 132 was never true
133 # Copy base cube for correct coordinate information.
134 other_tmp = base.copy()
135 # Flip the data and place in the copied cube.
136 other_tmp.data = np.flip(
137 other.data, other.coord(other_lat_name).cube_dims(other)
138 )
139 # Use original name and units from the other cube.
140 other_tmp.rename(other.name())
141 other_tmp.units = other.units
142 # Replace the cube.
143 other = other_tmp
145 # Equalise attributes so we can merge.
146 fully_equalise_attributes(CubeList([base, other]))
147 logger.debug("Base: %s\nOther: %s", base, other)
149 return base, other
152def _resolve_preserve_dims(
153 cube: Cube,
154 data_array: xr.DataArray,
155 preserved_coordinates: list[str] | str | None,
156) -> list[str] | None:
157 """Resolve preserve coordinates to xarray dimension names.
159 The ``scores`` package expects preserve dimensions to match xarray
160 dimension names. In Iris data, commonly used coordinates such as ``time``
161 may be auxiliary coordinates attached to a differently named dimension
162 (e.g. ``dim0``). This helper maps coordinate names to their underlying
163 dimension names and helps to convert from iris to xarray coordinate dimension names.
164 """
165 if preserved_coordinates is None:
166 return None
168 coord_names = (
169 [preserved_coordinates]
170 if isinstance(preserved_coordinates, str)
171 else preserved_coordinates
172 )
173 preserve_dims: list[str] = []
175 for coord_name in coord_names:
176 # Already an xarray dimension name.
177 if coord_name in data_array.dims:
178 if coord_name not in preserve_dims: 178 ↛ 180line 178 didn't jump to line 180 because the condition on line 178 was always true
179 preserve_dims.append(coord_name)
180 continue
182 # Otherwise, map coordinate name to dimension index/indices.
183 try:
184 dim_indices = cube.coord_dims(coord_name)
185 except iris.exceptions.CoordinateNotFoundError:
186 # Keep original name so scores raises a clear error for unknown keys.
187 if coord_name not in preserve_dims:
188 preserve_dims.append(coord_name)
189 continue
191 for dim_index in dim_indices:
192 dim_name = data_array.dims[dim_index]
193 if dim_name not in preserve_dims:
194 preserve_dims.append(dim_name)
196 return preserve_dims
199def scores_rmse_model_obs(
200 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
201):
202 r"""Calculate the Root Mean Square Error (RMSE) using scores.
204 Acts as a wrapper around the RMSE calculation from ``scores`` ([scoresa]_, [scoresb]_).
205 It is calculated as
207 .. math:: RMSE = \sqrt{\frac{1}{N} \Sigma(forecast - observations)^2}
209 Parameters
210 ----------
211 cubes: iris.cube.CubeList
212 A CubeList containing an observation cube and at least one model cube.
213 preserved_coordinates: list[str] | str | None, default is None.
214 The coordinates that you wish to preserve in the calculaiton of the
215 RMSE. For example if you want a map of each time you can preserve
216 ["time","grid_latitude", "grid_longitude"] or if you want a time series
217 you can preserve ["time"], if you want to collapse to a single value
218 use `None`. The default is `None`.
220 Returns
221 -------
222 scores_cube: iris.cube.Cube
223 A cube containing the RMSE between the models and observation cube.
224 """
225 rmse_cubes = CubeList()
226 model_list = CubeList()
228 for cb in cubes:
229 if "observed" in cb.long_name:
230 observed = cb
231 else:
232 model_list.append(cb)
234 for model in model_list:
235 input_cubelist = CubeList()
236 input_cubelist.append(observed)
237 input_cubelist.append(model)
238 rmse = scores_rmse(
239 input_cubelist, preserved_coordinates, obs_model_comparison=True
240 )
241 model_name = model.attributes["model_name"]
242 rmse.attributes["model_name"] = model_name
243 rmse_cubes.append(rmse)
245 return rmse_cubes
248def scores_mae_model_obs(
249 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
250):
251 r"""Calculate the Mean Absolute Error (MAE) using scores.
253 Acts as a wrapper around the MAE calculation from ``scores`` ([scoresa]_, [scoresb]_).
255 Parameters
256 ----------
257 cubes: iris.cube.CubeList
258 A CubeList containing an observation cube and at least one model cube.
259 preserved_coordinates: list[str] | str | None, default is None.
260 The coordinates that you wish to preserve in the calculaiton of the
261 MAE. For example if you want a map of each time you can preserve
262 ["time","latitude", "longitude"] or if you want a time series
263 you can preserve ["time"], if you want to collapse to a single value
264 use `None`. The default is `None`.
266 Returns
267 -------
268 scores_cube: iris.cube.Cube
269 A cube containing the MAE between the models and observation cube.
270 """
271 mae_cubes = CubeList()
272 model_list = CubeList()
274 for cb in cubes:
275 if "observed" in cb.long_name:
276 observed = cb
277 else:
278 model_list.append(cb)
280 for model in model_list:
281 input_cubelist = CubeList()
282 input_cubelist.append(observed)
283 input_cubelist.append(model)
284 mae = scores_mae(
285 input_cubelist, preserved_coordinates, obs_model_comparison=True
286 )
287 model_name = model.attributes["model_name"]
288 mae.attributes["model_name"] = model_name
289 mae_cubes.append(mae)
291 return mae_cubes
294def scores_rmse(
295 cubes: CubeList,
296 preserved_coordinates: list[str] | str | None = None,
297 obs_model_comparison: bool = False,
298):
299 r"""Calculate the Root Mean Square Error (RMSE) using scores.
301 Acts as a wrapper around the RMSE calculation from ``scores`` ([scoresa]_, [scoresb]_).
302 It is calculated as
304 .. math:: RMSE = \sqrt{\frac{1}{N} \Sigma(forecast - observations)^2}
306 Parameters
307 ----------
308 cubes: iris.cube.CubeList
309 A CubeList containing exactly two cubes: a base and an "other" model,
310 this can be an analysis and the model.
311 preserved_coordinates: list[str] | str | None, default is None.
312 The coordinates (or xarray dimension names) that you wish to preserve in the calculaiton of the
313 RMSE. For example if you want a map of each time you can preserve
314 ["time","grid_latitude", "grid_longitude"] or if you want a time series
315 you can preserve ["time"], if you want to collapse to a single value
316 use `None`. The default is `None`.
318 Returns
319 -------
320 scores_cube: iris.cube.Cube
321 A cube containing the RMSE between the base and other cube.
322 """
323 if obs_model_comparison:
324 for cb in cubes:
325 if "observed" in cb.long_name:
326 base = cb
327 else:
328 other = cb
329 else:
330 base, other = _sort_cubes_for_verification(cubes)
332 # Copy the coordinates of the input cubes.
333 other_xr = xr.DataArray.from_iris(other)
334 base_xr = xr.DataArray.from_iris(base)
335 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
337 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
338 # apply scores, and then transform it back.
339 scores_cube = xr.DataArray.to_iris(
340 scores.continuous.rmse(
341 other_xr,
342 base_xr,
343 preserve_dims=preserve_dims,
344 )
345 )
347 # If time is aggregated out, attach a scalar time coordinate with bounds
348 # so plotting can display the aggregated period in the title.
349 try:
350 if not scores_cube.coords("time"):
351 base_time = base.coord("time")
352 time_vals = (
353 base_time.bounds.flatten()
354 if base_time.has_bounds()
355 else base_time.points
356 )
357 t_start = float(time_vals[0])
358 t_end = float(time_vals[-1])
359 t_mid = 0.5 * (t_start + t_end)
361 scores_cube.add_aux_coord(
362 iris.coords.AuxCoord(
363 t_mid,
364 standard_name=base_time.standard_name,
365 long_name=base_time.long_name,
366 var_name=base_time.var_name,
367 units=base_time.units,
368 bounds=np.array([t_start, t_end]),
369 attributes=base_time.attributes.copy(),
370 )
371 )
372 except iris.exceptions.CoordinateNotFoundError:
373 pass
375 scores_cube.rename(f"RMSE_of_{base.name()}")
376 # if preserved_coordinates == ["grid_latitude", "grid_longitude"]:
377 # scores_cube.add_aux_coord(time_coord)
378 return scores_cube
381def scores_mae(
382 cubes: CubeList,
383 preserved_coordinates: list[str] | str | None = None,
384 obs_model_comparison: bool = False,
385):
386 r"""Calculate the Mean Absolute Error (MAE) using scores.
388 Acts as a wrapper around the MAE calculation from ``scores`` ([scoresa]_, [scoresb]_).
390 Parameters
391 ----------
392 cubes: iris.cube.CubeList
393 A CubeList containing exactly two cubes: a base and an "other" model,
394 this can be an analysis and the model.
395 preserved_coordinates: list[str] | str | None, default is None.
396 The coordinates that you wish to preserve in the calculaiton of the
397 MAE. For example if you want a map of each time you can preserve
398 ["time","grid_latitude", "grid_longitude"] or if you want a time series
399 you can preserve ["time"], if you want to collapse to a single value
400 use `None`. The default is `None`.
402 Returns
403 -------
404 scores_cube: iris.cube.Cube
405 A cube containing the MAE between the base and other cube.
406 """
407 if obs_model_comparison:
408 for cb in cubes:
409 if "observed" in cb.long_name:
410 base = cb
411 else:
412 other = cb
413 else:
414 base, other = _sort_cubes_for_verification(cubes)
416 # Copy the coordinates of the input cubes.
417 other_xr = xr.DataArray.from_iris(other)
418 base_xr = xr.DataArray.from_iris(base)
419 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
421 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
422 # apply scores, and then transform it back.
423 scores_cube = xr.DataArray.to_iris(
424 scores.continuous.mae(
425 other_xr,
426 base_xr,
427 preserve_dims=preserve_dims,
428 )
429 )
431 # If time is aggregated out, attach a scalar time coordinate with bounds
432 # so plotting can display the aggregated period in the title.
433 try:
434 if not scores_cube.coords("time"):
435 base_time = base.coord("time")
436 time_vals = (
437 base_time.bounds.flatten()
438 if base_time.has_bounds()
439 else base_time.points
440 )
441 t_start = float(time_vals[0])
442 t_end = float(time_vals[-1])
443 t_mid = 0.5 * (t_start + t_end)
445 scores_cube.add_aux_coord(
446 iris.coords.AuxCoord(
447 t_mid,
448 standard_name=base_time.standard_name,
449 long_name=base_time.long_name,
450 var_name=base_time.var_name,
451 units=base_time.units,
452 bounds=np.array([t_start, t_end]),
453 attributes=base_time.attributes.copy(),
454 )
455 )
456 except iris.exceptions.CoordinateNotFoundError:
457 pass
459 scores_cube.rename(f"MAE_of_{base.name()}")
460 return scores_cube
463def scores_additive_bias(
464 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
465):
466 r"""Calculate the Additive Bias (Mean Error) using scores.
468 Acts as a wrapper around the ME calculation from ``scores`` ([scoresa]_, [scoresb]_).
470 Parameters
471 ----------
472 cubes: iris.cube.CubeList
473 A CubeList containing exactly two cubes: a base and an "other" model,
474 this can be an analysis and the model.
475 preserved_coordinates: list[str] | str | None, default is None.
476 The coordinates that you wish to preserve in the calculaiton of the
477 ME. For example if you want a map of each time you can preserve
478 ["time","grid_latitude", "grid_longitude"] or if you want a time series
479 you can preserve ["time"], if you want to collapse to a single value
480 use `None`. The default is `None`.
482 Returns
483 -------
484 scores_cube: iris.cube.Cube
485 A cube containing the ME between the base and other cube.
486 """
487 base, other = _sort_cubes_for_verification(cubes)
489 # Copy the coordinates of the input cubes.
490 other_xr = xr.DataArray.from_iris(other)
491 base_xr = xr.DataArray.from_iris(base)
492 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
494 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
495 # apply scores, and then transform it back.
496 scores_cube = xr.DataArray.to_iris(
497 scores.continuous.additive_bias(
498 other_xr,
499 base_xr,
500 preserve_dims=preserve_dims,
501 )
502 )
504 # If time is aggregated out, attach a scalar time coordinate with bounds
505 # so plotting can display the aggregated period in the title.
506 try:
507 if not scores_cube.coords("time"): 507 ↛ 531line 507 didn't jump to line 531 because the condition on line 507 was always true
508 base_time = base.coord("time")
509 time_vals = (
510 base_time.bounds.flatten()
511 if base_time.has_bounds()
512 else base_time.points
513 )
514 t_start = float(time_vals[0])
515 t_end = float(time_vals[-1])
516 t_mid = 0.5 * (t_start + t_end)
518 scores_cube.add_aux_coord(
519 iris.coords.AuxCoord(
520 t_mid,
521 standard_name=base_time.standard_name,
522 long_name=base_time.long_name,
523 var_name=base_time.var_name,
524 units=base_time.units,
525 bounds=np.array([t_start, t_end]),
526 attributes=base_time.attributes.copy(),
527 )
528 )
529 except iris.exceptions.CoordinateNotFoundError:
530 pass
531 scores_cube.rename(f"Additive_Bias_of_{base.name()}")
532 return scores_cube
535def scores_correlation_pearsonr(
536 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
537):
538 r"""Calculate the Pearson's Correlation (PC) coefficient using scores.
540 Acts as a wrapper around the PC calculation from ``scores`` ([scoresa]_, [scoresb]_).
542 Parameters
543 ----------
544 cubes: iris.cube.CubeList
545 A CubeList containing exactly two cubes: a base and an "other" model,
546 this can be an analysis and the model.
547 preserved_coordinates: list[str] | str | None, default is None.
548 The coordinates that you wish to preserve in the calculation of the
549 PC. For example if you want a map of each time you can preserve
550 ["time","grid_latitude", "grid_longitude"] or if you want a time series
551 you can preserve ["time"], if you want to collapse to a single value
552 use `None`. The default is `None`.
554 Returns
555 -------
556 scores_cube: iris.cube.Cube
557 A cube containing the PC between the base and other cube.
558 """
559 base, other = _sort_cubes_for_verification(cubes)
561 # Copy the coordinates of the input cubes.
562 other_xr = xr.DataArray.from_iris(other)
563 base_xr = xr.DataArray.from_iris(base)
564 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
566 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
567 # apply scores, and then transform it back.
568 scores_cube = xr.DataArray.to_iris(
569 scores.continuous.correlation.pearsonr(
570 other_xr,
571 base_xr,
572 preserve_dims=preserve_dims,
573 )
574 )
576 # If time is aggregated out, attach a scalar time coordinate with bounds
577 # so plotting can display the aggregated period in the title.
578 try:
579 if not scores_cube.coords("time"): 579 ↛ 604line 579 didn't jump to line 604 because the condition on line 579 was always true
580 base_time = base.coord("time")
581 time_vals = (
582 base_time.bounds.flatten()
583 if base_time.has_bounds()
584 else base_time.points
585 )
586 t_start = float(time_vals[0])
587 t_end = float(time_vals[-1])
588 t_mid = 0.5 * (t_start + t_end)
590 scores_cube.add_aux_coord(
591 iris.coords.AuxCoord(
592 t_mid,
593 standard_name=base_time.standard_name,
594 long_name=base_time.long_name,
595 var_name=base_time.var_name,
596 units=base_time.units,
597 bounds=np.array([t_start, t_end]),
598 attributes=base_time.attributes.copy(),
599 )
600 )
601 except iris.exceptions.CoordinateNotFoundError:
602 pass
604 scores_cube.rename(f"Pearson_Correlation_of_{base.name()}")
605 return scores_cube
608def scores_crps_for_ensemble(
609 cubes: Cube | CubeList, method: str = "ecdf", control_member: int = 0
610) -> iris.Constraint:
611 r"""Calculate the CRPS for an ensemble.
613 Acts as a wrapper around the crps_for_ensemble from ``scores`` ([scoresa]_, [scoresb]_).
615 Lower CRPS values are better (implies experiment distribution is closer to control distribution/observations),
616 larger values are worse (implies distributions are dissimilar).
617 It is applicable across time and spatial scales as the focus is on the distribution of the values.
618 Default method is ecdf. ecdf is exact value from the empirical distributions,
619 whereas fair produces an approximated value based on a random sample of the underlying distribution.
621 See [CRPS]_ for further information.
623 Parameters
624 ----------
625 cubes: iris.cube.Cube
626 A Cube containing ensembles data
628 Returns
629 -------
630 crps: iris.cube.Cube
631 A cube containing the crps between the ensemble members and the control
632 """
633 if control_member != 0:
634 logger.warning("control member is usual 0")
636 if control_member not in cubes.coords("realization")[0].points:
637 new_control_member = cubes.coords("realization")[0].points[0]
638 logger.warning(
639 f"control member value {control_member} out of bounds, defaulting to control member={new_control_member}"
640 )
641 control_member = new_control_member
643 if cubes.coord("time").shape[0] == 1:
644 raise ValueError("Cube has only one time point.")
646 if cubes.coord("realization").shape[0] < 3:
647 raise ValueError("Cube should have one control member and at least two members")
649 ctrl = cubes.extract(generate_realization_constraint([control_member]))
650 ens_mem = cubes.extract(
651 generate_remove_single_ensemble_member_constraint(control_member)
652 )
654 # Realising the data in advance provides a large speedup
655 _ = ctrl.data
656 _ = ens_mem.data
657 del _
659 ctrl = xr.DataArray.from_iris(ctrl)
660 ens_mem = xr.DataArray.from_iris(ens_mem)
662 crps = xr.DataArray.to_iris(
663 scores.probability.crps_for_ensemble(
664 ens_mem,
665 ctrl,
666 ensemble_member_dim="realization",
667 method=method,
668 preserve_dims="time",
669 )
670 )
672 crps.rename(f"CRPS_of_{cubes[0].name()}")
673 _realization_callback(crps)
674 return crps
677def scores_pod_model_obs(
678 cubes: CubeList,
679 preserved_coordinates: list[str] | str | None,
680 threshold: str,
681 op_func: str,
682):
683 r"""
684 Compute the Probability of Detection (POD) score using Scores ([scoresa]_ [scoresb]_).
686 Parameters
687 ----------
688 cubes: iris.cube.CubeList
689 An iris cubelist containing model(s) and an observation cube.
690 preserved_coordinates: list | str | None
691 An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
692 then preserving coordinate 'time' will produce a probability of detection score for each timeslice (shape time). If None,
693 then it will return a single value score for all times/point locations.
694 threshold: str
695 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).
696 op_func: str
697 A string either containing 'lt' for less than or 'gt for greater than, to determine how the threshold is applied to the data
698 to generate the mask.
700 Returns
701 -------
702 cube: iris.cube
703 An iris cube, containing the probability of detection score for further plotting.
705 Notes
706 -----
707 The probability of detection calculates the proportion of observed events that meet a threshold that were correctly forecast by the model.
708 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
709 295k, that would be a positive hit. It does not take into account how far above/below a threshold a model forecasts.
711 It is calculated as .. math:: POD = \frac{true positives}{true positives + false negatives}
713 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.
715 POD produces a range of 0 to 1, where 1 is a perfect score.
716 """
717 # Split out model(s) and obs
718 models = CubeList()
719 for c in cubes:
720 if "observed" in c.long_name:
721 observed = c
722 else:
723 models.append(c)
725 # Setup cubelist to store results
726 scores_results = iris.cube.CubeList()
728 # Setup operators greater than, less than.
729 ops = {
730 "gt": operator.gt,
731 "lt": operator.lt,
732 }
734 try:
735 op = ops[op_func]
736 except KeyError as err:
737 raise ValueError(f"Operator {op_func} not supported.") from err
739 for model in models:
740 # Convert obs cubes to xarray and resolve preserved dimensions.
741 other_xr = xr.DataArray.from_iris(model)
742 base_xr = xr.DataArray.from_iris(observed)
743 preserve_dims = _resolve_preserve_dims(
744 observed, other_xr, preserved_coordinates
745 )
747 # Create event operator object using threshold and operator direction.
748 event_operator = scores.categorical.ThresholdEventOperator(
749 default_event_threshold=float(threshold), default_op_fn=op
750 )
752 # Generate binary fields using the event operator.
753 forecast_binary, observed_binary = event_operator.make_event_tables(
754 other_xr, base_xr
755 )
757 # Create binary contigency manager, as per Scores API, using transform to preserve preserve_dims
758 contingency_manager = scores.categorical.BinaryContingencyManager(
759 forecast_binary, observed_binary
760 ).transform(preserve_dims=preserve_dims)
762 # Get POD from the contigency manager, and convert back to an iris cube.
763 scores_cube = xr.DataArray.to_iris(
764 contingency_manager.probability_of_detection()
765 )
767 # Rename cube so it plots correctly alongside correcting cube units.
768 scores_cube.rename(
769 f"Probability_Of_Detection_{op_func}_{threshold}_{observed.name()}"
770 )
771 scores_cube.units = "1"
772 scores_cube.attributes["model_name"] = model.attributes["model_name"]
774 scores_results.append(scores_cube)
776 return scores_results