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