Coverage for src/CSET/operators/scoreswrappers.py: 83%
167 statements
« prev ^ index » next coverage.py v7.15.3, created at 2026-08-05 16:16 +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
19import iris
20import iris.exceptions
21import numpy as np
22import scores
23import scores.continuous
24import scores.probability
25import xarray as xr
26from iris.cube import Cube, CubeList
27from iris.util import reverse
29from CSET._common import is_increasing
30from CSET.operators._utils import fully_equalise_attributes, get_cube_yxcoordname
31from CSET.operators.constraints import (
32 generate_realization_constraint,
33 generate_remove_single_ensemble_member_constraint,
34)
35from CSET.operators.misc import _extract_common_time_points
36from CSET.operators.read import _realization_callback
37from CSET.operators.regrid import regrid_onto_cube
39logger = logging.getLogger(__name__)
42def _sort_cubes_for_verification(cubes: CubeList):
43 """Prepare cubes ready for verification in scores.
45 Parameters
46 ----------
47 cubes: iris.cube.CubeList
48 A CubeList of exact 2 cubes, one from each model.
50 Returns
51 -------
52 base: iris.cube.Cube
53 The cube from the "analysis" in the same format as the other model.
54 other: iris.cube.Cube
55 The cube from the model in the same format as the base model.
57 Raises
58 ------
59 ValueError: "cubes should contain exactly 2 cubes."
60 If any other number of cubes are present.
62 Notes
63 -----
64 This operator is used for sorting the data into the correct format. It
65 is likely going to need to be refactored out of CSET and perhaps moved into
66 `CSET._utils` given common code between here and `misc.difference`.
67 """
68 # Set cubes into correct format using code from difference operator
69 if len(cubes) != 2:
70 raise ValueError("cubes should contain exactly 2 cubes.")
71 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
72 other: Cube = cubes.extract_cube(
73 iris.Constraint(
74 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
75 )
76 )
78 # If cubes contain a pressure coordinate, ensure it is increasing.
79 for cube in cubes:
80 try:
81 if len(cube.coord("pressure").points) > 2 and not is_increasing( 81 ↛ 84line 81 didn't jump to line 84 because the condition on line 81 was never true
82 cube.coord("pressure").points
83 ):
84 reverse(cube, "pressure")
86 except iris.exceptions.CoordinateNotFoundError:
87 pass
89 # Extract just common time points.
90 base, other = _extract_common_time_points(base, other)
92 # Get spatial coord names.
93 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
94 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
96 # Ensure cubes to compare are on common differencing grid.
97 # This is triggered if either
98 # i) latitude and longitude shapes are not the same. Note grid points
99 # are not compared directly as these can differ through rounding
100 # errors.
101 # ii) or variables are known to often sit on different grid staggering
102 # in different models (e.g. cell center vs cell edge), as is the case
103 # for UM and LFRic comparisons.
104 # In future greater choice of regridding method might be applied depending
105 # on variable type. Linear regridding can in general be appropriate for smooth
106 # variables. Care should be taken with interpretation of differences
107 # given this dependency on regridding.
108 if (
109 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
110 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
111 ) or (
112 base.long_name
113 in [
114 "eastward_wind_at_10m",
115 "northward_wind_at_10m",
116 "northward_wind_at_cell_centres",
117 "eastward_wind_at_cell_centres",
118 "zonal_wind_at_pressure_levels",
119 "meridional_wind_at_pressure_levels",
120 "potential_vorticity_at_pressure_levels",
121 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
122 ]
123 ):
124 logger.debug("Linear regridding base cube to other grid to compute differences")
125 base = regrid_onto_cube(base, other, method="Linear")
127 # Figure out if we are comparing between UM and LFRic; flip array if so.
128 base_lat_direction = is_increasing(base.coord(base_lat_name).points)
129 other_lat_direction = is_increasing(other.coord(other_lat_name).points)
130 if base_lat_direction != other_lat_direction: 130 ↛ 132line 130 didn't jump to line 132 because the condition on line 130 was never true
131 # Copy base cube for correct coordinate information.
132 other_tmp = base.copy()
133 # Flip the data and place in the copied cube.
134 other_tmp.data = np.flip(
135 other.data, other.coord(other_lat_name).cube_dims(other)
136 )
137 # Use original name and units from the other cube.
138 other_tmp.rename(other.name())
139 other_tmp.units = other.units
140 # Replace the cube.
141 other = other_tmp
143 # Equalise attributes so we can merge.
144 fully_equalise_attributes(CubeList([base, other]))
145 logger.debug("Base: %s\nOther: %s", base, other)
147 return base, other
150def _resolve_preserve_dims(
151 cube: Cube,
152 data_array: xr.DataArray,
153 preserved_coordinates: list[str] | str | None,
154) -> list[str] | None:
155 """Resolve preserve coordinates to xarray dimension names.
157 The ``scores`` package expects preserve dimensions to match xarray
158 dimension names. In Iris data, commonly used coordinates such as ``time``
159 may be auxiliary coordinates attached to a differently named dimension
160 (e.g. ``dim0``). This helper maps coordinate names to their underlying
161 dimension names and helps to convert from iris to xarray coordinate dimension names.
162 """
163 if preserved_coordinates is None:
164 return None
166 coord_names = (
167 [preserved_coordinates]
168 if isinstance(preserved_coordinates, str)
169 else preserved_coordinates
170 )
171 preserve_dims: list[str] = []
173 for coord_name in coord_names:
174 # Already an xarray dimension name.
175 if coord_name in data_array.dims: 175 ↛ 181line 175 didn't jump to line 181 because the condition on line 175 was always true
176 if coord_name not in preserve_dims: 176 ↛ 178line 176 didn't jump to line 178 because the condition on line 176 was always true
177 preserve_dims.append(coord_name)
178 continue
180 # Otherwise, map coordinate name to dimension index/indices.
181 try:
182 dim_indices = cube.coord_dims(coord_name)
183 except iris.exceptions.CoordinateNotFoundError:
184 # Keep original name so scores raises a clear error for unknown keys.
185 if coord_name not in preserve_dims:
186 preserve_dims.append(coord_name)
187 continue
189 for dim_index in dim_indices:
190 dim_name = data_array.dims[dim_index]
191 if dim_name not in preserve_dims:
192 preserve_dims.append(dim_name)
194 return preserve_dims
197def _collapse_ensemble_mean(data_array: xr.DataArray) -> xr.DataArray:
198 """Collapse a realization dimension to its mean when present."""
199 if "realization" in data_array.dims:
200 return data_array.mean(dim="realization", keep_attrs=True)
202 return data_array
205def scores_rmse(cubes: CubeList, preserved_coordinates: list[str] | str | None = None):
206 r"""Calculate the Root Mean Square Error (RMSE) using scores.
208 Acts as a wrapper around the RMSE calculation from ``scores`` ([scoresa]_, [scoresb]_).
209 It is calculated as
211 .. math:: RMSE = \sqrt{\frac{1}{N} \Sigma(forecast - observations)^2}
213 Parameters
214 ----------
215 cubes: iris.cube.CubeList
216 A CubeList containing exactly two cubes: a base and an "other" model,
217 this can be an analysis and the model.
218 preserved_coordinates: list[str] | str | None, default is None.
219 The coordinates (or xarray dimension names) that you wish to preserve in the calculation of the
220 RMSE. For example if you want a map of each time you can preserve
221 ["time","grid_latitude", "grid_longitude"] or if you want a time series
222 you can preserve ["time"], if you want to collapse to a single value
223 use `None`. The default is `None`.
224 If an ensemble realization dimension is present, it is collapsed to the
225 ensemble mean before the RMSE is calculated.
227 Returns
228 -------
229 scores_cube: iris.cube.Cube
230 A cube containing the RMSE between the base and other cube.
232 References
233 ----------
234 .. [scoresa] Leeuwenburg, T., Loveday, N., Ebert, E. E., Cook, H.,
235 Khanarmuei, M., Taggart, R. J., Ramanathan, N., Carroll, M., Chong, S.,
236 Griffiths, A., & Sharples, J. (2024) "scores: A Python package for
237 verifying and evaluating models and predictions with xarray". Journal
238 of Open Source Software, vol. 9, 6889. doi: 10.21105/joss.06889
240 .. [scoresb] Leeuwenburg, T., Loveday, N., Ramanathan, N., Chong, S.,
241 Taggart, R. J., Shrestha, D., Khanarmuei, M., Cook, H., Bluett, L., Ebert,
242 E. E., Carroll, M., Trotta, B., Bishop, S., Squire, D. T., Griffiths, A.,
243 Pagano, T. C., Fisher, A. J., Mandelbaum, T., Jinghan, F., … Smallwood, J.
244 (2026) "scores: Metrics for the verification, evaluation and optimisation of
245 forecasts, predictions or models (2.5.0)". Zenodo. doi: 10.5281/zenodo.18638494
246 """
247 base, other = _sort_cubes_for_verification(cubes)
249 is_ensemble_mean = (
250 "realization" in xr.DataArray.from_iris(base).dims
251 or "realization" in xr.DataArray.from_iris(other).dims
252 )
254 # if ensemble data then calculate the ensemble mean first before calculating the RMSE
255 other_xr = _collapse_ensemble_mean(xr.DataArray.from_iris(other))
256 base_xr = _collapse_ensemble_mean(xr.DataArray.from_iris(base))
257 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
259 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
260 # apply scores, and then transform it back.
261 scores_cube = xr.DataArray.to_iris(
262 scores.continuous.rmse(
263 other_xr,
264 base_xr,
265 preserve_dims=preserve_dims,
266 )
267 )
269 # If time is aggregated out, attach a scalar time coordinate with bounds
270 # so plotting can display the aggregated period in the title.
271 try:
272 if not scores_cube.coords("time"):
273 base_time = base.coord("time")
274 time_vals = (
275 base_time.bounds.flatten()
276 if base_time.has_bounds()
277 else base_time.points
278 )
279 t_start = float(time_vals[0])
280 t_end = float(time_vals[-1])
281 t_mid = 0.5 * (t_start + t_end)
283 scores_cube.add_aux_coord(
284 iris.coords.AuxCoord(
285 t_mid,
286 standard_name=base_time.standard_name,
287 long_name=base_time.long_name,
288 var_name=base_time.var_name,
289 units=base_time.units,
290 bounds=np.array([t_start, t_end]),
291 attributes=base_time.attributes.copy(),
292 )
293 )
294 except iris.exceptions.CoordinateNotFoundError:
295 pass
296 scores_cube.rename(f"RMSE_of_{base.name()}")
297 # if preserved_coordinates == ["grid_latitude", "grid_longitude"]:
298 # scores_cube.add_aux_coord(time_coord)
300 # if ensemble add ensemble attribute
301 if is_ensemble_mean:
302 scores_cube.attributes["cset_ensemble_mean"] = "true"
303 return scores_cube
306def scores_mae(cubes: CubeList, preserved_coordinates: list[str] | str | None = None):
307 r"""Calculate the Mean Absolute Error (MAE) using scores.
309 Acts as a wrapper around the MAE calculation from ``scores`` ([scoresa]_, [scoresb]_).
311 Parameters
312 ----------
313 cubes: iris.cube.CubeList
314 A CubeList containing exactly two cubes: a base and an "other" model,
315 this can be an analysis and the model.
316 preserved_coordinates: list[str] | str | None, default is None.
317 The coordinates that you wish to preserve in the calculaiton of the
318 MAE. For example if you want a map of each time you can preserve
319 ["time","grid_latitude", "grid_longitude"] or if you want a time series
320 you can preserve ["time"], if you want to collapse to a single value
321 use `None`. The default is `None`.
323 Returns
324 -------
325 scores_cube: iris.cube.Cube
326 A cube containing the MAE between the base and other cube.
328 References
329 ----------
330 .. [scoresa] Leeuwenburg, T., Loveday, N., Ebert, E. E., Cook, H.,
331 Khanarmuei, M., Taggart, R. J., Ramanathan, N., Carroll, M., Chong, S.,
332 Griffiths, A., & Sharples, J. (2024) "scores: A Python package for
333 verifying and evaluating models and predictions with xarray". Journal
334 of Open Source Software, vol. 9, 6889. doi: 10.21105/joss.06889
336 .. [scoresb] Leeuwenburg, T., Loveday, N., Ramanathan, N., Chong, S.,
337 Taggart, R. J., Shrestha, D., Khanarmuei, M., Cook, H., Bluett, L., Ebert,
338 E. E., Carroll, M., Trotta, B., Bishop, S., Squire, D. T., Griffiths, A.,
339 Pagano, T. C., Fisher, A. J., Mandelbaum, T., Jinghan, F., … Smallwood, J.
340 (2026) "scores: Metrics for the verification, evaluation and optimisation of
341 forecasts, predictions or models (2.5.0)". Zenodo. doi: 10.5281/zenodo.18638494
342 """
343 base, other = _sort_cubes_for_verification(cubes)
345 # Copy the coordinates of the input cubes.
346 other_xr = xr.DataArray.from_iris(other)
347 base_xr = xr.DataArray.from_iris(base)
348 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
350 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
351 # apply scores, and then transform it back.
352 scores_cube = xr.DataArray.to_iris(
353 scores.continuous.mae(
354 other_xr,
355 base_xr,
356 preserve_dims=preserve_dims,
357 )
358 )
360 # If time is aggregated out, attach a scalar time coordinate with bounds
361 # so plotting can display the aggregated period in the title.
362 try:
363 if not scores_cube.coords("time"): 363 ↛ 387line 363 didn't jump to line 387 because the condition on line 363 was always true
364 base_time = base.coord("time")
365 time_vals = (
366 base_time.bounds.flatten()
367 if base_time.has_bounds()
368 else base_time.points
369 )
370 t_start = float(time_vals[0])
371 t_end = float(time_vals[-1])
372 t_mid = 0.5 * (t_start + t_end)
374 scores_cube.add_aux_coord(
375 iris.coords.AuxCoord(
376 t_mid,
377 standard_name=base_time.standard_name,
378 long_name=base_time.long_name,
379 var_name=base_time.var_name,
380 units=base_time.units,
381 bounds=np.array([t_start, t_end]),
382 attributes=base_time.attributes.copy(),
383 )
384 )
385 except iris.exceptions.CoordinateNotFoundError:
386 pass
387 scores_cube.rename(f"MAE_of_{base.name()}")
388 return scores_cube
391def scores_additive_bias(
392 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
393):
394 r"""Calculate the Additive Bias (Mean Error) using scores.
396 Acts as a wrapper around the ME calculation from ``scores`` ([scoresa]_, [scoresb]_).
398 Parameters
399 ----------
400 cubes: iris.cube.CubeList
401 A CubeList containing exactly two cubes: a base and an "other" model,
402 this can be an analysis and the model.
403 preserved_coordinates: list[str] | str | None, default is None.
404 The coordinates that you wish to preserve in the calculaiton of the
405 ME. For example if you want a map of each time you can preserve
406 ["time","grid_latitude", "grid_longitude"] or if you want a time series
407 you can preserve ["time"], if you want to collapse to a single value
408 use `None`. The default is `None`.
410 Returns
411 -------
412 scores_cube: iris.cube.Cube
413 A cube containing the ME between the base and other cube.
415 References
416 ----------
417 .. [scoresa] Leeuwenburg, T., Loveday, N., Ebert, E. E., Cook, H.,
418 Khanarmuei, M., Taggart, R. J., Ramanathan, N., Carroll, M., Chong, S.,
419 Griffiths, A., & Sharples, J. (2024) "scores: A Python package for
420 verifying and evaluating models and predictions with xarray". Journal
421 of Open Source Software, vol. 9, 6889. doi: 10.21105/joss.06889
423 .. [scoresb] Leeuwenburg, T., Loveday, N., Ramanathan, N., Chong, S.,
424 Taggart, R. J., Shrestha, D., Khanarmuei, M., Cook, H., Bluett, L., Ebert,
425 E. E., Carroll, M., Trotta, B., Bishop, S., Squire, D. T., Griffiths, A.,
426 Pagano, T. C., Fisher, A. J., Mandelbaum, T., Jinghan, F., … Smallwood, J.
427 (2026) "scores: Metrics for the verification, evaluation and optimisation of
428 forecasts, predictions or models (2.5.0)". Zenodo. doi: 10.5281/zenodo.18638494
429 """
430 base, other = _sort_cubes_for_verification(cubes)
432 # Copy the coordinates of the input cubes.
433 other_xr = xr.DataArray.from_iris(other)
434 base_xr = xr.DataArray.from_iris(base)
435 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
437 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
438 # apply scores, and then transform it back.
439 scores_cube = xr.DataArray.to_iris(
440 scores.continuous.additive_bias(
441 other_xr,
442 base_xr,
443 preserve_dims=preserve_dims,
444 )
445 )
447 # If time is aggregated out, attach a scalar time coordinate with bounds
448 # so plotting can display the aggregated period in the title.
449 try:
450 if not scores_cube.coords("time"): 450 ↛ 474line 450 didn't jump to line 474 because the condition on line 450 was always true
451 base_time = base.coord("time")
452 time_vals = (
453 base_time.bounds.flatten()
454 if base_time.has_bounds()
455 else base_time.points
456 )
457 t_start = float(time_vals[0])
458 t_end = float(time_vals[-1])
459 t_mid = 0.5 * (t_start + t_end)
461 scores_cube.add_aux_coord(
462 iris.coords.AuxCoord(
463 t_mid,
464 standard_name=base_time.standard_name,
465 long_name=base_time.long_name,
466 var_name=base_time.var_name,
467 units=base_time.units,
468 bounds=np.array([t_start, t_end]),
469 attributes=base_time.attributes.copy(),
470 )
471 )
472 except iris.exceptions.CoordinateNotFoundError:
473 pass
474 scores_cube.rename(f"Additive_Bias_of_{base.name()}")
475 return scores_cube
478def scores_correlation_pearsonr(
479 cubes: CubeList, preserved_coordinates: list[str] | str | None = None
480):
481 r"""Calculate the Pearson's Correlation (PC) coefficient using scores.
483 Acts as a wrapper around the PC calculation from ``scores`` ([scoresa]_, [scoresb]_).
485 Parameters
486 ----------
487 cubes: iris.cube.CubeList
488 A CubeList containing exactly two cubes: a base and an "other" model,
489 this can be an analysis and the model.
490 preserved_coordinates: list[str] | str | None, default is None.
491 The coordinates that you wish to preserve in the calculation of the
492 PC. For example if you want a map of each time you can preserve
493 ["time","grid_latitude", "grid_longitude"] or if you want a time series
494 you can preserve ["time"], if you want to collapse to a single value
495 use `None`. The default is `None`.
497 Returns
498 -------
499 scores_cube: iris.cube.Cube
500 A cube containing the PC between the base and other cube.
502 References
503 ----------
504 .. [scoresa] Leeuwenburg, T., Loveday, N., Ebert, E. E., Cook, H.,
505 Khanarmuei, M., Taggart, R. J., Ramanathan, N., Carroll, M., Chong, S.,
506 Griffiths, A., & Sharples, J. (2024) "scores: A Python package for
507 verifying and evaluating models and predictions with xarray". Journal
508 of Open Source Software, vol. 9, 6889. doi: 10.21105/joss.06889
510 .. [scoresb] Leeuwenburg, T., Loveday, N., Ramanathan, N., Chong, S.,
511 Taggart, R. J., Shrestha, D., Khanarmuei, M., Cook, H., Bluett, L., Ebert,
512 E. E., Carroll, M., Trotta, B., Bishop, S., Squire, D. T., Griffiths, A.,
513 Pagano, T. C., Fisher, A. J., Mandelbaum, T., Jinghan, F., … Smallwood, J.
514 (2026) "scores: Metrics for the verification, evaluation and optimisation of
515 forecasts, predictions or models (2.5.0)". Zenodo. doi: 10.5281/zenodo.18638494
516 """
517 base, other = _sort_cubes_for_verification(cubes)
519 # Copy the coordinates of the input cubes.
520 other_xr = xr.DataArray.from_iris(other)
521 base_xr = xr.DataArray.from_iris(base)
522 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
524 # Scores operates on xarray data arrays, so we transform the iris cube into an array,
525 # apply scores, and then transform it back.
526 scores_cube = xr.DataArray.to_iris(
527 scores.continuous.correlation.pearsonr(
528 other_xr,
529 base_xr,
530 preserve_dims=preserve_dims,
531 )
532 )
534 # If time is aggregated out, attach a scalar time coordinate with bounds
535 # so plotting can display the aggregated period in the title.
536 try:
537 if not scores_cube.coords("time"): 537 ↛ 562line 537 didn't jump to line 562 because the condition on line 537 was always true
538 base_time = base.coord("time")
539 time_vals = (
540 base_time.bounds.flatten()
541 if base_time.has_bounds()
542 else base_time.points
543 )
544 t_start = float(time_vals[0])
545 t_end = float(time_vals[-1])
546 t_mid = 0.5 * (t_start + t_end)
548 scores_cube.add_aux_coord(
549 iris.coords.AuxCoord(
550 t_mid,
551 standard_name=base_time.standard_name,
552 long_name=base_time.long_name,
553 var_name=base_time.var_name,
554 units=base_time.units,
555 bounds=np.array([t_start, t_end]),
556 attributes=base_time.attributes.copy(),
557 )
558 )
559 except iris.exceptions.CoordinateNotFoundError:
560 pass
562 scores_cube.rename(f"Pearson_Correlation_of_{base.name()}")
563 return scores_cube
566def scores_crps_for_ensemble(
567 cubes: Cube | CubeList, method: str = "ecdf", control_member: int = 0
568) -> iris.Constraint:
569 r"""Calculate the CRPS for an ensemble.
571 Acts as a wrapper around the crps_for_ensemble from ``scores`` ([scores_a]_, [scores_b]_).
573 Lower CRPS values are better (implies experiment distribution is closer to control distribution/observations),
574 larger values are worse (implies distributions are dissimilar).
575 It is applicable across time and spatial scales as the focus is on the distribution of the values.
576 Default method is ecdf. ecdf is exact value from the empirical distributions,
577 whereas fair produces an approximated value based on a random sample of the underlying distribution.
579 See [CRPS] for further information.
581 Parameters
582 ----------
583 cubes: iris.cube.Cube
584 A Cube containing ensembles data
586 Returns
587 -------
588 crps: iris.cube.Cube
589 A cube containing the crps between the ensemble members and the control
591 References
592 ----------
593 .. [scores_a] Leeuwenburg, T., Loveday, N., Ebert, E. E., Cook, H.,
594 Khanarmuei, M., Taggart, R. J., Ramanathan, N., Carroll, M., Chong, S.,
595 Griffiths, A., & Sharples, J. (2024) "scores: A Python package for
596 verifying and evaluating models and predictions with xarray". Journal
597 of Open Source Software, vol. 9, 6889. doi: 10.21105/joss.06889
599 .. [scores_b] Leeuwenburg, T., Loveday, N., Ramanathan, N., Chong, S.,
600 Taggart, R. J., Shrestha, D., Khanarmuei, M., Cook, H., Bluett, L., Ebert,
601 E. E., Carroll, M., Trotta, B., Bishop, S., Squire, D. T., Griffiths, A.,
602 Pagano, T. C., Fisher, A. J., Mandelbaum, T., Jinghan, F., … Smallwood, J.
603 (2026) "scores: Metrics for the verification, evaluation and optimisation of
604 forecasts, predictions or models (2.5.0)". Zenodo. doi: 10.5281/zenodo.18638494
606 .. [CRPS]
607 Hersbach, H., 2000: Decomposition of the Continuous Ranked
608 Probability Score for Ensemble Prediction Systems. Wea.
609 Forecasting, 15, 559–570, https://doi.org/10.1175/1520-0434(2000)015<0559:DOTCRP>2.0.CO;2.
610 """
611 if control_member != 0:
612 logger.warning("control member is usual 0")
614 if control_member not in cubes.coords("realization")[0].points:
615 new_control_member = cubes.coords("realization")[0].points[0]
616 logger.warning(
617 f"control member value {control_member} out of bounds, defaulting to control member={new_control_member}"
618 )
619 control_member = new_control_member
621 if cubes.coord("time").shape[0] == 1:
622 raise ValueError("Cube has only one time point.")
624 if cubes.coord("realization").shape[0] < 3:
625 raise ValueError("Cube should have one control member and at least two members")
627 ctrl = cubes.extract(generate_realization_constraint([control_member]))
628 ens_mem = cubes.extract(
629 generate_remove_single_ensemble_member_constraint(control_member)
630 )
632 # Realising the data in advance provides a large speedup
633 _ = ctrl.data
634 _ = ens_mem.data
635 del _
637 ctrl = xr.DataArray.from_iris(ctrl)
638 ens_mem = xr.DataArray.from_iris(ens_mem)
640 crps = xr.DataArray.to_iris(
641 scores.probability.crps_for_ensemble(
642 ens_mem,
643 ctrl,
644 ensemble_member_dim="realization",
645 method=method,
646 preserve_dims="time",
647 )
648 )
650 crps.rename(f"CRPS_of_{cubes[0].name()}")
651 _realization_callback(crps)
652 return crps