Coverage for src/CSET/operators/scoreswrappers.py: 89%
226 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-19 12:33 +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 fully_equalise_attributes, get_cube_yxcoordname
34from CSET.operators.constraints import (
35 generate_realization_constraint,
36 generate_remove_single_ensemble_member_constraint,
37)
38from CSET.operators.misc import _extract_common_time_points
39from CSET.operators.read import _realization_callback
40from CSET.operators.regrid import regrid_onto_cube
42logger = logging.getLogger(__name__)
45def scores_rmse(
46 cubes: CubeList,
47 preserved_coordinates: list[str] | str | None = None,
48) -> CubeList:
49 r"""Calculate the Root Mean Square Error (RMSE) using scores.
51 Acts as a wrapper around the RMSE calculation from ``scores`` ([scoresa]_, [scoresb]_).
52 It is calculated as
54 .. math:: RMSE = \sqrt{\frac{1}{N} \Sigma(forecast - observations)^2}
56 Parameters
57 ----------
58 cubes: iris.cube.CubeList
59 A CubeList containing exactly two cubes: a base and an "other" model,
60 this can be an analysis and the model.
61 preserved_coordinates: list[str] | str | None, default is None.
62 The coordinates (or xarray dimension names) that you wish to preserve in the calculaiton of the
63 RMSE. For example if you want a map of each time you can preserve
64 ["time","grid_latitude", "grid_longitude"] or if you want a time series
65 you can preserve ["time"], if you want to collapse to a single value
66 use `None`. The default is `None`.
68 Returns
69 -------
70 scores_cubelist: iris.cube.CubeList
71 A cubelist containing the RMSE between the base and other cube.
72 """
73 scores_cubelist = CubeList()
75 base, others = _split_base_and_other(cubes)
77 for other in others:
78 base, other = _process_cubes_for_verification(base, other)
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 base, other = _process_cubes_for_verification(base, other)
117 scores_cube = _make_scores_cube(base, other, "mae", preserved_coordinates)
118 scores_cubelist.append(scores_cube)
120 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
123def scores_additive_bias(
124 cubes: CubeList,
125 preserved_coordinates: list[str] | str | None = None,
126) -> CubeList:
127 r"""Calculate the Additive Bias (Mean Error) using scores.
129 Acts as a wrapper around the ME calculation from ``scores`` ([scoresa]_, [scoresb]_).
131 Parameters
132 ----------
133 cubes: iris.cube.CubeList
134 A CubeList containing exactly two cubes: a base and an "other" model,
135 this can be an analysis and the model.
136 preserved_coordinates: list[str] | str | None, default is None.
137 The coordinates that you wish to preserve in the calculaiton of the
138 ME. For example if you want a map of each time you can preserve
139 ["time","grid_latitude", "grid_longitude"] or if you want a time series
140 you can preserve ["time"], if you want to collapse to a single value
141 use `None`. The default is `None`.
143 Returns
144 -------
145 scores_cubelist: iris.cube.CubeList
146 A cubelist containing the ME between the base and other cube(s).
147 """
148 scores_cubelist = CubeList()
149 base, others = _split_base_and_other(cubes)
151 for other in others:
152 base, other = _process_cubes_for_verification(base, other)
154 scores_cube = _make_scores_cube(
155 base, other, "additive_bias", preserved_coordinates
156 )
158 scores_cubelist.append(scores_cube)
160 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
163def scores_correlation_pearsonr(
164 cubes: CubeList,
165 preserved_coordinates: list[str] | str | None = None,
166) -> CubeList:
167 r"""Calculate the Pearson's Correlation (PC) coefficient using scores.
169 Acts as a wrapper around the PC calculation from ``scores`` ([scoresa]_, [scoresb]_).
171 Parameters
172 ----------
173 cubes: iris.cube.CubeList
174 A CubeList containing exactly two cubes: a base and an "other" model,
175 this can be an analysis and the model.
176 preserved_coordinates: list[str] | str | None, default is None.
177 The coordinates that you wish to preserve in the calculation of the
178 PC. For example if you want a map of each time you can preserve
179 ["time","grid_latitude", "grid_longitude"] or if you want a time series
180 you can preserve ["time"], if you want to collapse to a single value
181 use `None`. The default is `None`.
183 Returns
184 -------
185 scores_cubelist: iris.cube.CubeList
186 A cubelist containing the PC between the base and other cube(s).
187 """
188 scores_cubelist = CubeList()
189 base, others = _split_base_and_other(cubes)
191 for other in others:
192 base, other = _process_cubes_for_verification(base, other)
194 scores_cube = _make_scores_cube(
195 base, other, "pearson_correlation", preserved_coordinates
196 )
198 scores_cubelist.append(scores_cube)
200 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
203def scores_crps_for_ensemble(
204 cubes: Cube | CubeList, method: str = "ecdf", control_member: int = 0
205) -> Cube:
206 r"""Calculate the CRPS for an ensemble.
208 Acts as a wrapper around the crps_for_ensemble from ``scores`` ([scoresa]_, [scoresb]_).
210 Lower CRPS values are better (implies experiment distribution is closer to control distribution/observations),
211 larger values are worse (implies distributions are dissimilar).
212 It is applicable across time and spatial scales as the focus is on the distribution of the values.
213 Default method is ecdf. ecdf is exact value from the empirical distributions,
214 whereas fair produces an approximated value based on a random sample of the underlying distribution.
216 See [CRPS]_ for further information.
218 Parameters
219 ----------
220 cubes: iris.cube.Cube
221 A Cube containing ensembles data
223 method: str ["ecfd" or "fair"]
224 Determines the method to use for calculating the CRPS. Defaults to "ecdf".
226 control_member: int
227 What the realisation the control member of the ensemble is. Defaults to 0. Sometimes this is 1.
229 Returns
230 -------
231 crps: iris.cube.Cube
232 A cube containing the crps between the ensemble members and the control
233 """
234 if control_member != 0:
235 logger.warning("control member is usual 0")
237 if control_member not in cubes.coords("realization")[0].points:
238 new_control_member = cubes.coords("realization")[0].points[0]
239 logger.warning(
240 f"control member value {control_member} out of bounds, defaulting to control member={new_control_member}"
241 )
242 control_member = new_control_member
244 if cubes.coord("time").shape[0] == 1:
245 raise ValueError("Cube has only one time point.")
247 if cubes.coord("realization").shape[0] < 3:
248 raise ValueError("Cube should have one control member and at least two members")
250 ctrl = cubes.extract(generate_realization_constraint([control_member]))
251 ens_mem = cubes.extract(
252 generate_remove_single_ensemble_member_constraint(control_member)
253 )
255 # Realising the data in advance provides a large speedup
256 _ = ctrl.data
257 _ = ens_mem.data
258 del _
260 ctrl = xr.DataArray.from_iris(ctrl)
261 ens_mem = xr.DataArray.from_iris(ens_mem)
263 crps = xr.DataArray.to_iris(
264 scores.probability.crps_for_ensemble(
265 ens_mem,
266 ctrl,
267 ensemble_member_dim="realization",
268 method=method,
269 preserve_dims="time",
270 )
271 )
273 crps.rename(f"CRPS_of_{cubes[0].name()}")
274 _realization_callback(crps)
275 return crps
278def _scores_categorical_metric(
279 cubes: CubeList,
280 preserved_coordinates: list[str] | str | None,
281 threshold: str,
282 op_func: str,
283 metric: str,
284) -> CubeList:
285 """
286 Prepare cubes for computing categorical metrics using Scores.
288 Parameters
289 ----------
290 cubes: iris.cube.CubeList
291 An iris cubelist containing model(s) and an observation cube.
292 preserved_coordinates: list | str | None
293 An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
294 then preserving coordinate 'time' will produce a probability of detection score for each timeslice (shape time). If None,
295 then it will return a single value score for all times/point locations.
296 threshold: str
297 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).
298 op_func: str
299 A string either containing 'lt' for less than or 'gt for greater than, to determine how the threshold is applied to the data
300 to generate the mask.
301 metric: str
302 The scores metric to compute.
304 Returns
305 -------
306 scores_results: iris.cube.CubeList
307 An iris cubelist, containing the scores metric for each model for further plotting.
309 """
310 # Split obs/models
311 models = CubeList()
312 for c in cubes:
313 if "observed" in c.long_name:
314 observed = c
315 else:
316 models.append(c)
318 ops = {
319 "gt": operator.gt,
320 "lt": operator.lt,
321 }
323 # Check if this exists.
324 try:
325 op = ops[op_func]
326 except KeyError as err:
327 raise ValueError(f"Operator {op_func} not supported.") from err
329 scores_results = CubeList()
331 for model in models:
332 # Convert obs cubes to xarray and resolve preserved dimensions.
333 other_xr = xr.DataArray.from_iris(model)
334 base_xr = xr.DataArray.from_iris(observed)
335 preserve_dims = _resolve_preserve_dims(
336 observed, other_xr, preserved_coordinates
337 )
339 # Create event operator object using threshold and operator direction.
340 event_operator = scores.categorical.ThresholdEventOperator(
341 default_event_threshold=float(threshold), default_op_fn=op
342 )
344 # Generate binary fields using the event operator.
345 forecast_binary, observed_binary = event_operator.make_event_tables(
346 other_xr, base_xr
347 )
349 # Create binary contigency manager, as per Scores API, using transform to preserve preserve_dims
350 contingency_manager = scores.categorical.BinaryContingencyManager(
351 forecast_binary, observed_binary
352 ).transform(preserve_dims=preserve_dims)
354 # Compute required categorical score.
355 if metric == "pod":
356 result = contingency_manager.probability_of_detection()
357 name = "Probability_Of_Detection"
359 elif metric == "ets":
360 result = contingency_manager.equitable_threat_score()
361 name = "Equitable_Threat_Score"
363 elif metric == "fb":
364 result = contingency_manager.frequency_bias()
365 name = "Frequency_Bias"
367 elif metric == "pfd": 367 ↛ 372line 367 didn't jump to line 372 because the condition on line 367 was always true
368 result = contingency_manager.probability_of_false_detection()
369 name = "Probability_Of_False_Detection"
371 else:
372 raise ValueError(f"Unknown metric {metric}")
374 scores_cube = xr.DataArray.to_iris(result)
375 _fix_spatial_coord(cubes[0], scores_cube)
376 scores_cube.rename(f"{name}_{op_func}_{threshold}_{observed.name()}")
377 scores_cube.units = "1"
378 scores_cube.attributes["model_name"] = model.attributes["model_name"]
380 scores_results.append(scores_cube)
382 return scores_results
385def scores_pod(
386 cubes,
387 preserved_coordinates,
388 threshold,
389 op_func,
390) -> CubeList:
391 r"""
392 Compute the Probability of Detection (POD) score using Scores ([scoresa]_ [scoresb]_).
394 Parameters
395 ----------
396 cubes: iris.cube.CubeList
397 An iris cubelist containing model(s) and an observation cube.
398 preserved_coordinates: list | str | None
399 An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
400 then preserving coordinate 'time' will produce the equitable threat score for each timeslice (shape time). If None,
401 then it will return a single value score for all times/point locations.
402 threshold: str
403 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).
404 op_func: str
405 A string either containing 'lt' for less than or 'gt for greater than, to determine how the threshold is applied to the data
406 to generate the mask.
408 Returns
409 -------
410 iris.cube.CubeList
411 An iris cubelist, containing the probability of detection score for each model for further plotting.
414 Notes
415 -----
416 The probability of detection calculates the proportion of observed events that meet a threshold that were correctly forecast by the model.
417 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
418 295k, that would be a positive hit. It does not take into account how far above/below a threshold a model forecasts.
420 It is calculated as .. math:: POD = \frac{true positives}{true positives + false negatives}
422 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.
424 POD produces a range of 0 to 1, where 1 is a perfect score.
425 """
426 return _scores_categorical_metric(
427 cubes,
428 preserved_coordinates,
429 threshold,
430 op_func,
431 "pod",
432 )
435def scores_ets(
436 cubes: CubeList,
437 preserved_coordinates: list[str] | str | None,
438 threshold: str,
439 op_func: str,
440) -> CubeList:
441 r"""
442 Compute the Equitable Threat Score (ETS) score using Scores ([scoresa]_ [scoresb]_).
444 Parameters
445 ----------
446 cubes: iris.cube.CubeList
447 An iris cubelist containing model(s) and an observation cube.
448 preserved_coordinates: list | str | None
449 An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
450 then preserving coordinate 'time' will produce the equitable threat score for each timeslice (shape time). If None,
451 then it will return a single value score for all times/point locations.
452 threshold: str
453 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).
454 op_func: str
455 A string either containing 'lt' for less than or 'gt for greater than, to determine how the threshold is applied to the data
456 to generate the mask.
458 Returns
459 -------
460 iris.cube.CubeList
461 An iris cubelist, containing the probability of detection score for each model for further plotting.
463 Notes
464 -----
465 The Equitable Threat Score (ETS) evaluates the accuracy of forecasts for events that meet a specified threshold,
466 hile accounting for correct forecasts that could occur purely by chance. Unlike the Probability of Detection (POD),
467 ETS considers hits, misses, and false alarms, providing a more balanced assessment of forecast skill.
469 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
470 would be counted as a hit. ETS adjusts the total number of hits by removing the number of hits expected due to random chance.
472 It is calculated as:
474 .. math::
476 ETS = \frac{hits - hits_{random}}
477 {hits + misses + false\ alarms - hits_{random}}
479 where
481 hits_{random} = \frac{(hits + misses)(hits + false\ alarms)}{total count}
483 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
484 random chance.
485 """
486 return _scores_categorical_metric(
487 cubes,
488 preserved_coordinates,
489 threshold,
490 op_func,
491 "ets",
492 )
495def scores_pfd(
496 cubes: CubeList,
497 preserved_coordinates: list[str] | str | None,
498 threshold: str,
499 op_func: str,
500) -> CubeList:
501 r"""
502 Compute the Probability of False Detection (PFD) score using Scores ([scoresa]_ [scoresb]_).
504 Parameters
505 ----------
506 cubes: iris.cube.CubeList
507 An iris cubelist containing model(s) and an observation cube.
508 preserved_coordinates: list | str | None
509 An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
510 then preserving coordinate 'time' will produce the probability of false detection score for each timeslice (shape time). If None,
511 then it will return a single value score for all times/point locations.
512 threshold: str
513 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).
514 op_func: str
515 A string either containing 'lt' for less than or 'gt' for greater than, to determine how the threshold is applied to the data
516 to generate the mask.
518 Returns
519 -------
520 iris.cube.CubeList
521 An iris cubelist, containing the probability of false detection score for each model for further plotting.
523 Notes
524 -----
525 The Probability of False Detection (PFD) measures the proportion of observed non-events
526 that were incorrectly forecast as events. It is a measure of the false alarm rate and
527 provides information on how often a forecast system predicts threshold exceedances when
528 none actually occurred.
530 For example, if the threshold is 290 K and op_func is gt (greater than), an observation
531 of 288 K and a forecast of 295 K would be considered a false alarm, since the model
532 predicted an event but the observed value did not exceed the threshold.
534 It is calculated as:
536 .. math::
538 POFD = \frac{false\ alarms}
539 {false\ alarms + true\ negatives}
541 where
543 false alarms
544 Number of occasions where an event was forecast but did not occur.
546 true negatives
547 Number of occasions where neither the forecast nor the observations
548 indicated an event.
550 PFD ranges from 0 to 1, where 0 indicates a perfect score with no false alarms,
551 and 1 indicates that every observed non-event was incorrectly forecast as an event.
553 Lower values are therefore better.
554 """
555 return _scores_categorical_metric(
556 cubes,
557 preserved_coordinates,
558 threshold,
559 op_func,
560 "pfd",
561 )
564def scores_frequency_bias(
565 cubes: CubeList,
566 preserved_coordinates: list[str] | str | None,
567 threshold: str,
568 op_func: str,
569) -> CubeList:
570 r"""
571 Compute the Frequency Bias (FB) score using Scores ([scoresa]_ [scoresb]_).
573 Parameters
574 ----------
575 cubes: iris.cube.CubeList
576 An iris cubelist containing model(s) and an observation cube.
577 preserved_coordinates: list | str | None
578 An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
579 then preserving coordinate 'time' will produce the frequency bias score for each timeslice (shape time). If None,
580 then it will return a single value score for all times/point locations.
581 threshold: str
582 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).
583 op_func: str
584 A string either containing 'lt' for less than or 'gt' for greater than, to determine how the threshold is applied to the data
585 to generate the mask.
587 Returns
588 -------
589 iris.cube.CubeList
590 An iris cube, containing the frequency bias score for each model for further plotting.
592 Notes
593 -----
594 The Frequency Bias (FB) measures whether a forecasting system predicts an
595 event too frequently or too infrequently compared to observations. Unlike
596 metrics such as the Probability of Detection (POD) or Equitable Threat Score (ETS),
597 Frequency Bias does not assess the accuracy of forecast locations or timings,
598 only the overall frequency with which events are forecast.
600 For example, if the threshold is 290 K and op_func is gt (greater than), and
601 events exceeding this threshold are forecast twice as often as they are observed,
602 the frequency bias would be approximately 2. Conversely, if events are forecast
603 only half as often as they occur, the frequency bias would be approximately 0.5.
605 It is calculated as:
607 .. math::
609 Frequency\ Bias = \frac{hits + false\ alarms}
610 {hits + misses}
612 where
614 hits
615 Number of occasions where an event was forecast and observed.
617 false alarms
618 Number of occasions where an event was forecast but not observed.
620 misses
621 Number of occasions where an event was observed but not forecast.
623 Frequency Bias ranges from 0 to infinity, where:
625 * 1 indicates the forecast predicts events at the correct frequency.
626 * Greater than 1 indicates overforecasting of events.
627 * Less than 1 indicates underforecasting of events.
629 A perfect frequency bias score is therefore 1, although a value of 1 does not
630 necessarily imply a skillful forecast, as hits and false alarms may compensate
631 for one another.
632 """
633 return _scores_categorical_metric(
634 cubes,
635 preserved_coordinates,
636 threshold,
637 op_func,
638 "fb",
639 )
642def _make_scores_cube(
643 base: Cube, other: Cube, metric: str, preserved_coordinates: list[str]
644) -> Cube:
645 r"""Make the scores cube using the given scores metric.
647 Parameters
648 ----------
649 base: iris.cube.Cube
650 The cube from the "analysis" or observed cube.
651 other: iris.cube.Cube
652 The cube from the model.
654 metric: str
655 The scores metric to compute.
657 preserved_coordinates: list[str]
658 The list of preserved coordinates given by a user.
660 Returns
661 -------
662 scores_cube: iris.cube.Cube
663 The cube containing the calculated scores metric.
666 """
667 other_xr = xr.DataArray.from_iris(other)
668 base_xr = xr.DataArray.from_iris(base)
669 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
671 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
672 # apply scores, and then transform it back.
673 if metric == "rmse":
674 scores_cube = xr.DataArray.to_iris(
675 scores.continuous.rmse(other_xr, base_xr, preserve_dims=preserve_dims)
676 )
677 scores_cube.rename(f"RMSE_of_{base.name()}")
678 elif metric == "mae":
679 scores_cube = xr.DataArray.to_iris(
680 scores.continuous.mae(other_xr, base_xr, preserve_dims=preserve_dims)
681 )
682 scores_cube.rename(f"MAE_of_{base.name()}")
683 elif metric == "additive_bias":
684 scores_cube = xr.DataArray.to_iris(
685 scores.continuous.additive_bias(
686 other_xr, base_xr, preserve_dims=preserve_dims
687 )
688 )
689 scores_cube.rename(f"Additive_Bias_of_{base.name()}")
690 elif metric == "pearson_correlation": 690 ↛ 698line 690 didn't jump to line 698 because the condition on line 690 was always true
691 scores_cube = xr.DataArray.to_iris(
692 scores.continuous.correlation.pearsonr(
693 other_xr, base_xr, preserve_dims=preserve_dims
694 )
695 )
696 scores_cube.rename(f"Pearson_Correlation_of_{base.name()}")
697 else:
698 raise ValueError(f"Scores Unknown metric: {metric}")
700 _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 _ensure_increasing_pressure_coordinates(cubes: CubeList) -> CubeList:
744 r"""Ensure the pressure coordinate is increasing.
746 Parameters
747 ----------
748 cubes: iris.cube.CubeList
749 A CubeList of n cubes
751 Returns
752 -------
753 Cubes: iris.cube.CubeList
754 The original cube list but where each cube is ensured to have an increasing pressure coordinate.
755 """
756 for cube in cubes:
757 try:
758 if len(cube.coord("pressure").points) > 2 and not is_increasing(
759 cube.coord("pressure").points
760 ):
761 reverse(cube, "pressure")
763 except iris.exceptions.CoordinateNotFoundError:
764 pass
767def _process_cubes_for_verification(base: Cube, other: Cube) -> tuple[Cube, Cube]:
768 r"""Prepare cubes ready for verification in scores.
770 Parameters
771 ----------
772 base: iris.cube.Cube
773 The cube from the "analysis" or observed cube.
774 other: iris.cube.Cube
775 The cube from the model.
777 Returns
778 -------
779 base: iris.cube.Cube
780 The cube from the "analysis" or observed cube in the same format as the other model.
781 other: iris.cube.Cube
782 The cube from the model in the same format as the base model.
784 Notes
785 -----
786 This operator is used for sorting the data into the correct format. It
787 is likely going to need to be refactored out of CSET and perhaps moved into
788 `CSET._utils` given common code between here and `misc.difference`.
789 """
790 # Set cubes into correct format using code from difference operator
792 # Extract just common time points.
793 other_model_name = other.attributes["model_name"]
795 base, other = _extract_common_time_points(base, other)
797 if base.coords("pressure") and other.coords("pressure"):
798 common_pressure = np.intersect1d(
799 base.coord("pressure").points,
800 other.coord("pressure").points,
801 )
803 pressure_constraint = iris.Constraint(pressure=common_pressure.tolist())
805 base = base.extract(pressure_constraint)
806 other = other.extract(pressure_constraint)
808 # Get spatial coord names.
809 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
810 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
812 # Ensure cubes to compare are on common differencing grid.
813 # This is triggered if either
814 # i) latitude and longitude shapes are not the same. Note grid points
815 # are not compared directly as these can differ through rounding
816 # errors.
817 # ii) or variables are known to often sit on different grid staggering
818 # in different models (e.g. cell center vs cell edge), as is the case
819 # for UM and LFRic comparisons.
820 # In future greater choice of regridding method might be applied depending
821 # on variable type. Linear regridding can in general be appropriate for smooth
822 # variables. Care should be taken with interpretation of differences
823 # given this dependency on regridding.
824 if (
825 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
826 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
827 ) or (
828 base.long_name
829 in [
830 "eastward_wind_at_10m",
831 "northward_wind_at_10m",
832 "northward_wind_at_cell_centres",
833 "eastward_wind_at_cell_centres",
834 "zonal_wind_at_pressure_levels",
835 "meridional_wind_at_pressure_levels",
836 "potential_vorticity_at_pressure_levels",
837 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
838 ]
839 ):
840 logger.debug("Linear regridding base cube to other grid to compute differences")
841 base = regrid_onto_cube(base, other, method="Linear")
843 # Figure out if we are comparing between UM and LFRic; flip array if so.
844 base_lat_direction = is_increasing(base.coord(base_lat_name).points)
845 other_lat_direction = is_increasing(other.coord(other_lat_name).points)
846 if base_lat_direction != other_lat_direction: 846 ↛ 848line 846 didn't jump to line 848 because the condition on line 846 was never true
847 # Copy base cube for correct coordinate information.
848 other_tmp = base.copy()
849 # Flip the data and place in the copied cube.
850 other_tmp.data = np.flip(
851 other.data, other.coord(other_lat_name).cube_dims(other)
852 )
853 # Use original name and units from the other cube.
854 other_tmp.rename(other.name())
855 other_tmp.units = other.units
856 # Replace the cube.
857 other = other_tmp
859 # Equalise attributes so we can merge.
860 fully_equalise_attributes(CubeList([base, other]))
862 other.attributes["model_name"] = other_model_name
863 logger.debug("Base: %s\nOther: %s", base, other)
865 return base, other
868def _resolve_preserve_dims(
869 cube: Cube,
870 data_array: xr.DataArray,
871 preserved_coordinates: list[str] | str | None,
872) -> list[str] | None:
873 r"""Resolve preserve coordinates to xarray dimension names.
875 The ``scores`` package expects preserve dimensions to match xarray
876 dimension names. In Iris data, commonly used coordinates such as ``time``
877 may be auxiliary coordinates attached to a differently named dimension
878 (e.g. ``dim0``). This helper maps coordinate names to their underlying
879 dimension names and helps to convert from iris to xarray coordinate dimension names.
881 Parameters
882 ----------
883 cube: iris.cube.Cube
884 The ccomparison model cube.
885 data_array: xr.DataArray
886 The comparison model cube, but in xarray format.
887 preserved_coordinates: list[str]
888 The list of preserved coordinates given by a user.
891 Returns
892 -------
893 preserved_dims : list[str]
894 List of preserved dimension names in xarray convention.
896 """
897 if preserved_coordinates is None:
898 return None
900 coord_names = (
901 [preserved_coordinates]
902 if isinstance(preserved_coordinates, str)
903 else preserved_coordinates
904 )
905 preserve_dims: list[str] = []
907 for coord_name in coord_names:
908 # Already an xarray dimension name.
909 if coord_name in data_array.dims:
910 if coord_name not in preserve_dims: 910 ↛ 912line 910 didn't jump to line 912 because the condition on line 910 was always true
911 preserve_dims.append(coord_name)
912 continue
914 # Otherwise, map coordinate name to dimension index/indices.
915 try:
916 dim_indices = cube.coord_dims(coord_name)
917 except iris.exceptions.CoordinateNotFoundError:
918 # Keep original name so scores raises a clear error for unknown keys.
919 if coord_name not in preserve_dims:
920 preserve_dims.append(coord_name)
921 continue
923 for dim_index in dim_indices:
924 dim_name = data_array.dims[dim_index]
925 if dim_name not in preserve_dims:
926 preserve_dims.append(dim_name)
928 return preserve_dims
931def _attach_scaler_time_coord_maybe(scores_cube: Cube, base: Cube) -> None:
932 r"""Attaches scaler time coordinate if time is aggregated out.
934 In place function that attaches a scaler time coordinate to scores_cube
935 if time is aggregated out so plotting can display the aggregated period in the title.
937 Parameters
938 ----------
939 scores_cube: iris.cube.Cube
940 The calculated scores cube.
941 base: iris.cube.Cube
942 The cube from the "analysis" or observed cube.
944 Returns
945 -------
946 None
948 """
949 try:
950 if not scores_cube.coords("time"):
951 base_time = base.coord("time")
952 time_vals = (
953 base_time.bounds.flatten()
954 if base_time.has_bounds()
955 else base_time.points
956 )
957 t_start = float(time_vals[0])
958 t_end = float(time_vals[-1])
959 t_mid = 0.5 * (t_start + t_end)
961 scores_cube.add_aux_coord(
962 iris.coords.AuxCoord(
963 t_mid,
964 standard_name=base_time.standard_name,
965 long_name=base_time.long_name,
966 var_name=base_time.var_name,
967 units=base_time.units,
968 bounds=np.array([t_start, t_end]),
969 attributes=base_time.attributes.copy(),
970 )
971 )
972 except iris.exceptions.CoordinateNotFoundError:
973 pass
976def _split_base_and_other(cubes: CubeList):
977 r"""Split the cube into base and other cubes.
979 Split depends on whether there
980 is an observed cube in the cubes. If there is an observed cube,
981 then 'base' is the observed cube, if not then 'base' is the comparison
982 model cube.
984 Parameters
985 ----------
986 cubes: iris.cube.CubeList
987 Cubes to split into base and other cubes.
989 Returns
990 -------
991 tuple
992 A tuple containing a base cube, and other cube/cubelist.
994 """
995 obs_cube = [cb for cb in cubes if "observed" in (cb.long_name or "")]
996 if obs_cube:
997 if len(obs_cube) > 1: 997 ↛ 998line 997 didn't jump to line 998 because the condition on line 997 was never true
998 raise ValueError(
999 f"Expected exactly one 'observed' cube, found {len(obs_cube)}"
1000 )
1001 base = obs_cube[0]
1002 others = CubeList(cb for cb in cubes if cb is not base)
1003 return base, others
1005 return _sort_cube_into_base_and_other(cubes)