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