Coverage for src/CSET/operators/scoreswrappers.py: 89%
225 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-15 09:49 +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()
78 base, others = _split_base_and_other(cubes)
80 for other in others:
81 base, other = _process_cubes_for_verification(base, other)
83 scores_cube = _make_scores_cube(base, other, "rmse", preserved_coordinates)
85 scores_cubelist.append(scores_cube)
87 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
90def scores_mae(
91 cubes: CubeList,
92 preserved_coordinates: list[str] | str | None = None,
93) -> CubeList:
94 r"""Calculate the Mean Absolute Error (MAE) using scores.
96 Acts as a wrapper around the MAE calculation from ``scores`` ([scoresa]_, [scoresb]_).
98 Parameters
99 ----------
100 cubes: iris.cube.CubeList
101 A CubeList containing exactly two cubes: a base and an "other" model,
102 this can be an analysis and the model.
103 preserved_coordinates: list[str] | str | None, default is None.
104 The coordinates that you wish to preserve in the calculaiton of the
105 MAE. For example if you want a map of each time you can preserve
106 ["time","grid_latitude", "grid_longitude"] or if you want a time series
107 you can preserve ["time"], if you want to collapse to a single value
108 use `None`. The default is `None`.
110 Returns
111 -------
112 scores_cubelist: iris.cube.CubeList
113 A cubelist containing the MAE between the base and other cube(s).
114 """
115 scores_cubelist = CubeList()
116 base, others = _split_base_and_other(cubes)
118 for other in others:
119 base, other = _process_cubes_for_verification(base, other)
120 scores_cube = _make_scores_cube(base, other, "mae", preserved_coordinates)
121 scores_cubelist.append(scores_cube)
123 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
126def scores_additive_bias(
127 cubes: CubeList,
128 preserved_coordinates: list[str] | str | None = None,
129) -> CubeList:
130 r"""Calculate the Additive Bias (Mean Error) using scores.
132 Acts as a wrapper around the ME calculation from ``scores`` ([scoresa]_, [scoresb]_).
134 Parameters
135 ----------
136 cubes: iris.cube.CubeList
137 A CubeList containing exactly two cubes: a base and an "other" model,
138 this can be an analysis and the model.
139 preserved_coordinates: list[str] | str | None, default is None.
140 The coordinates that you wish to preserve in the calculaiton of the
141 ME. For example if you want a map of each time you can preserve
142 ["time","grid_latitude", "grid_longitude"] or if you want a time series
143 you can preserve ["time"], if you want to collapse to a single value
144 use `None`. The default is `None`.
146 Returns
147 -------
148 scores_cubelist: iris.cube.CubeList
149 A cubelist containing the ME between the base and other cube(s).
150 """
151 scores_cubelist = CubeList()
152 base, others = _split_base_and_other(cubes)
154 for other in others:
155 base, other = _process_cubes_for_verification(base, other)
157 scores_cube = _make_scores_cube(
158 base, other, "additive_bias", preserved_coordinates
159 )
161 scores_cubelist.append(scores_cube)
163 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
166def scores_correlation_pearsonr(
167 cubes: CubeList,
168 preserved_coordinates: list[str] | str | None = None,
169) -> CubeList:
170 r"""Calculate the Pearson's Correlation (PC) coefficient using scores.
172 Acts as a wrapper around the PC calculation from ``scores`` ([scoresa]_, [scoresb]_).
174 Parameters
175 ----------
176 cubes: iris.cube.CubeList
177 A CubeList containing exactly two cubes: a base and an "other" model,
178 this can be an analysis and the model.
179 preserved_coordinates: list[str] | str | None, default is None.
180 The coordinates that you wish to preserve in the calculation of the
181 PC. For example if you want a map of each time you can preserve
182 ["time","grid_latitude", "grid_longitude"] or if you want a time series
183 you can preserve ["time"], if you want to collapse to a single value
184 use `None`. The default is `None`.
186 Returns
187 -------
188 scores_cubelist: iris.cube.CubeList
189 A cubelist containing the PC between the base and other cube(s).
190 """
191 scores_cubelist = CubeList()
192 base, others = _split_base_and_other(cubes)
194 for other in others:
195 base, other = _process_cubes_for_verification(base, other)
197 scores_cube = _make_scores_cube(
198 base, other, "pearson_correlation", preserved_coordinates
199 )
201 scores_cubelist.append(scores_cube)
203 return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist
206def scores_crps_for_ensemble(
207 cubes: Cube | CubeList, method: str = "ecdf", control_member: int = 0
208) -> Cube:
209 r"""Calculate the CRPS for an ensemble.
211 Acts as a wrapper around the crps_for_ensemble from ``scores`` ([scoresa]_, [scoresb]_).
213 Lower CRPS values are better (implies experiment distribution is closer to control distribution/observations),
214 larger values are worse (implies distributions are dissimilar).
215 It is applicable across time and spatial scales as the focus is on the distribution of the values.
216 Default method is ecdf. ecdf is exact value from the empirical distributions,
217 whereas fair produces an approximated value based on a random sample of the underlying distribution.
219 See [CRPS]_ for further information.
221 Parameters
222 ----------
223 cubes: iris.cube.Cube
224 A Cube containing ensembles data
226 method: str ["ecfd" or "fair"]
227 Determines the method to use for calculating the CRPS. Defaults to "ecdf".
229 control_member: int
230 What the realisation the control member of the ensemble is. Defaults to 0. Sometimes this is 1.
232 Returns
233 -------
234 crps: iris.cube.Cube
235 A cube containing the crps between the ensemble members and the control
236 """
237 if control_member != 0:
238 logger.warning("control member is usual 0")
240 if control_member not in cubes.coords("realization")[0].points:
241 new_control_member = cubes.coords("realization")[0].points[0]
242 logger.warning(
243 f"control member value {control_member} out of bounds, defaulting to control member={new_control_member}"
244 )
245 control_member = new_control_member
247 if cubes.coord("time").shape[0] == 1:
248 raise ValueError("Cube has only one time point.")
250 if cubes.coord("realization").shape[0] < 3:
251 raise ValueError("Cube should have one control member and at least two members")
253 ctrl = cubes.extract(generate_realization_constraint([control_member]))
254 ens_mem = cubes.extract(
255 generate_remove_single_ensemble_member_constraint(control_member)
256 )
258 # Realising the data in advance provides a large speedup
259 _ = ctrl.data
260 _ = ens_mem.data
261 del _
263 ctrl = xr.DataArray.from_iris(ctrl)
264 ens_mem = xr.DataArray.from_iris(ens_mem)
266 crps = xr.DataArray.to_iris(
267 scores.probability.crps_for_ensemble(
268 ens_mem,
269 ctrl,
270 ensemble_member_dim="realization",
271 method=method,
272 preserve_dims="time",
273 )
274 )
276 crps.rename(f"CRPS_of_{cubes[0].name()}")
277 _realization_callback(crps)
278 return crps
281def _scores_categorical_metric(
282 cubes: CubeList,
283 preserved_coordinates: list[str] | str | None,
284 threshold: str,
285 op_func: str,
286 metric: str,
287) -> CubeList:
288 """
289 Prepare cubes for computing categorical metrics using Scores.
291 Parameters
292 ----------
293 cubes: iris.cube.CubeList
294 An iris cubelist containing model(s) and an observation cube.
295 preserved_coordinates: list | str | None
296 An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
297 then preserving coordinate 'time' will produce a probability of detection score for each timeslice (shape time). If None,
298 then it will return a single value score for all times/point locations.
299 threshold: str
300 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).
301 op_func: str
302 A string either containing 'lt' for less than or 'gt for greater than, to determine how the threshold is applied to the data
303 to generate the mask.
304 metric: str
305 The scores metric to compute.
307 Returns
308 -------
309 scores_results: iris.cube.CubeList
310 An iris cubelist, containing the scores metric for each model for further plotting.
312 """
313 # Split obs/models
314 models = CubeList()
315 for c in cubes:
316 if "observed" in c.long_name:
317 observed = c
318 else:
319 models.append(c)
321 ops = {
322 "gt": operator.gt,
323 "lt": operator.lt,
324 }
326 # Check if this exists.
327 try:
328 op = ops[op_func]
329 except KeyError as err:
330 raise ValueError(f"Operator {op_func} not supported.") from err
332 scores_results = CubeList()
334 for model in models:
335 # Convert obs cubes to xarray and resolve preserved dimensions.
336 other_xr = xr.DataArray.from_iris(model)
337 base_xr = xr.DataArray.from_iris(observed)
338 preserve_dims = _resolve_preserve_dims(
339 observed, other_xr, preserved_coordinates
340 )
342 # Create event operator object using threshold and operator direction.
343 event_operator = scores.categorical.ThresholdEventOperator(
344 default_event_threshold=float(threshold), default_op_fn=op
345 )
347 # Generate binary fields using the event operator.
348 forecast_binary, observed_binary = event_operator.make_event_tables(
349 other_xr, base_xr
350 )
352 # Create binary contigency manager, as per Scores API, using transform to preserve preserve_dims
353 contingency_manager = scores.categorical.BinaryContingencyManager(
354 forecast_binary, observed_binary
355 ).transform(preserve_dims=preserve_dims)
357 # Compute required categorical score.
358 if metric == "pod":
359 result = contingency_manager.probability_of_detection()
360 name = "Probability_Of_Detection"
362 elif metric == "ets":
363 result = contingency_manager.equitable_threat_score()
364 name = "Equitable_Threat_Score"
366 elif metric == "fb":
367 result = contingency_manager.frequency_bias()
368 name = "Frequency_Bias"
370 elif metric == "pfd": 370 ↛ 375line 370 didn't jump to line 375 because the condition on line 370 was always true
371 result = contingency_manager.probability_of_false_detection()
372 name = "Probability_Of_False_Detection"
374 else:
375 raise ValueError(f"Unknown metric {metric}")
377 scores_cube = xr.DataArray.to_iris(result)
378 _fix_spatial_coord(cubes[0], scores_cube)
379 scores_cube.rename(f"{name}_{op_func}_{threshold}_{observed.name()}")
380 scores_cube.units = "1"
381 scores_cube.attributes["model_name"] = model.attributes["model_name"]
383 scores_results.append(scores_cube)
385 return scores_results
388def scores_pod(
389 cubes,
390 preserved_coordinates,
391 threshold,
392 op_func,
393) -> CubeList:
394 r"""
395 Compute the Probability of Detection (POD) score using Scores ([scoresa]_ [scoresb]_).
397 Parameters
398 ----------
399 cubes: iris.cube.CubeList
400 An iris cubelist containing model(s) and an observation cube.
401 preserved_coordinates: list | str | None
402 An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
403 then preserving coordinate 'time' will produce the equitable threat score for each timeslice (shape time). If None,
404 then it will return a single value score for all times/point locations.
405 threshold: str
406 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).
407 op_func: str
408 A string either containing 'lt' for less than or 'gt for greater than, to determine how the threshold is applied to the data
409 to generate the mask.
411 Returns
412 -------
413 iris.cube.CubeList
414 An iris cubelist, containing the probability of detection score for each model for further plotting.
417 Notes
418 -----
419 The probability of detection calculates the proportion of observed events that meet a threshold that were correctly forecast by the model.
420 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
421 295k, that would be a positive hit. It does not take into account how far above/below a threshold a model forecasts.
423 It is calculated as .. math:: POD = \frac{true positives}{true positives + false negatives}
425 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.
427 POD produces a range of 0 to 1, where 1 is a perfect score.
428 """
429 return _scores_categorical_metric(
430 cubes,
431 preserved_coordinates,
432 threshold,
433 op_func,
434 "pod",
435 )
438def scores_ets(
439 cubes: CubeList,
440 preserved_coordinates: list[str] | str | None,
441 threshold: str,
442 op_func: str,
443) -> CubeList:
444 r"""
445 Compute the Equitable Threat Score (ETS) score using Scores ([scoresa]_ [scoresb]_).
447 Parameters
448 ----------
449 cubes: iris.cube.CubeList
450 An iris cubelist containing model(s) and an observation cube.
451 preserved_coordinates: list | str | None
452 An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
453 then preserving coordinate 'time' will produce the equitable threat score for each timeslice (shape time). If None,
454 then it will return a single value score for all times/point locations.
455 threshold: str
456 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).
457 op_func: str
458 A string either containing 'lt' for less than or 'gt for greater than, to determine how the threshold is applied to the data
459 to generate the mask.
461 Returns
462 -------
463 iris.cube.CubeList
464 An iris cubelist, containing the probability of detection score for each model for further plotting.
466 Notes
467 -----
468 The Equitable Threat Score (ETS) evaluates the accuracy of forecasts for events that meet a specified threshold,
469 hile accounting for correct forecasts that could occur purely by chance. Unlike the Probability of Detection (POD),
470 ETS considers hits, misses, and false alarms, providing a more balanced assessment of forecast skill.
472 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
473 would be counted as a hit. ETS adjusts the total number of hits by removing the number of hits expected due to random chance.
475 It is calculated as:
477 .. math::
479 ETS = \frac{hits - hits_{random}}
480 {hits + misses + false\ alarms - hits_{random}}
482 where
484 hits_{random} = \frac{(hits + misses)(hits + false\ alarms)}{total count}
486 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
487 random chance.
488 """
489 return _scores_categorical_metric(
490 cubes,
491 preserved_coordinates,
492 threshold,
493 op_func,
494 "ets",
495 )
498def scores_pfd(
499 cubes: CubeList,
500 preserved_coordinates: list[str] | str | None,
501 threshold: str,
502 op_func: str,
503) -> CubeList:
504 r"""
505 Compute the Probability of False Detection (PFD) score using Scores ([scoresa]_ [scoresb]_).
507 Parameters
508 ----------
509 cubes: iris.cube.CubeList
510 An iris cubelist containing model(s) and an observation cube.
511 preserved_coordinates: list | str | None
512 An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
513 then preserving coordinate 'time' will produce the probability of false detection score for each timeslice (shape time). If None,
514 then it will return a single value score for all times/point locations.
515 threshold: str
516 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).
517 op_func: str
518 A string either containing 'lt' for less than or 'gt' for greater than, to determine how the threshold is applied to the data
519 to generate the mask.
521 Returns
522 -------
523 iris.cube.CubeList
524 An iris cubelist, containing the probability of false detection score for each model for further plotting.
526 Notes
527 -----
528 The Probability of False Detection (PFD) measures the proportion of observed non-events
529 that were incorrectly forecast as events. It is a measure of the false alarm rate and
530 provides information on how often a forecast system predicts threshold exceedances when
531 none actually occurred.
533 For example, if the threshold is 290 K and op_func is gt (greater than), an observation
534 of 288 K and a forecast of 295 K would be considered a false alarm, since the model
535 predicted an event but the observed value did not exceed the threshold.
537 It is calculated as:
539 .. math::
541 POFD = \frac{false\ alarms}
542 {false\ alarms + true\ negatives}
544 where
546 false alarms
547 Number of occasions where an event was forecast but did not occur.
549 true negatives
550 Number of occasions where neither the forecast nor the observations
551 indicated an event.
553 PFD ranges from 0 to 1, where 0 indicates a perfect score with no false alarms,
554 and 1 indicates that every observed non-event was incorrectly forecast as an event.
556 Lower values are therefore better.
557 """
558 return _scores_categorical_metric(
559 cubes,
560 preserved_coordinates,
561 threshold,
562 op_func,
563 "pfd",
564 )
567def scores_frequency_bias(
568 cubes: CubeList,
569 preserved_coordinates: list[str] | str | None,
570 threshold: str,
571 op_func: str,
572) -> CubeList:
573 r"""
574 Compute the Frequency Bias (FB) score using Scores ([scoresa]_ [scoresb]_).
576 Parameters
577 ----------
578 cubes: iris.cube.CubeList
579 An iris cubelist containing model(s) and an observation cube.
580 preserved_coordinates: list | str | None
581 An object containing which coordinates to preserve in the computation. For example, if cubes contain shape time, point location,
582 then preserving coordinate 'time' will produce the frequency bias score for each timeslice (shape time). If None,
583 then it will return a single value score for all times/point locations.
584 threshold: str
585 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).
586 op_func: str
587 A string either containing 'lt' for less than or 'gt' for greater than, to determine how the threshold is applied to the data
588 to generate the mask.
590 Returns
591 -------
592 iris.cube.CubeList
593 An iris cube, containing the frequency bias score for each model for further plotting.
595 Notes
596 -----
597 The Frequency Bias (FB) measures whether a forecasting system predicts an
598 event too frequently or too infrequently compared to observations. Unlike
599 metrics such as the Probability of Detection (POD) or Equitable Threat Score (ETS),
600 Frequency Bias does not assess the accuracy of forecast locations or timings,
601 only the overall frequency with which events are forecast.
603 For example, if the threshold is 290 K and op_func is gt (greater than), and
604 events exceeding this threshold are forecast twice as often as they are observed,
605 the frequency bias would be approximately 2. Conversely, if events are forecast
606 only half as often as they occur, the frequency bias would be approximately 0.5.
608 It is calculated as:
610 .. math::
612 Frequency\ Bias = \frac{hits + false\ alarms}
613 {hits + misses}
615 where
617 hits
618 Number of occasions where an event was forecast and observed.
620 false alarms
621 Number of occasions where an event was forecast but not observed.
623 misses
624 Number of occasions where an event was observed but not forecast.
626 Frequency Bias ranges from 0 to infinity, where:
628 * 1 indicates the forecast predicts events at the correct frequency.
629 * Greater than 1 indicates overforecasting of events.
630 * Less than 1 indicates underforecasting of events.
632 A perfect frequency bias score is therefore 1, although a value of 1 does not
633 necessarily imply a skillful forecast, as hits and false alarms may compensate
634 for one another.
635 """
636 return _scores_categorical_metric(
637 cubes,
638 preserved_coordinates,
639 threshold,
640 op_func,
641 "fb",
642 )
645def _make_scores_cube(
646 base: Cube, other: Cube, metric: str, preserved_coordinates: list[str]
647) -> Cube:
648 r"""Make the scores cube using the given scores metric.
650 Parameters
651 ----------
652 base: iris.cube.Cube
653 The cube from the "analysis" or observed cube.
654 other: iris.cube.Cube
655 The cube from the model.
657 metric: str
658 The scores metric to compute.
660 preserved_coordinates: list[str]
661 The list of preserved coordinates given by a user.
663 Returns
664 -------
665 scores_cube: iris.cube.Cube
666 The cube containing the calculated scores metric.
669 """
670 other_xr = xr.DataArray.from_iris(other)
671 base_xr = xr.DataArray.from_iris(base)
672 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
674 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
675 # apply scores, and then transform it back.
676 if metric == "rmse":
677 scores_cube = xr.DataArray.to_iris(
678 scores.continuous.rmse(other_xr, base_xr, preserve_dims=preserve_dims)
679 )
680 scores_cube.rename(f"RMSE_of_{base.name()}")
681 elif metric == "mae":
682 scores_cube = xr.DataArray.to_iris(
683 scores.continuous.mae(other_xr, base_xr, preserve_dims=preserve_dims)
684 )
685 scores_cube.rename(f"MAE_of_{base.name()}")
686 elif metric == "additive_bias":
687 scores_cube = xr.DataArray.to_iris(
688 scores.continuous.additive_bias(
689 other_xr, base_xr, preserve_dims=preserve_dims
690 )
691 )
692 scores_cube.rename(f"Additive_Bias_of_{base.name()}")
693 elif metric == "pearson_correlation": 693 ↛ 701line 693 didn't jump to line 701 because the condition on line 693 was always true
694 scores_cube = xr.DataArray.to_iris(
695 scores.continuous.correlation.pearsonr(
696 other_xr, base_xr, preserve_dims=preserve_dims
697 )
698 )
699 scores_cube.rename(f"Pearson_Correlation_of_{base.name()}")
700 else:
701 raise ValueError(f"Scores Unknown metric: {metric}")
703 _fix_spatial_coord(base, scores_cube)
704 _attach_scaler_time_coord_maybe(scores_cube, base)
705 model_name = other.attributes["model_name"]
706 scores_cube.attributes["model_name"] = model_name
707 return scores_cube
710def _fix_spatial_coord(base_cube: Cube, scores_cube: Cube):
711 if not scores_cube.coords("grid_latitude") or not base_cube.coords("grid_latitude"):
712 return
713 if not scores_cube.coord("grid_latitude").coord_system: 713 ↛ exitline 713 didn't return from function '_fix_spatial_coord' because the condition on line 713 was always true
714 for xy_name in ("grid_latitude", "grid_longitude"):
715 scores_cube.coord(xy_name).coord_system = base_cube.coord(
716 xy_name
717 ).coord_system
720def _sort_cube_into_base_and_other(cubes: CubeList) -> tuple[Cube, CubeList]:
721 r"""Sorts cube into base and other models.
723 Parameters
724 ----------
725 cubes: iris.cube.CubeList
726 A CubeList of multiple cubes. One base cube and other model cubes.
728 Returns
729 -------
730 base: iris.cube.Cube
731 The cube from the "analysis" or observed cube in the same format as the other model.
732 others: iris.cube.CubeList
733 The cube list of containing the cube(s) from the model in the same format as the base model.
735 """
736 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
737 others: CubeList = cubes.extract(
738 iris.Constraint(
739 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
740 )
741 )
743 return base, others
746def _ensure_increasing_pressure_coordinates(cubes: CubeList) -> CubeList:
747 r"""Ensure the pressure coordinate is increasing.
749 Parameters
750 ----------
751 cubes: iris.cube.CubeList
752 A CubeList of n cubes
754 Returns
755 -------
756 Cubes: iris.cube.CubeList
757 The original cube list but where each cube is ensured to have an increasing pressure coordinate.
758 """
759 for cube in cubes:
760 try:
761 if len(cube.coord("pressure").points) > 2 and not is_increasing(
762 cube.coord("pressure").points
763 ):
764 reverse(cube, "pressure")
766 except iris.exceptions.CoordinateNotFoundError:
767 pass
770def _process_cubes_for_verification(base: Cube, other: Cube) -> tuple[Cube, Cube]:
771 r"""Prepare cubes ready for verification in scores.
773 Parameters
774 ----------
775 base: iris.cube.Cube
776 The cube from the "analysis" or observed cube.
777 other: iris.cube.Cube
778 The cube from the model.
780 Returns
781 -------
782 base: iris.cube.Cube
783 The cube from the "analysis" or observed cube in the same format as the other model.
784 other: iris.cube.Cube
785 The cube from the model in the same format as the base model.
787 Notes
788 -----
789 This operator is used for sorting the data into the correct format. It
790 is likely going to need to be refactored out of CSET and perhaps moved into
791 `CSET._utils` given common code between here and `misc.difference`.
792 """
793 # raises value error if cube time points don't match
794 ensure_common_time_points(base, other)
796 # Set cubes into correct format using code from difference operator
798 other_model_name = other.attributes["model_name"]
800 if base.coords("pressure") and other.coords("pressure"):
801 common_pressure = np.intersect1d(
802 base.coord("pressure").points,
803 other.coord("pressure").points,
804 )
806 pressure_constraint = iris.Constraint(pressure=common_pressure.tolist())
808 base = base.extract(pressure_constraint)
809 other = other.extract(pressure_constraint)
811 # Get spatial coord names.
812 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
813 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
815 # Ensure cubes to compare are on common differencing grid.
816 # This is triggered if either
817 # i) latitude and longitude shapes are not the same. Note grid points
818 # are not compared directly as these can differ through rounding
819 # errors.
820 # ii) or variables are known to often sit on different grid staggering
821 # in different models (e.g. cell center vs cell edge), as is the case
822 # for UM and LFRic comparisons.
823 # In future greater choice of regridding method might be applied depending
824 # on variable type. Linear regridding can in general be appropriate for smooth
825 # variables. Care should be taken with interpretation of differences
826 # given this dependency on regridding.
827 if (
828 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
829 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
830 ) or (
831 base.long_name
832 in [
833 "eastward_wind_at_10m",
834 "northward_wind_at_10m",
835 "northward_wind_at_cell_centres",
836 "eastward_wind_at_cell_centres",
837 "zonal_wind_at_pressure_levels",
838 "meridional_wind_at_pressure_levels",
839 "potential_vorticity_at_pressure_levels",
840 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
841 ]
842 ):
843 logger.debug("Linear regridding base cube to other grid to compute differences")
844 base = regrid_onto_cube(base, other, method="Linear")
846 # Figure out if we are comparing between UM and LFRic; flip array if so.
847 base_lat_direction = is_increasing(base.coord(base_lat_name).points)
848 other_lat_direction = is_increasing(other.coord(other_lat_name).points)
849 if base_lat_direction != other_lat_direction: 849 ↛ 851line 849 didn't jump to line 851 because the condition on line 849 was never true
850 # Copy base cube for correct coordinate information.
851 other_tmp = base.copy()
852 # Flip the data and place in the copied cube.
853 other_tmp.data = np.flip(
854 other.data, other.coord(other_lat_name).cube_dims(other)
855 )
856 # Use original name and units from the other cube.
857 other_tmp.rename(other.name())
858 other_tmp.units = other.units
859 # Replace the cube.
860 other = other_tmp
862 # Equalise attributes so we can merge.
863 fully_equalise_attributes(CubeList([base, other]))
865 other.attributes["model_name"] = other_model_name
866 logger.debug("Base: %s\nOther: %s", base, other)
868 return base, other
871def _resolve_preserve_dims(
872 cube: Cube,
873 data_array: xr.DataArray,
874 preserved_coordinates: list[str] | str | None,
875) -> list[str] | None:
876 r"""Resolve preserve coordinates to xarray dimension names.
878 The ``scores`` package expects preserve dimensions to match xarray
879 dimension names. In Iris data, commonly used coordinates such as ``time``
880 may be auxiliary coordinates attached to a differently named dimension
881 (e.g. ``dim0``). This helper maps coordinate names to their underlying
882 dimension names and helps to convert from iris to xarray coordinate dimension names.
884 Parameters
885 ----------
886 cube: iris.cube.Cube
887 The ccomparison model cube.
888 data_array: xr.DataArray
889 The comparison model cube, but in xarray format.
890 preserved_coordinates: list[str]
891 The list of preserved coordinates given by a user.
894 Returns
895 -------
896 preserved_dims : list[str]
897 List of preserved dimension names in xarray convention.
899 """
900 if preserved_coordinates is None:
901 return None
903 coord_names = (
904 [preserved_coordinates]
905 if isinstance(preserved_coordinates, str)
906 else preserved_coordinates
907 )
908 preserve_dims: list[str] = []
910 for coord_name in coord_names:
911 # Already an xarray dimension name.
912 if coord_name in data_array.dims:
913 if coord_name not in preserve_dims: 913 ↛ 915line 913 didn't jump to line 915 because the condition on line 913 was always true
914 preserve_dims.append(coord_name)
915 continue
917 # Otherwise, map coordinate name to dimension index/indices.
918 try:
919 dim_indices = cube.coord_dims(coord_name)
920 except iris.exceptions.CoordinateNotFoundError:
921 # Keep original name so scores raises a clear error for unknown keys.
922 if coord_name not in preserve_dims:
923 preserve_dims.append(coord_name)
924 continue
926 for dim_index in dim_indices:
927 dim_name = data_array.dims[dim_index]
928 if dim_name not in preserve_dims:
929 preserve_dims.append(dim_name)
931 return preserve_dims
934def _attach_scaler_time_coord_maybe(scores_cube: Cube, base: Cube) -> None:
935 r"""Attaches scaler time coordinate if time is aggregated out.
937 In place function that attaches a scaler time coordinate to scores_cube
938 if time is aggregated out so plotting can display the aggregated period in the title.
940 Parameters
941 ----------
942 scores_cube: iris.cube.Cube
943 The calculated scores cube.
944 base: iris.cube.Cube
945 The cube from the "analysis" or observed cube.
947 Returns
948 -------
949 None
951 """
952 try:
953 if not scores_cube.coords("time"):
954 base_time = base.coord("time")
955 time_vals = (
956 base_time.bounds.flatten()
957 if base_time.has_bounds()
958 else base_time.points
959 )
960 t_start = float(time_vals[0])
961 t_end = float(time_vals[-1])
962 t_mid = 0.5 * (t_start + t_end)
964 scores_cube.add_aux_coord(
965 iris.coords.AuxCoord(
966 t_mid,
967 standard_name=base_time.standard_name,
968 long_name=base_time.long_name,
969 var_name=base_time.var_name,
970 units=base_time.units,
971 bounds=np.array([t_start, t_end]),
972 attributes=base_time.attributes.copy(),
973 )
974 )
975 except iris.exceptions.CoordinateNotFoundError:
976 pass
979def _split_base_and_other(cubes: CubeList):
980 r"""Split the cube into base and other cubes.
982 Split depends on whether there
983 is an observed cube in the cubes. If there is an observed cube,
984 then 'base' is the observed cube, if not then 'base' is the comparison
985 model cube.
987 Parameters
988 ----------
989 cubes: iris.cube.CubeList
990 Cubes to split into base and other cubes.
992 Returns
993 -------
994 tuple
995 A tuple containing a base cube, and other cube/cubelist.
997 """
998 obs_cube = [cb for cb in cubes if "observed" in (cb.long_name or "")]
999 if obs_cube:
1000 if len(obs_cube) > 1: 1000 ↛ 1001line 1000 didn't jump to line 1001 because the condition on line 1000 was never true
1001 raise ValueError(
1002 f"Expected exactly one 'observed' cube, found {len(obs_cube)}"
1003 )
1004 base = obs_cube[0]
1005 others = CubeList(cb for cb in cubes if cb is not base)
1006 return base, others
1008 return _sort_cube_into_base_and_other(cubes)