Coverage for src/CSET/operators/scoreswrappers.py: 91%
208 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-20 10:04 +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_rmse(
249 cubes: CubeList,
250 preserved_coordinates: list[str] | str | None = None,
251 obs_model_comparison: bool = False,
252):
253 r"""Calculate the Root Mean Square Error (RMSE) using scores.
255 Acts as a wrapper around the RMSE calculation from ``scores`` ([scoresa]_, [scoresb]_).
256 It is calculated as
258 .. math:: RMSE = \sqrt{\frac{1}{N} \Sigma(forecast - observations)^2}
260 Parameters
261 ----------
262 cubes: iris.cube.CubeList
263 A CubeList containing exactly two cubes: a base and an "other" model,
264 this can be an analysis and the model.
265 preserved_coordinates: list[str] | str | None, default is None.
266 The coordinates (or xarray dimension names) that you wish to preserve in the calculaiton of the
267 RMSE. For example if you want a map of each time you can preserve
268 ["time","grid_latitude", "grid_longitude"] or if you want a time series
269 you can preserve ["time"], if you want to collapse to a single value
270 use `None`. The default is `None`.
272 Returns
273 -------
274 scores_cube: iris.cube.Cube
275 A cube containing the RMSE between the base and other cube.
276 """
277 if obs_model_comparison:
278 for cb in cubes:
279 if "observed" in cb.long_name:
280 base = cb
281 else:
282 other = cb
283 else:
284 base, other = _sort_cubes_for_verification(cubes)
286 # Copy the coordinates of the input cubes.
287 other_xr = xr.DataArray.from_iris(other)
288 base_xr = xr.DataArray.from_iris(base)
289 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
291 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
292 # apply scores, and then transform it back.
293 scores_cube = xr.DataArray.to_iris(
294 scores.continuous.rmse(
295 other_xr,
296 base_xr,
297 preserve_dims=preserve_dims,
298 )
299 )
301 # If time is aggregated out, attach a scalar time coordinate with bounds
302 # so plotting can display the aggregated period in the title.
303 try:
304 if not scores_cube.coords("time"):
305 base_time = base.coord("time")
306 time_vals = (
307 base_time.bounds.flatten()
308 if base_time.has_bounds()
309 else base_time.points
310 )
311 t_start = float(time_vals[0])
312 t_end = float(time_vals[-1])
313 t_mid = 0.5 * (t_start + t_end)
315 scores_cube.add_aux_coord(
316 iris.coords.AuxCoord(
317 t_mid,
318 standard_name=base_time.standard_name,
319 long_name=base_time.long_name,
320 var_name=base_time.var_name,
321 units=base_time.units,
322 bounds=np.array([t_start, t_end]),
323 attributes=base_time.attributes.copy(),
324 )
325 )
326 except iris.exceptions.CoordinateNotFoundError:
327 pass
329 scores_cube.rename(f"RMSE_of_{base.name()}")
330 # if preserved_coordinates == ["grid_latitude", "grid_longitude"]:
331 # scores_cube.add_aux_coord(time_coord)
332 return scores_cube
335def scores_mae(cubes: CubeList, preserved_coordinates: list[str] | str | None = None):
336 r"""Calculate the Mean Absolute Error (MAE) using scores.
338 Acts as a wrapper around the MAE calculation from ``scores`` ([scoresa]_, [scoresb]_).
340 Parameters
341 ----------
342 cubes: iris.cube.CubeList
343 A CubeList containing exactly two cubes: a base and an "other" model,
344 this can be an analysis and the model.
345 preserved_coordinates: list[str] | str | None, default is None.
346 The coordinates that you wish to preserve in the calculaiton of the
347 MAE. For example if you want a map of each time you can preserve
348 ["time","grid_latitude", "grid_longitude"] or if you want a time series
349 you can preserve ["time"], if you want to collapse to a single value
350 use `None`. The default is `None`.
352 Returns
353 -------
354 scores_cube: iris.cube.Cube
355 A cube containing the MAE between the base and other cube.
356 """
357 base, other = _sort_cubes_for_verification(cubes)
359 # Copy the coordinates of the input cubes.
360 other_xr = xr.DataArray.from_iris(other)
361 base_xr = xr.DataArray.from_iris(base)
362 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
364 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
365 # apply scores, and then transform it back.
366 scores_cube = xr.DataArray.to_iris(
367 scores.continuous.mae(
368 other_xr,
369 base_xr,
370 preserve_dims=preserve_dims,
371 )
372 )
374 # If time is aggregated out, attach a scalar time coordinate with bounds
375 # so plotting can display the aggregated period in the title.
376 try:
377 if not scores_cube.coords("time"): 377 ↛ 402line 377 didn't jump to line 402 because the condition on line 377 was always true
378 base_time = base.coord("time")
379 time_vals = (
380 base_time.bounds.flatten()
381 if base_time.has_bounds()
382 else base_time.points
383 )
384 t_start = float(time_vals[0])
385 t_end = float(time_vals[-1])
386 t_mid = 0.5 * (t_start + t_end)
388 scores_cube.add_aux_coord(
389 iris.coords.AuxCoord(
390 t_mid,
391 standard_name=base_time.standard_name,
392 long_name=base_time.long_name,
393 var_name=base_time.var_name,
394 units=base_time.units,
395 bounds=np.array([t_start, t_end]),
396 attributes=base_time.attributes.copy(),
397 )
398 )
399 except iris.exceptions.CoordinateNotFoundError:
400 pass
402 scores_cube.rename(f"MAE_of_{base.name()}")
403 return scores_cube
406def scores_additive_bias(
407 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
408):
409 r"""Calculate the Additive Bias (Mean Error) using scores.
411 Acts as a wrapper around the ME calculation from ``scores`` ([scoresa]_, [scoresb]_).
413 Parameters
414 ----------
415 cubes: iris.cube.CubeList
416 A CubeList containing exactly two cubes: a base and an "other" model,
417 this can be an analysis and the model.
418 preserved_coordinates: list[str] | str | None, default is None.
419 The coordinates that you wish to preserve in the calculaiton of the
420 ME. For example if you want a map of each time you can preserve
421 ["time","grid_latitude", "grid_longitude"] or if you want a time series
422 you can preserve ["time"], if you want to collapse to a single value
423 use `None`. The default is `None`.
425 Returns
426 -------
427 scores_cube: iris.cube.Cube
428 A cube containing the ME between the base and other cube.
429 """
430 base, other = _sort_cubes_for_verification(cubes)
432 # Copy the coordinates of the input cubes.
433 other_xr = xr.DataArray.from_iris(other)
434 base_xr = xr.DataArray.from_iris(base)
435 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
437 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
438 # apply scores, and then transform it back.
439 scores_cube = xr.DataArray.to_iris(
440 scores.continuous.additive_bias(
441 other_xr,
442 base_xr,
443 preserve_dims=preserve_dims,
444 )
445 )
447 # If time is aggregated out, attach a scalar time coordinate with bounds
448 # so plotting can display the aggregated period in the title.
449 try:
450 if not scores_cube.coords("time"): 450 ↛ 474line 450 didn't jump to line 474 because the condition on line 450 was always true
451 base_time = base.coord("time")
452 time_vals = (
453 base_time.bounds.flatten()
454 if base_time.has_bounds()
455 else base_time.points
456 )
457 t_start = float(time_vals[0])
458 t_end = float(time_vals[-1])
459 t_mid = 0.5 * (t_start + t_end)
461 scores_cube.add_aux_coord(
462 iris.coords.AuxCoord(
463 t_mid,
464 standard_name=base_time.standard_name,
465 long_name=base_time.long_name,
466 var_name=base_time.var_name,
467 units=base_time.units,
468 bounds=np.array([t_start, t_end]),
469 attributes=base_time.attributes.copy(),
470 )
471 )
472 except iris.exceptions.CoordinateNotFoundError:
473 pass
474 scores_cube.rename(f"Additive_Bias_of_{base.name()}")
475 return scores_cube
478def scores_correlation_pearsonr(
479 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
480):
481 r"""Calculate the Pearson's Correlation (PC) coefficient using scores.
483 Acts as a wrapper around the PC calculation from ``scores`` ([scoresa]_, [scoresb]_).
485 Parameters
486 ----------
487 cubes: iris.cube.CubeList
488 A CubeList containing exactly two cubes: a base and an "other" model,
489 this can be an analysis and the model.
490 preserved_coordinates: list[str] | str | None, default is None.
491 The coordinates that you wish to preserve in the calculation of the
492 PC. For example if you want a map of each time you can preserve
493 ["time","grid_latitude", "grid_longitude"] or if you want a time series
494 you can preserve ["time"], if you want to collapse to a single value
495 use `None`. The default is `None`.
497 Returns
498 -------
499 scores_cube: iris.cube.Cube
500 A cube containing the PC between the base and other cube.
501 """
502 base, other = _sort_cubes_for_verification(cubes)
504 # Copy the coordinates of the input cubes.
505 other_xr = xr.DataArray.from_iris(other)
506 base_xr = xr.DataArray.from_iris(base)
507 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
509 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
510 # apply scores, and then transform it back.
511 scores_cube = xr.DataArray.to_iris(
512 scores.continuous.correlation.pearsonr(
513 other_xr,
514 base_xr,
515 preserve_dims=preserve_dims,
516 )
517 )
519 # If time is aggregated out, attach a scalar time coordinate with bounds
520 # so plotting can display the aggregated period in the title.
521 try:
522 if not scores_cube.coords("time"): 522 ↛ 547line 522 didn't jump to line 547 because the condition on line 522 was always true
523 base_time = base.coord("time")
524 time_vals = (
525 base_time.bounds.flatten()
526 if base_time.has_bounds()
527 else base_time.points
528 )
529 t_start = float(time_vals[0])
530 t_end = float(time_vals[-1])
531 t_mid = 0.5 * (t_start + t_end)
533 scores_cube.add_aux_coord(
534 iris.coords.AuxCoord(
535 t_mid,
536 standard_name=base_time.standard_name,
537 long_name=base_time.long_name,
538 var_name=base_time.var_name,
539 units=base_time.units,
540 bounds=np.array([t_start, t_end]),
541 attributes=base_time.attributes.copy(),
542 )
543 )
544 except iris.exceptions.CoordinateNotFoundError:
545 pass
547 scores_cube.rename(f"Pearson_Correlation_of_{base.name()}")
548 return scores_cube
551def scores_crps_for_ensemble(
552 cubes: Cube | CubeList, method: str = "ecdf", control_member: int = 0
553) -> iris.Constraint:
554 r"""Calculate the CRPS for an ensemble.
556 Acts as a wrapper around the crps_for_ensemble from ``scores`` ([scoresa]_, [scoresb]_).
558 Lower CRPS values are better (implies experiment distribution is closer to control distribution/observations),
559 larger values are worse (implies distributions are dissimilar).
560 It is applicable across time and spatial scales as the focus is on the distribution of the values.
561 Default method is ecdf. ecdf is exact value from the empirical distributions,
562 whereas fair produces an approximated value based on a random sample of the underlying distribution.
564 See [CRPS]_ for further information.
566 Parameters
567 ----------
568 cubes: iris.cube.Cube
569 A Cube containing ensembles data
571 Returns
572 -------
573 crps: iris.cube.Cube
574 A cube containing the crps between the ensemble members and the control
575 """
576 if control_member != 0:
577 logger.warning("control member is usual 0")
579 if control_member not in cubes.coords("realization")[0].points:
580 new_control_member = cubes.coords("realization")[0].points[0]
581 logger.warning(
582 f"control member value {control_member} out of bounds, defaulting to control member={new_control_member}"
583 )
584 control_member = new_control_member
586 if cubes.coord("time").shape[0] == 1:
587 raise ValueError("Cube has only one time point.")
589 if cubes.coord("realization").shape[0] < 3:
590 raise ValueError("Cube should have one control member and at least two members")
592 ctrl = cubes.extract(generate_realization_constraint([control_member]))
593 ens_mem = cubes.extract(
594 generate_remove_single_ensemble_member_constraint(control_member)
595 )
597 # Realising the data in advance provides a large speedup
598 _ = ctrl.data
599 _ = ens_mem.data
600 del _
602 ctrl = xr.DataArray.from_iris(ctrl)
603 ens_mem = xr.DataArray.from_iris(ens_mem)
605 crps = xr.DataArray.to_iris(
606 scores.probability.crps_for_ensemble(
607 ens_mem,
608 ctrl,
609 ensemble_member_dim="realization",
610 method=method,
611 preserve_dims="time",
612 )
613 )
615 crps.rename(f"CRPS_of_{cubes[0].name()}")
616 _realization_callback(crps)
617 return crps
620def scores_pod_model_obs(
621 cubes: CubeList,
622 preserved_coordinates: list[str] | str | None,
623 threshold: str,
624 op_func: str,
625):
626 r"""
627 Compute the Probability of Detection (POD) score using Scores ([scoresa]_ [scoresb]_).
629 Parameters
630 ----------
631 cubes: iris.cube.CubeList
632 An iris cubelist containing model(s) and an observation cube.
633 preserved_coordinates: list | str | None
634 An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
635 then preserving coordinate 'time' will produce a probability of detection score for each timeslice (shape time). If None,
636 then it will return a single value score for all times/point locations.
637 threshold: str
638 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).
639 op_func: str
640 A string either containing 'lt' for less than or 'gt for greater than, to determine how the threshold is applied to the data
641 to generate the mask.
643 Returns
644 -------
645 cube: iris.cube
646 An iris cube, containing the probability of detection score for further plotting.
648 Notes
649 -----
650 The probability of detection calculates the proportion of observed events that meet a threshold that were correctly forecast by the model.
651 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
652 295k, that would be a positive hit. It does not take into account how far above/below a threshold a model forecasts.
654 It is calculated as .. math:: POD = \frac{true positives}{true positives + false negatives}
656 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.
658 POD produces a range of 0 to 1, where 1 is a perfect score.
659 """
660 # Split out model(s) and obs
661 models = CubeList()
662 for c in cubes:
663 if "observed" in c.long_name:
664 observed = c
665 else:
666 models.append(c)
668 # Setup cubelist to store results
669 scores_results = iris.cube.CubeList()
671 # Setup operators greater than, less than.
672 ops = {
673 "gt": operator.gt,
674 "lt": operator.lt,
675 }
677 try:
678 op = ops[op_func]
679 except KeyError as err:
680 raise ValueError(f"Operator {op_func} not supported.") from err
682 for model in models:
683 # Convert obs cubes to xarray and resolve preserved dimensions.
684 other_xr = xr.DataArray.from_iris(model)
685 base_xr = xr.DataArray.from_iris(observed)
686 preserve_dims = _resolve_preserve_dims(
687 observed, other_xr, preserved_coordinates
688 )
690 # Create event operator object using threshold and operator direction.
691 event_operator = scores.categorical.ThresholdEventOperator(
692 default_event_threshold=float(threshold), default_op_fn=op
693 )
695 # Generate binary fields using the event operator.
696 forecast_binary, observed_binary = event_operator.make_event_tables(
697 other_xr, base_xr
698 )
700 # Create binary contigency manager, as per Scores API, using transform to preserve preserve_dims
701 contingency_manager = scores.categorical.BinaryContingencyManager(
702 forecast_binary, observed_binary
703 ).transform(preserve_dims=preserve_dims)
705 # Get POD from the contigency manager, and convert back to an iris cube.
706 scores_cube = xr.DataArray.to_iris(
707 contingency_manager.probability_of_detection()
708 )
710 # Rename cube so it plots correctly alongside correcting cube units.
711 scores_cube.rename(
712 f"Probability_Of_Detection_{op_func}_{threshold}_{observed.name()}"
713 )
714 scores_cube.units = "1"
715 scores_cube.attributes["model_name"] = model.attributes["model_name"]
717 scores_results.append(scores_cube)
719 return scores_results