Coverage for src/CSET/operators/scoreswrappers.py: 89%
223 statements
« prev ^ index » next coverage.py v7.17.2, created at 2026-11-11 15:11 +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.coords
22import iris.exceptions
23import numpy as np
24import scores
25import scores.categorical
26import scores.continuous
27import scores.probability
28import xarray as xr
29from iris.cube import Cube, CubeList
30from iris.util import reverse
32from CSET._common import is_increasing
33from CSET.operators._utils import (
34 ensure_common_time_points,
35 fully_equalise_attributes,
36 get_cube_yxcoordname,
37)
38from CSET.operators.constraints import (
39 generate_realization_constraint,
40 generate_remove_single_ensemble_member_constraint,
41)
42from CSET.operators.read import _realization_callback
43from CSET.operators.regrid import regrid_onto_cube
45logger = logging.getLogger(__name__)
48def scores_rmse(
49 cubes: CubeList,
50 preserved_coordinates: list[str] | str | None = None,
51) -> CubeList:
52 r"""Calculate the Root Mean Square Error (RMSE) using scores.
54 Acts as a wrapper around the RMSE calculation from ``scores`` ([scoresa]_, [scoresb]_).
55 It is calculated as
57 .. math:: RMSE = \sqrt{\frac{1}{N} \Sigma(forecast - observations)^2}
59 Parameters
60 ----------
61 cubes: iris.cube.CubeList
62 A CubeList containing exactly two cubes: a base and an "other" model,
63 this can be an analysis and the model.
64 preserved_coordinates: list[str] | str | None, default is None.
65 The coordinates (or xarray dimension names) that you wish to preserve in the calculaiton of the
66 RMSE. For example if you want a map of each time you can preserve
67 ["time","grid_latitude", "grid_longitude"] or if you want a time series
68 you can preserve ["time"], if you want to collapse to a single value
69 use `None`. The default is `None`.
71 Returns
72 -------
73 scores_cubelist: iris.cube.CubeList
74 A cubelist containing the RMSE between the base and other cube.
75 """
76 scores_cubelist = CubeList()
77 base, others = _split_base_and_other(cubes)
79 for other in others:
80 scores_cube = _make_scores_cube(base, other, "rmse", preserved_coordinates)
82 scores_cubelist.append(scores_cube)
84 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
87def scores_mae(
88 cubes: CubeList,
89 preserved_coordinates: list[str] | str | None = None,
90) -> CubeList:
91 r"""Calculate the Mean Absolute Error (MAE) using scores.
93 Acts as a wrapper around the MAE calculation from ``scores`` ([scoresa]_, [scoresb]_).
95 Parameters
96 ----------
97 cubes: iris.cube.CubeList
98 A CubeList containing exactly two cubes: a base and an "other" model,
99 this can be an analysis and the model.
100 preserved_coordinates: list[str] | str | None, default is None.
101 The coordinates that you wish to preserve in the calculaiton of the
102 MAE. For example if you want a map of each time you can preserve
103 ["time","grid_latitude", "grid_longitude"] or if you want a time series
104 you can preserve ["time"], if you want to collapse to a single value
105 use `None`. The default is `None`.
107 Returns
108 -------
109 scores_cubelist: iris.cube.CubeList
110 A cubelist containing the MAE between the base and other cube(s).
111 """
112 scores_cubelist = CubeList()
113 base, others = _split_base_and_other(cubes)
115 for other in others:
116 scores_cube = _make_scores_cube(base, other, "mae", preserved_coordinates)
117 scores_cubelist.append(scores_cube)
119 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
122def scores_additive_bias(
123 cubes: CubeList,
124 preserved_coordinates: list[str] | str | None = None,
125) -> CubeList:
126 r"""Calculate the Additive Bias (Mean Error) using scores.
128 Acts as a wrapper around the ME calculation from ``scores`` ([scoresa]_, [scoresb]_).
130 Parameters
131 ----------
132 cubes: iris.cube.CubeList
133 A CubeList containing exactly two cubes: a base and an "other" model,
134 this can be an analysis and the model.
135 preserved_coordinates: list[str] | str | None, default is None.
136 The coordinates that you wish to preserve in the calculaiton of the
137 ME. For example if you want a map of each time you can preserve
138 ["time","grid_latitude", "grid_longitude"] or if you want a time series
139 you can preserve ["time"], if you want to collapse to a single value
140 use `None`. The default is `None`.
142 Returns
143 -------
144 scores_cubelist: iris.cube.CubeList
145 A cubelist containing the ME between the base and other cube(s).
146 """
147 scores_cubelist = CubeList()
148 base, others = _split_base_and_other(cubes)
150 for other in others:
151 scores_cube = _make_scores_cube(
152 base, other, "additive_bias", preserved_coordinates
153 )
155 scores_cubelist.append(scores_cube)
157 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
160def scores_correlation_pearsonr(
161 cubes: CubeList,
162 preserved_coordinates: list[str] | str | None = None,
163) -> CubeList:
164 r"""Calculate the Pearson's Correlation (PC) coefficient using scores.
166 Acts as a wrapper around the PC calculation from ``scores`` ([scoresa]_, [scoresb]_).
168 Parameters
169 ----------
170 cubes: iris.cube.CubeList
171 A CubeList containing exactly two cubes: a base and an "other" model,
172 this can be an analysis and the model.
173 preserved_coordinates: list[str] | str | None, default is None.
174 The coordinates that you wish to preserve in the calculation of the
175 PC. For example if you want a map of each time you can preserve
176 ["time","grid_latitude", "grid_longitude"] or if you want a time series
177 you can preserve ["time"], if you want to collapse to a single value
178 use `None`. The default is `None`.
180 Returns
181 -------
182 scores_cubelist: iris.cube.CubeList
183 A cubelist containing the PC between the base and other cube(s).
184 """
185 scores_cubelist = CubeList()
186 base, others = _split_base_and_other(cubes)
188 for other in others:
189 scores_cube = _make_scores_cube(
190 base, other, "pearson_correlation", preserved_coordinates
191 )
193 scores_cubelist.append(scores_cube)
195 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
198def scores_crps_for_ensemble(
199 cubes: Cube | CubeList, method: str = "ecdf", control_member: int = 0
200) -> Cube:
201 r"""Calculate the CRPS for an ensemble.
203 Acts as a wrapper around the crps_for_ensemble from ``scores`` ([scoresa]_, [scoresb]_).
205 Lower CRPS values are better (implies experiment distribution is closer to control distribution/observations),
206 larger values are worse (implies distributions are dissimilar).
207 It is applicable across time and spatial scales as the focus is on the distribution of the values.
208 Default method is ecdf. ecdf is exact value from the empirical distributions,
209 whereas fair produces an approximated value based on a random sample of the underlying distribution.
211 See [CRPS]_ for further information.
213 Parameters
214 ----------
215 cubes: iris.cube.Cube
216 A Cube containing ensembles data
218 method: str ["ecfd" or "fair"]
219 Determines the method to use for calculating the CRPS. Defaults to "ecdf".
221 control_member: int
222 What the realisation the control member of the ensemble is. Defaults to 0. Sometimes this is 1.
224 Returns
225 -------
226 crps: iris.cube.Cube
227 A cube containing the crps between the ensemble members and the control
228 """
229 if control_member != 0:
230 logger.warning("control member is usual 0")
232 if control_member not in cubes.coords("realization")[0].points:
233 new_control_member = cubes.coords("realization")[0].points[0]
234 logger.warning(
235 f"control member value {control_member} out of bounds, defaulting to control member={new_control_member}"
236 )
237 control_member = new_control_member
239 if cubes.coord("time").shape[0] == 1:
240 raise ValueError("Cube has only one time point.")
242 if cubes.coord("realization").shape[0] < 3:
243 raise ValueError("Cube should have one control member and at least two members")
245 ctrl = cubes.extract(generate_realization_constraint([control_member]))
246 ens_mem = cubes.extract(
247 generate_remove_single_ensemble_member_constraint(control_member)
248 )
250 # Realising the data in advance provides a large speedup
251 _ = ctrl.data
252 _ = ens_mem.data
253 del _
255 ctrl = xr.DataArray.from_iris(ctrl)
256 ens_mem = xr.DataArray.from_iris(ens_mem)
258 crps = xr.DataArray.to_iris(
259 scores.probability.crps_for_ensemble(
260 ens_mem,
261 ctrl,
262 ensemble_member_dim="realization",
263 method=method,
264 preserve_dims="time",
265 )
266 )
268 crps.rename(f"CRPS_of_{cubes[0].name()}")
269 _realization_callback(crps)
270 return crps
273def _scores_categorical_metric(
274 cubes: CubeList,
275 preserved_coordinates: list[str] | str | None,
276 threshold: str,
277 op_func: str,
278 metric: str,
279) -> CubeList:
280 """
281 Prepare cubes for computing categorical metrics using Scores.
283 Parameters
284 ----------
285 cubes: iris.cube.CubeList
286 An iris cubelist containing model(s) and an observation cube.
287 preserved_coordinates: list | str | None
288 An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
289 then preserving coordinate 'time' will produce a probability of detection score for each timeslice (shape time). If None,
290 then it will return a single value score for all times/point locations.
291 threshold: str
292 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).
293 op_func: str
294 A string either containing 'lt' for less than or 'gt for greater than, to determine how the threshold is applied to the data
295 to generate the mask.
296 metric: str
297 The scores metric to compute.
299 Returns
300 -------
301 scores_results: iris.cube.CubeList
302 An iris cubelist, containing the scores metric for each model for further plotting.
304 """
305 # Split obs/models
306 models = CubeList()
307 for c in cubes:
308 if "observed" in c.long_name:
309 observed = c
310 else:
311 models.append(c)
313 ops = {
314 "gt": operator.gt,
315 "lt": operator.lt,
316 }
318 # Check if this exists.
319 try:
320 op = ops[op_func]
321 except KeyError as err:
322 raise ValueError(f"Operator {op_func} not supported.") from err
324 scores_results = CubeList()
326 for model in models:
327 # Convert obs cubes to xarray and resolve preserved dimensions.
328 other_xr = xr.DataArray.from_iris(model)
329 base_xr = xr.DataArray.from_iris(observed)
330 preserve_dims = _resolve_preserve_dims(
331 observed, other_xr, preserved_coordinates
332 )
334 # Create event operator object using threshold and operator direction.
335 event_operator = scores.categorical.ThresholdEventOperator(
336 default_event_threshold=float(threshold), default_op_fn=op
337 )
339 # Generate binary fields using the event operator.
340 forecast_binary, observed_binary = event_operator.make_event_tables(
341 other_xr, base_xr
342 )
344 # Create binary contigency manager, as per Scores API, using transform to preserve preserve_dims
345 contingency_manager = scores.categorical.BinaryContingencyManager(
346 forecast_binary, observed_binary
347 ).transform(preserve_dims=preserve_dims)
349 # Compute required categorical score.
350 if metric == "pod":
351 result = contingency_manager.probability_of_detection()
352 name = "Probability_Of_Detection"
354 elif metric == "ets":
355 result = contingency_manager.equitable_threat_score()
356 name = "Equitable_Threat_Score"
358 elif metric == "fb":
359 result = contingency_manager.frequency_bias()
360 name = "Frequency_Bias"
362 elif metric == "pfd": 362 ↛ 367line 362 didn't jump to line 367 because the condition on line 362 was always true
363 result = contingency_manager.probability_of_false_detection()
364 name = "Probability_Of_False_Detection"
366 else:
367 raise ValueError(f"Unknown metric {metric}")
369 scores_cube = xr.DataArray.to_iris(result)
370 _fix_spatial_coord(cubes[0], scores_cube)
371 scores_cube.rename(f"{name}_{op_func}_{threshold}_{observed.name()}")
372 scores_cube.units = "1"
373 scores_cube.attributes["model_name"] = model.attributes["model_name"]
375 scores_results.append(scores_cube)
377 return scores_results
380def scores_pod(
381 cubes,
382 preserved_coordinates,
383 threshold,
384 op_func,
385) -> CubeList:
386 r"""
387 Compute the Probability of Detection (POD) score using Scores ([scoresa]_ [scoresb]_).
389 Parameters
390 ----------
391 cubes: iris.cube.CubeList
392 An iris cubelist containing model(s) and an observation cube.
393 preserved_coordinates: list | str | None
394 An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
395 then preserving coordinate 'time' will produce the equitable threat score for each timeslice (shape time). If None,
396 then it will return a single value score for all times/point locations.
397 threshold: str
398 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).
399 op_func: str
400 A string either containing 'lt' for less than or 'gt for greater than, to determine how the threshold is applied to the data
401 to generate the mask.
403 Returns
404 -------
405 iris.cube.CubeList
406 An iris cubelist, containing the probability of detection score for each model for further plotting.
409 Notes
410 -----
411 The probability of detection calculates the proportion of observed events that meet a threshold that were correctly forecast by the model.
412 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
413 295k, that would be a positive hit. It does not take into account how far above/below a threshold a model forecasts.
415 It is calculated as .. math:: POD = \frac{true positives}{true positives + false negatives}
417 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.
419 POD produces a range of 0 to 1, where 1 is a perfect score.
420 """
421 return _scores_categorical_metric(
422 cubes,
423 preserved_coordinates,
424 threshold,
425 op_func,
426 "pod",
427 )
430def scores_ets(
431 cubes: CubeList,
432 preserved_coordinates: list[str] | str | None,
433 threshold: str,
434 op_func: str,
435) -> CubeList:
436 r"""
437 Compute the Equitable Threat Score (ETS) score using Scores ([scoresa]_ [scoresb]_).
439 Parameters
440 ----------
441 cubes: iris.cube.CubeList
442 An iris cubelist containing model(s) and an observation cube.
443 preserved_coordinates: list | str | None
444 An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
445 then preserving coordinate 'time' will produce the equitable threat score for each timeslice (shape time). If None,
446 then it will return a single value score for all times/point locations.
447 threshold: str
448 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).
449 op_func: str
450 A string either containing 'lt' for less than or 'gt for greater than, to determine how the threshold is applied to the data
451 to generate the mask.
453 Returns
454 -------
455 iris.cube.CubeList
456 An iris cubelist, containing the probability of detection score for each model for further plotting.
458 Notes
459 -----
460 The Equitable Threat Score (ETS) evaluates the accuracy of forecasts for events that meet a specified threshold,
461 hile accounting for correct forecasts that could occur purely by chance. Unlike the Probability of Detection (POD),
462 ETS considers hits, misses, and false alarms, providing a more balanced assessment of forecast skill.
464 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
465 would be counted as a hit. ETS adjusts the total number of hits by removing the number of hits expected due to random chance.
467 It is calculated as:
469 .. math::
471 ETS = \frac{hits - hits_{random}}
472 {hits + misses + false\ alarms - hits_{random}}
474 where
476 hits_{random} = \frac{(hits + misses)(hits + false\ alarms)}{total count}
478 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
479 random chance.
480 """
481 return _scores_categorical_metric(
482 cubes,
483 preserved_coordinates,
484 threshold,
485 op_func,
486 "ets",
487 )
490def scores_pfd(
491 cubes: CubeList,
492 preserved_coordinates: list[str] | str | None,
493 threshold: str,
494 op_func: str,
495) -> CubeList:
496 r"""
497 Compute the Probability of False Detection (PFD) score using Scores ([scoresa]_ [scoresb]_).
499 Parameters
500 ----------
501 cubes: iris.cube.CubeList
502 An iris cubelist containing model(s) and an observation cube.
503 preserved_coordinates: list | str | None
504 An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
505 then preserving coordinate 'time' will produce the probability of false detection score for each timeslice (shape time). If None,
506 then it will return a single value score for all times/point locations.
507 threshold: str
508 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).
509 op_func: str
510 A string either containing 'lt' for less than or 'gt' for greater than, to determine how the threshold is applied to the data
511 to generate the mask.
513 Returns
514 -------
515 iris.cube.CubeList
516 An iris cubelist, containing the probability of false detection score for each model for further plotting.
518 Notes
519 -----
520 The Probability of False Detection (PFD) measures the proportion of observed non-events
521 that were incorrectly forecast as events. It is a measure of the false alarm rate and
522 provides information on how often a forecast system predicts threshold exceedances when
523 none actually occurred.
525 For example, if the threshold is 290 K and op_func is gt (greater than), an observation
526 of 288 K and a forecast of 295 K would be considered a false alarm, since the model
527 predicted an event but the observed value did not exceed the threshold.
529 It is calculated as:
531 .. math::
533 POFD = \frac{false\ alarms}
534 {false\ alarms + true\ negatives}
536 where
538 false alarms
539 Number of occasions where an event was forecast but did not occur.
541 true negatives
542 Number of occasions where neither the forecast nor the observations
543 indicated an event.
545 PFD ranges from 0 to 1, where 0 indicates a perfect score with no false alarms,
546 and 1 indicates that every observed non-event was incorrectly forecast as an event.
548 Lower values are therefore better.
549 """
550 return _scores_categorical_metric(
551 cubes,
552 preserved_coordinates,
553 threshold,
554 op_func,
555 "pfd",
556 )
559def scores_frequency_bias(
560 cubes: CubeList,
561 preserved_coordinates: list[str] | str | None,
562 threshold: str,
563 op_func: str,
564) -> CubeList:
565 r"""
566 Compute the Frequency Bias (FB) score using Scores ([scoresa]_ [scoresb]_).
568 Parameters
569 ----------
570 cubes: iris.cube.CubeList
571 An iris cubelist containing model(s) and an observation cube.
572 preserved_coordinates: list | str | None
573 An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
574 then preserving coordinate 'time' will produce the frequency bias score for each timeslice (shape time). If None,
575 then it will return a single value score for all times/point locations.
576 threshold: str
577 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).
578 op_func: str
579 A string either containing 'lt' for less than or 'gt' for greater than, to determine how the threshold is applied to the data
580 to generate the mask.
582 Returns
583 -------
584 iris.cube.CubeList
585 An iris cube, containing the frequency bias score for each model for further plotting.
587 Notes
588 -----
589 The Frequency Bias (FB) measures whether a forecasting system predicts an
590 event too frequently or too infrequently compared to observations. Unlike
591 metrics such as the Probability of Detection (POD) or Equitable Threat Score (ETS),
592 Frequency Bias does not assess the accuracy of forecast locations or timings,
593 only the overall frequency with which events are forecast.
595 For example, if the threshold is 290 K and op_func is gt (greater than), and
596 events exceeding this threshold are forecast twice as often as they are observed,
597 the frequency bias would be approximately 2. Conversely, if events are forecast
598 only half as often as they occur, the frequency bias would be approximately 0.5.
600 It is calculated as:
602 .. math::
604 Frequency\ Bias = \frac{hits + false\ alarms}
605 {hits + misses}
607 where
609 hits
610 Number of occasions where an event was forecast and observed.
612 false alarms
613 Number of occasions where an event was forecast but not observed.
615 misses
616 Number of occasions where an event was observed but not forecast.
618 Frequency Bias ranges from 0 to infinity, where:
620 * 1 indicates the forecast predicts events at the correct frequency.
621 * Greater than 1 indicates overforecasting of events.
622 * Less than 1 indicates underforecasting of events.
624 A perfect frequency bias score is therefore 1, although a value of 1 does not
625 necessarily imply a skillful forecast, as hits and false alarms may compensate
626 for one another.
627 """
628 return _scores_categorical_metric(
629 cubes,
630 preserved_coordinates,
631 threshold,
632 op_func,
633 "fb",
634 )
637def _make_scores_cube(
638 base: Cube, other: Cube, metric: str, preserved_coordinates: list[str]
639) -> Cube:
640 r"""Make the scores cube using the given scores metric.
642 Parameters
643 ----------
644 base: iris.cube.Cube
645 The cube from the "analysis" or observed cube.
646 other: iris.cube.Cube
647 The cube from the model.
649 metric: str
650 The scores metric to compute.
652 preserved_coordinates: list[str]
653 The list of preserved coordinates given by a user.
655 Returns
656 -------
657 scores_cube: iris.cube.Cube
658 The cube containing the calculated scores metric.
661 """
662 if not base.long_name or not "observed" in base.long_name:
663 base, other = _process_cubes_for_verification(base, other)
665 other_xr = xr.DataArray.from_iris(other)
666 base_xr = xr.DataArray.from_iris(base)
667 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
669 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
670 # apply scores, and then transform it back.
672 if metric == "rmse":
673 scores_cube = xr.DataArray.to_iris(
674 scores.continuous.rmse(other_xr, base_xr, preserve_dims=preserve_dims)
675 )
676 scores_cube.rename(f"RMSE_of_{base.name()}")
677 elif metric == "mae":
678 scores_cube = xr.DataArray.to_iris(
679 scores.continuous.mae(other_xr, base_xr, preserve_dims=preserve_dims)
680 )
681 scores_cube.rename(f"MAE_of_{base.name()}")
682 elif metric == "additive_bias":
683 scores_cube = xr.DataArray.to_iris(
684 scores.continuous.additive_bias(
685 other_xr, base_xr, preserve_dims=preserve_dims
686 )
687 )
688 scores_cube.rename(f"Additive_Bias_of_{base.name()}")
689 elif metric == "pearson_correlation": 689 ↛ 697line 689 didn't jump to line 697 because the condition on line 689 was always true
690 scores_cube = xr.DataArray.to_iris(
691 scores.continuous.correlation.pearsonr(
692 other_xr, base_xr, preserve_dims=preserve_dims
693 )
694 )
695 scores_cube.rename(f"Pearson_Correlation_of_{base.name()}")
696 else:
697 raise ValueError(f"Scores Unknown metric: {metric}")
699 _fix_spatial_coord(base, scores_cube)
701 _attach_scaler_time_coord_maybe(scores_cube, base)
702 model_name = other.attributes["model_name"]
703 scores_cube.attributes["model_name"] = model_name
704 return scores_cube
707def _fix_spatial_coord(base_cube: Cube, scores_cube: Cube):
708 if not scores_cube.coords("grid_latitude") or not base_cube.coords("grid_latitude"):
709 return
710 if not scores_cube.coord("grid_latitude").coord_system: 710 ↛ exitline 710 didn't return from function '_fix_spatial_coord' because the condition on line 710 was always true
711 for xy_name in ("grid_latitude", "grid_longitude"):
712 scores_cube.coord(xy_name).coord_system = base_cube.coord(
713 xy_name
714 ).coord_system
717def _sort_cube_into_base_and_other(cubes: CubeList) -> tuple[Cube, CubeList]:
718 r"""Sorts cube into base and other models.
720 Parameters
721 ----------
722 cubes: iris.cube.CubeList
723 A CubeList of multiple cubes. One base cube and other model cubes.
725 Returns
726 -------
727 base: iris.cube.Cube
728 The cube from the "analysis" or observed cube in the same format as the other model.
729 others: iris.cube.CubeList
730 The cube list of containing the cube(s) from the model in the same format as the base model.
732 """
733 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
734 others: CubeList = cubes.extract(
735 iris.Constraint(
736 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
737 )
738 )
740 return base, others
743def _split_base_and_other(cubes: CubeList):
744 r"""Split the cube into base and other cubes.
746 Split depends on whether there
747 is an observed cube in the cubes. If there is an observed cube,
748 then 'base' is the observed cube, if not then 'base' is the comparison
749 model cube.
751 Parameters
752 ----------
753 cubes: iris.cube.CubeList
754 Cubes to split into base and other cubes.
756 Returns
757 -------
758 tuple
759 A tuple containing a base cube, and other cube/cubelist.
761 """
762 obs_cube = [cb for cb in cubes if "observed" in (cb.long_name or "")]
763 if obs_cube:
764 if len(obs_cube) > 1: 764 ↛ 765line 764 didn't jump to line 765 because the condition on line 764 was never true
765 raise ValueError(
766 f"Expected exactly one 'observed' cube, found {len(obs_cube)}"
767 )
768 base = obs_cube[0]
769 others = CubeList(cb for cb in cubes if cb is not base)
770 return base, others
772 return _sort_cube_into_base_and_other(cubes)
775def _ensure_increasing_pressure_coordinates(cubes: CubeList) -> CubeList:
776 r"""Ensure the pressure coordinate is increasing.
778 Parameters
779 ----------
780 cubes: iris.cube.CubeList
781 A CubeList of n cubes
783 Returns
784 -------
785 Cubes: iris.cube.CubeList
786 The original cube list but where each cube is ensured to have an increasing pressure coordinate.
787 """
788 for cube in cubes:
789 try:
790 if len(cube.coord("pressure").points) > 2 and not is_increasing(
791 cube.coord("pressure").points
792 ):
793 reverse(cube, "pressure")
795 except iris.exceptions.CoordinateNotFoundError:
796 pass
799def _process_cubes_for_verification(base: Cube, other: Cube) -> tuple[Cube, Cube]:
800 r"""Prepare cubes ready for verification in scores.
802 Parameters
803 ----------
804 base: iris.cube.Cube
805 The cube from the "analysis" or observed cube.
806 other: iris.cube.Cube
807 The cube from the model.
809 Returns
810 -------
811 base: iris.cube.Cube
812 The cube from the "analysis" or observed cube in the same format as the other model.
813 other: iris.cube.Cube
814 The cube from the model in the same format as the base model.
816 Notes
817 -----
818 This operator is used for sorting the data into the correct format. It
819 is likely going to need to be refactored out of CSET and perhaps moved into
820 `CSET._utils` given common code between here and `misc.difference`.
821 """
822 # raises value error if cube time points don't match
823 ensure_common_time_points(base, other)
825 # Set cubes into correct format using code from difference operator
827 other_model_name = other.attributes["model_name"]
829 if base.coords("pressure") and other.coords("pressure"):
830 common_pressure = np.intersect1d(
831 base.coord("pressure").points,
832 other.coord("pressure").points,
833 )
835 pressure_constraint = iris.Constraint(pressure=common_pressure.tolist())
837 base = base.extract(pressure_constraint)
838 other = other.extract(pressure_constraint)
840 # Get spatial coord names.
841 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
842 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
844 # Ensure cubes to compare are on common differencing grid.
845 # This is triggered if either
846 # i) latitude and longitude shapes are not the same. Note grid points
847 # are not compared directly as these can differ through rounding
848 # errors.
849 # ii) or variables are known to often sit on different grid staggering
850 # in different models (e.g. cell center vs cell edge), as is the case
851 # for UM and LFRic comparisons.
852 # In future greater choice of regridding method might be applied depending
853 # on variable type. Linear regridding can in general be appropriate for smooth
854 # variables. Care should be taken with interpretation of differences
855 # given this dependency on regridding.
857 if (
858 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
859 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
860 ) or (
861 base.long_name
862 in [
863 "eastward_wind_at_10m",
864 "northward_wind_at_10m",
865 "northward_wind_at_cell_centres",
866 "eastward_wind_at_cell_centres",
867 "zonal_wind_at_pressure_levels",
868 "meridional_wind_at_pressure_levels",
869 "potential_vorticity_at_pressure_levels",
870 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
871 ]
872 ):
873 logger.debug("Linear regridding base cube to other grid to compute differences")
874 base = regrid_onto_cube(base, other, method="Linear")
876 # Figure out if we are comparing between UM and LFRic; flip array if so.
877 base_lat_direction = is_increasing(base.coord(base_lat_name).points)
878 other_lat_direction = is_increasing(other.coord(other_lat_name).points)
879 if base_lat_direction != other_lat_direction: 879 ↛ 881line 879 didn't jump to line 881 because the condition on line 879 was never true
880 # Copy base cube for correct coordinate information.
881 other_tmp = base.copy()
882 # Flip the data and place in the copied cube.
883 other_tmp.data = np.flip(
884 other.data, other.coord(other_lat_name).cube_dims(other)
885 )
886 # Use original name and units from the other cube.
887 other_tmp.rename(other.name())
888 other_tmp.units = other.units
889 # Replace the cube.
890 other = other_tmp
892 # Equalise attributes so we can merge.
893 fully_equalise_attributes(CubeList([base, other]))
895 other.attributes["model_name"] = other_model_name
896 logger.debug("Base: %s\nOther: %s", base, other)
898 return base, other
901def _resolve_preserve_dims(
902 cube: Cube,
903 data_array: xr.DataArray,
904 preserved_coordinates: list[str] | str | None,
905) -> list[str] | None:
906 r"""Resolve preserve coordinates to xarray dimension names.
908 The ``scores`` package expects preserve dimensions to match xarray
909 dimension names. In Iris data, commonly used coordinates such as ``time``
910 may be auxiliary coordinates attached to a differently named dimension
911 (e.g. ``dim0``). This helper maps coordinate names to their underlying
912 dimension names and helps to convert from iris to xarray coordinate dimension names.
914 Parameters
915 ----------
916 cube: iris.cube.Cube
917 The ccomparison model cube.
918 data_array: xr.DataArray
919 The comparison model cube, but in xarray format.
920 preserved_coordinates: list[str]
921 The list of preserved coordinates given by a user.
924 Returns
925 -------
926 preserved_dims : list[str]
927 List of preserved dimension names in xarray convention.
929 """
930 if preserved_coordinates is None:
931 return None
933 coord_names = (
934 [preserved_coordinates]
935 if isinstance(preserved_coordinates, str)
936 else preserved_coordinates
937 )
938 preserve_dims: list[str] = []
940 for coord_name in coord_names:
941 # Already an xarray dimension name.
942 if coord_name in data_array.dims:
943 if coord_name not in preserve_dims: 943 ↛ 945line 943 didn't jump to line 945 because the condition on line 943 was always true
944 preserve_dims.append(coord_name)
945 continue
947 # Otherwise, map coordinate name to dimension index/indices.
948 try:
949 dim_indices = cube.coord_dims(coord_name)
950 except iris.exceptions.CoordinateNotFoundError:
951 # Keep original name so scores raises a clear error for unknown keys.
952 if coord_name not in preserve_dims:
953 preserve_dims.append(coord_name)
954 continue
956 for dim_index in dim_indices:
957 dim_name = data_array.dims[dim_index]
958 if dim_name not in preserve_dims:
959 preserve_dims.append(dim_name)
961 return preserve_dims
964def _attach_scaler_time_coord_maybe(scores_cube: Cube, base: Cube) -> None:
965 r"""Attaches scaler time coordinate if time is aggregated out.
967 In place function that attaches a scaler time coordinate to scores_cube
968 if time is aggregated out so plotting can display the aggregated period in the title.
970 Parameters
971 ----------
972 scores_cube: iris.cube.Cube
973 The calculated scores cube.
974 base: iris.cube.Cube
975 The cube from the "analysis" or observed cube.
977 Returns
978 -------
979 None
981 """
982 try:
983 if not scores_cube.coords("time") and not scores_cube.coords("forecast_period"):
984 base_time = base.coord("time")
985 time_vals = (
986 base_time.bounds.flatten()
987 if base_time.has_bounds()
988 else base_time.points
989 )
990 t_start = float(time_vals[0])
991 t_end = float(time_vals[-1])
992 t_mid = 0.5 * (t_start + t_end)
994 scores_cube.add_aux_coord(
995 iris.coords.AuxCoord(
996 t_mid,
997 standard_name=base_time.standard_name,
998 long_name=base_time.long_name,
999 var_name=base_time.var_name,
1000 units=base_time.units,
1001 bounds=np.array([t_start, t_end]),
1002 attributes=base_time.attributes.copy(),
1003 )
1004 )
1006 except iris.exceptions.CoordinateNotFoundError:
1007 pass