Coverage for src/CSET/operators/scoreswrappers.py: 92%
229 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-17 13:24 +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_additive_bias_model_obs(
249 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
250):
251 r"""Calculate the Additive Bias (Mean Error) using scores.
253 Acts as a wrapper around the ME 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 ME. 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.CubeList
269 A cube list containing the ME between the models and observation cube.
270 """
271 additive_bias_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 additive_bias = scores_additive_bias(
285 input_cubelist, preserved_coordinates, obs_model_comparison=True
286 )
287 model_name = model.attributes["model_name"]
288 additive_bias.attributes["model_name"] = model_name
289 additive_bias_cubes.append(additive_bias)
291 return additive_bias_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(cubes: CubeList, preserved_coordinates: list[str] | str | None = None):
382 r"""Calculate the Mean Absolute Error (MAE) using scores.
384 Acts as a wrapper around the MAE calculation from ``scores`` ([scoresa]_, [scoresb]_).
386 Parameters
387 ----------
388 cubes: iris.cube.CubeList
389 A CubeList containing exactly two cubes: a base and an "other" model,
390 this can be an analysis and the model.
391 preserved_coordinates: list[str] | str | None, default is None.
392 The coordinates that you wish to preserve in the calculaiton of the
393 MAE. For example if you want a map of each time you can preserve
394 ["time","grid_latitude", "grid_longitude"] or if you want a time series
395 you can preserve ["time"], if you want to collapse to a single value
396 use `None`. The default is `None`.
398 Returns
399 -------
400 scores_cube: iris.cube.Cube
401 A cube containing the MAE between the base and other cube.
402 """
403 base, other = _sort_cubes_for_verification(cubes)
405 # Copy the coordinates of the input cubes.
406 other_xr = xr.DataArray.from_iris(other)
407 base_xr = xr.DataArray.from_iris(base)
408 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
410 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
411 # apply scores, and then transform it back.
412 scores_cube = xr.DataArray.to_iris(
413 scores.continuous.mae(
414 other_xr,
415 base_xr,
416 preserve_dims=preserve_dims,
417 )
418 )
420 # If time is aggregated out, attach a scalar time coordinate with bounds
421 # so plotting can display the aggregated period in the title.
422 try:
423 if not scores_cube.coords("time"): 423 ↛ 448line 423 didn't jump to line 448 because the condition on line 423 was always true
424 base_time = base.coord("time")
425 time_vals = (
426 base_time.bounds.flatten()
427 if base_time.has_bounds()
428 else base_time.points
429 )
430 t_start = float(time_vals[0])
431 t_end = float(time_vals[-1])
432 t_mid = 0.5 * (t_start + t_end)
434 scores_cube.add_aux_coord(
435 iris.coords.AuxCoord(
436 t_mid,
437 standard_name=base_time.standard_name,
438 long_name=base_time.long_name,
439 var_name=base_time.var_name,
440 units=base_time.units,
441 bounds=np.array([t_start, t_end]),
442 attributes=base_time.attributes.copy(),
443 )
444 )
445 except iris.exceptions.CoordinateNotFoundError:
446 pass
448 scores_cube.rename(f"MAE_of_{base.name()}")
449 return scores_cube
452def scores_additive_bias(
453 cubes: CubeList,
454 preserved_coordinates: list[str] | str | None = None,
455 obs_model_comparison: bool = False,
456):
457 r"""Calculate the Additive Bias (Mean Error) using scores.
459 Acts as a wrapper around the ME calculation from ``scores`` ([scoresa]_, [scoresb]_).
461 Parameters
462 ----------
463 cubes: iris.cube.CubeList
464 A CubeList containing exactly two cubes: a base and an "other" model,
465 this can be an analysis and the model.
466 preserved_coordinates: list[str] | str | None, default is None.
467 The coordinates that you wish to preserve in the calculaiton of the
468 ME. For example if you want a map of each time you can preserve
469 ["time","grid_latitude", "grid_longitude"] or if you want a time series
470 you can preserve ["time"], if you want to collapse to a single value
471 use `None`. The default is `None`.
473 Returns
474 -------
475 scores_cube: iris.cube.Cube
476 A cube containing the ME between the base and other cube.
477 """
478 if obs_model_comparison:
479 for cb in cubes:
480 if "observed" in cb.long_name:
481 base = cb
482 else:
483 other = cb
484 else:
485 base, other = _sort_cubes_for_verification(cubes)
487 # Copy the coordinates of the input cubes.
488 other_xr = xr.DataArray.from_iris(other)
489 base_xr = xr.DataArray.from_iris(base)
490 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
492 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
493 # apply scores, and then transform it back.
494 scores_cube = xr.DataArray.to_iris(
495 scores.continuous.additive_bias(
496 other_xr,
497 base_xr,
498 preserve_dims=preserve_dims,
499 )
500 )
502 # If time is aggregated out, attach a scalar time coordinate with bounds
503 # so plotting can display the aggregated period in the title.
504 try:
505 if not scores_cube.coords("time"):
506 base_time = base.coord("time")
507 time_vals = (
508 base_time.bounds.flatten()
509 if base_time.has_bounds()
510 else base_time.points
511 )
512 t_start = float(time_vals[0])
513 t_end = float(time_vals[-1])
514 t_mid = 0.5 * (t_start + t_end)
516 scores_cube.add_aux_coord(
517 iris.coords.AuxCoord(
518 t_mid,
519 standard_name=base_time.standard_name,
520 long_name=base_time.long_name,
521 var_name=base_time.var_name,
522 units=base_time.units,
523 bounds=np.array([t_start, t_end]),
524 attributes=base_time.attributes.copy(),
525 )
526 )
527 except iris.exceptions.CoordinateNotFoundError:
528 pass
529 scores_cube.rename(f"Additive_Bias_of_{base.name()}")
530 return scores_cube
533def scores_correlation_pearsonr(
534 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
535):
536 r"""Calculate the Pearson's Correlation (PC) coefficient using scores.
538 Acts as a wrapper around the PC calculation from ``scores`` ([scoresa]_, [scoresb]_).
540 Parameters
541 ----------
542 cubes: iris.cube.CubeList
543 A CubeList containing exactly two cubes: a base and an "other" model,
544 this can be an analysis and the model.
545 preserved_coordinates: list[str] | str | None, default is None.
546 The coordinates that you wish to preserve in the calculation of the
547 PC. For example if you want a map of each time you can preserve
548 ["time","grid_latitude", "grid_longitude"] or if you want a time series
549 you can preserve ["time"], if you want to collapse to a single value
550 use `None`. The default is `None`.
552 Returns
553 -------
554 scores_cube: iris.cube.Cube
555 A cube containing the PC between the base and other cube.
556 """
557 base, other = _sort_cubes_for_verification(cubes)
559 # Copy the coordinates of the input cubes.
560 other_xr = xr.DataArray.from_iris(other)
561 base_xr = xr.DataArray.from_iris(base)
562 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
564 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
565 # apply scores, and then transform it back.
566 scores_cube = xr.DataArray.to_iris(
567 scores.continuous.correlation.pearsonr(
568 other_xr,
569 base_xr,
570 preserve_dims=preserve_dims,
571 )
572 )
574 # If time is aggregated out, attach a scalar time coordinate with bounds
575 # so plotting can display the aggregated period in the title.
576 try:
577 if not scores_cube.coords("time"): 577 ↛ 602line 577 didn't jump to line 602 because the condition on line 577 was always true
578 base_time = base.coord("time")
579 time_vals = (
580 base_time.bounds.flatten()
581 if base_time.has_bounds()
582 else base_time.points
583 )
584 t_start = float(time_vals[0])
585 t_end = float(time_vals[-1])
586 t_mid = 0.5 * (t_start + t_end)
588 scores_cube.add_aux_coord(
589 iris.coords.AuxCoord(
590 t_mid,
591 standard_name=base_time.standard_name,
592 long_name=base_time.long_name,
593 var_name=base_time.var_name,
594 units=base_time.units,
595 bounds=np.array([t_start, t_end]),
596 attributes=base_time.attributes.copy(),
597 )
598 )
599 except iris.exceptions.CoordinateNotFoundError:
600 pass
602 scores_cube.rename(f"Pearson_Correlation_of_{base.name()}")
603 return scores_cube
606def scores_crps_for_ensemble(
607 cubes: Cube | CubeList, method: str = "ecdf", control_member: int = 0
608) -> iris.Constraint:
609 r"""Calculate the CRPS for an ensemble.
611 Acts as a wrapper around the crps_for_ensemble from ``scores`` ([scoresa]_, [scoresb]_).
613 Lower CRPS values are better (implies experiment distribution is closer to control distribution/observations),
614 larger values are worse (implies distributions are dissimilar).
615 It is applicable across time and spatial scales as the focus is on the distribution of the values.
616 Default method is ecdf. ecdf is exact value from the empirical distributions,
617 whereas fair produces an approximated value based on a random sample of the underlying distribution.
619 See [CRPS]_ for further information.
621 Parameters
622 ----------
623 cubes: iris.cube.Cube
624 A Cube containing ensembles data
626 Returns
627 -------
628 crps: iris.cube.Cube
629 A cube containing the crps between the ensemble members and the control
630 """
631 if control_member != 0:
632 logger.warning("control member is usual 0")
634 if control_member not in cubes.coords("realization")[0].points:
635 new_control_member = cubes.coords("realization")[0].points[0]
636 logger.warning(
637 f"control member value {control_member} out of bounds, defaulting to control member={new_control_member}"
638 )
639 control_member = new_control_member
641 if cubes.coord("time").shape[0] == 1:
642 raise ValueError("Cube has only one time point.")
644 if cubes.coord("realization").shape[0] < 3:
645 raise ValueError("Cube should have one control member and at least two members")
647 ctrl = cubes.extract(generate_realization_constraint([control_member]))
648 ens_mem = cubes.extract(
649 generate_remove_single_ensemble_member_constraint(control_member)
650 )
652 # Realising the data in advance provides a large speedup
653 _ = ctrl.data
654 _ = ens_mem.data
655 del _
657 ctrl = xr.DataArray.from_iris(ctrl)
658 ens_mem = xr.DataArray.from_iris(ens_mem)
660 crps = xr.DataArray.to_iris(
661 scores.probability.crps_for_ensemble(
662 ens_mem,
663 ctrl,
664 ensemble_member_dim="realization",
665 method=method,
666 preserve_dims="time",
667 )
668 )
670 crps.rename(f"CRPS_of_{cubes[0].name()}")
671 _realization_callback(crps)
672 return crps
675def scores_pod_model_obs(
676 cubes: CubeList,
677 preserved_coordinates: list[str] | str | None,
678 threshold: str,
679 op_func: str,
680):
681 r"""
682 Compute the Probability of Detection (POD) score using Scores ([scoresa]_ [scoresb]_).
684 Parameters
685 ----------
686 cubes: iris.cube.CubeList
687 An iris cubelist containing model(s) and an observation cube.
688 preserved_coordinates: list | str | None
689 An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
690 then preserving coordinate 'time' will produce a probability of detection score for each timeslice (shape time). If None,
691 then it will return a single value score for all times/point locations.
692 threshold: str
693 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).
694 op_func: str
695 A string either containing 'lt' for less than or 'gt for greater than, to determine how the threshold is applied to the data
696 to generate the mask.
698 Returns
699 -------
700 cube: iris.cube
701 An iris cube, containing the probability of detection score for further plotting.
703 Notes
704 -----
705 The probability of detection calculates the proportion of observed events that meet a threshold that were correctly forecast by the model.
706 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
707 295k, that would be a positive hit. It does not take into account how far above/below a threshold a model forecasts.
709 It is calculated as .. math:: POD = \frac{true positives}{true positives + false negatives}
711 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.
713 POD produces a range of 0 to 1, where 1 is a perfect score.
714 """
715 # Split out model(s) and obs
716 models = CubeList()
717 for c in cubes:
718 if "observed" in c.long_name:
719 observed = c
720 else:
721 models.append(c)
723 # Setup cubelist to store results
724 scores_results = iris.cube.CubeList()
726 # Setup operators greater than, less than.
727 ops = {
728 "gt": operator.gt,
729 "lt": operator.lt,
730 }
732 try:
733 op = ops[op_func]
734 except KeyError as err:
735 raise ValueError(f"Operator {op_func} not supported.") from err
737 for model in models:
738 # Convert obs cubes to xarray and resolve preserved dimensions.
739 other_xr = xr.DataArray.from_iris(model)
740 base_xr = xr.DataArray.from_iris(observed)
741 preserve_dims = _resolve_preserve_dims(
742 observed, other_xr, preserved_coordinates
743 )
745 # Create event operator object using threshold and operator direction.
746 event_operator = scores.categorical.ThresholdEventOperator(
747 default_event_threshold=float(threshold), default_op_fn=op
748 )
750 # Generate binary fields using the event operator.
751 forecast_binary, observed_binary = event_operator.make_event_tables(
752 other_xr, base_xr
753 )
755 # Create binary contigency manager, as per Scores API, using transform to preserve preserve_dims
756 contingency_manager = scores.categorical.BinaryContingencyManager(
757 forecast_binary, observed_binary
758 ).transform(preserve_dims=preserve_dims)
760 # Get POD from the contigency manager, and convert back to an iris cube.
761 scores_cube = xr.DataArray.to_iris(
762 contingency_manager.probability_of_detection()
763 )
765 # Rename cube so it plots correctly alongside correcting cube units.
766 scores_cube.rename(
767 f"Probability_Of_Detection_{op_func}_{threshold}_{observed.name()}"
768 )
769 scores_cube.units = "1"
770 scores_cube.attributes["model_name"] = model.attributes["model_name"]
772 scores_results.append(scores_cube)
774 return scores_results