Coverage for src/CSET/operators/scoreswrappers.py: 82%

1# © Crown copyright, Met Office (2022-2026) and CSET contributors.

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

7# http://www.apache.org/licenses/LICENSE-2.0

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")

126

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

142

143 # Equalise attributes so we can merge.

144 fully_equalise_attributes(CubeList([base, other]))

145 logger.debug("Base: %s\nOther: %s", base, other)

146

147 return base, other

148

149

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.

156

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

165

166 coord_names = (

167 [preserved_coordinates]

168 if isinstance(preserved_coordinates, str)

169 else preserved_coordinates

170 )

171 preserve_dims: list[str] = []

172

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

179

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

188

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)

193

194 return preserve_dims

195

196

197def scores_rmse(cubes: CubeList, preserved_coordinates: list[str] | str | None = None):

198 r"""Calculate the Root Mean Square Error (RMSE) using scores.

199

200 Acts as a wrapper around the RMSE calculation from ``scores`` ([scoresa]_, [scoresb]_).

201 It is calculated as

202

203 .. math:: RMSE = \sqrt{\frac{1}{N} \Sigma(forecast - observations)^2}

204

205 Parameters

206 ----------

207 cubes: iris.cube.CubeList

208 A CubeList containing exactly two cubes: a base and an "other" model,

209 this can be an analysis and the model.

210 preserved_coordinates: list[str] | str | None, default is None.

211 The coordinates (or xarray dimension names) that you wish to preserve in the calculaiton of the

212 RMSE. For example if you want a map of each time you can preserve

213 ["time","grid_latitude", "grid_longitude"] or if you want a time series

214 you can preserve ["time"], if you want to collapse to a single value

215 use `None`. The default is `None`.

216

217 Returns

218 -------

219 scores_cube: iris.cube.Cube

220 A cube containing the RMSE between the base and other cube.

221

222 References

223 ----------

224 .. [scoresa] Leeuwenburg, T., Loveday, N., Ebert, E. E., Cook, H.,

225 Khanarmuei, M., Taggart, R. J., Ramanathan, N., Carroll, M., Chong, S.,

226 Griffiths, A., & Sharples, J. (2024) "scores: A Python package for

227 verifying and evaluating models and predictions with xarray". Journal

228 of Open Source Software, vol. 9, 6889. doi: 10.21105/joss.06889

229

230 .. [scoresb] Leeuwenburg, T., Loveday, N., Ramanathan, N., Chong, S.,

231 Taggart, R. J., Shrestha, D., Khanarmuei, M., Cook, H., Bluett, L., Ebert,

232 E. E., Carroll, M., Trotta, B., Bishop, S., Squire, D. T., Griffiths, A.,

233 Pagano, T. C., Fisher, A. J., Mandelbaum, T., Jinghan, F., … Smallwood, J.

234 (2026) "scores: Metrics for the verification, evaluation and optimisation of

235 forecasts, predictions or models (2.5.0)". Zenodo. doi: 10.5281/zenodo.18638494

236 """

237 base, other = _sort_cubes_for_verification(cubes)

238

239 # Copy the coordinates of the input cubes.

240 other_xr = xr.DataArray.from_iris(other)

241 base_xr = xr.DataArray.from_iris(base)

242 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)

243

244 # Scores operates on xarray data arrays, so we transform the iris cube into an array,

245 # apply scores, and then transform it back.

246 scores_cube = xr.DataArray.to_iris(

247 scores.continuous.rmse(

248 other_xr,

249 base_xr,

250 preserve_dims=preserve_dims,

251 )

252 )

253

254 # If time is aggregated out, attach a scalar time coordinate with bounds

255 # so plotting can display the aggregated period in the title.

256 try:

257 if not scores_cube.coords("time"):

258 base_time = base.coord("time")

259 time_vals = (

260 base_time.bounds.flatten()

261 if base_time.has_bounds()

262 else base_time.points

263 )

264 t_start = float(time_vals[0])

265 t_end = float(time_vals[-1])

266 t_mid = 0.5 * (t_start + t_end)

267

268 scores_cube.add_aux_coord(

269 iris.coords.AuxCoord(

270 t_mid,

271 standard_name=base_time.standard_name,

272 long_name=base_time.long_name,

273 var_name=base_time.var_name,

274 units=base_time.units,

275 bounds=np.array([t_start, t_end]),

276 attributes=base_time.attributes.copy(),

277 )

278 )

279 except iris.exceptions.CoordinateNotFoundError:

280 pass

281

282 scores_cube.rename(f"RMSE_of_{base.name()}")

283 # if preserved_coordinates == ["grid_latitude", "grid_longitude"]:

284 # scores_cube.add_aux_coord(time_coord)

285 return scores_cube

286

287

288def scores_mae(cubes: CubeList, preserved_coordinates: list[str] | str | None = None):

289 r"""Calculate the Mean Absolute Error (MAE) using scores.

290

291 Acts as a wrapper around the MAE calculation from ``scores`` ([scoresa]_, [scoresb]_).

292

293 Parameters

294 ----------

295 cubes: iris.cube.CubeList

296 A CubeList containing exactly two cubes: a base and an "other" model,

297 this can be an analysis and the model.

298 preserved_coordinates: list[str] | str | None, default is None.

299 The coordinates that you wish to preserve in the calculaiton of the

300 MAE. For example if you want a map of each time you can preserve

301 ["time","grid_latitude", "grid_longitude"] or if you want a time series

302 you can preserve ["time"], if you want to collapse to a single value

303 use `None`. The default is `None`.

304

305 Returns

306 -------

307 scores_cube: iris.cube.Cube

308 A cube containing the MAE between the base and other cube.

309

310 References

311 ----------

312 .. [scoresa] Leeuwenburg, T., Loveday, N., Ebert, E. E., Cook, H.,

313 Khanarmuei, M., Taggart, R. J., Ramanathan, N., Carroll, M., Chong, S.,

314 Griffiths, A., & Sharples, J. (2024) "scores: A Python package for

315 verifying and evaluating models and predictions with xarray". Journal

316 of Open Source Software, vol. 9, 6889. doi: 10.21105/joss.06889

317

318 .. [scoresb] Leeuwenburg, T., Loveday, N., Ramanathan, N., Chong, S.,

319 Taggart, R. J., Shrestha, D., Khanarmuei, M., Cook, H., Bluett, L., Ebert,

320 E. E., Carroll, M., Trotta, B., Bishop, S., Squire, D. T., Griffiths, A.,

321 Pagano, T. C., Fisher, A. J., Mandelbaum, T., Jinghan, F., … Smallwood, J.

322 (2026) "scores: Metrics for the verification, evaluation and optimisation of

323 forecasts, predictions or models (2.5.0)". Zenodo. doi: 10.5281/zenodo.18638494

324 """

325 base, other = _sort_cubes_for_verification(cubes)

326

327 # Copy the coordinates of the input cubes.

328 other_xr = xr.DataArray.from_iris(other)

329 base_xr = xr.DataArray.from_iris(base)

330 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)

331

332 # Scores operates on xarray data arrays, so we transform the iris cube into an array,

333 # apply scores, and then transform it back.

334 scores_cube = xr.DataArray.to_iris(

335 scores.continuous.mae(

336 other_xr,

337 base_xr,

338 preserve_dims=preserve_dims,

339 )

340 )

341

342 # If time is aggregated out, attach a scalar time coordinate with bounds

343 # so plotting can display the aggregated period in the title.

344 try:

345 if not scores_cube.coords("time"): 345 ↛ 370line 345 didn't jump to line 370 because the condition on line 345 was always true

346 base_time = base.coord("time")

347 time_vals = (

348 base_time.bounds.flatten()

349 if base_time.has_bounds()

350 else base_time.points

351 )

352 t_start = float(time_vals[0])

353 t_end = float(time_vals[-1])

354 t_mid = 0.5 * (t_start + t_end)

355

356 scores_cube.add_aux_coord(

357 iris.coords.AuxCoord(

358 t_mid,

359 standard_name=base_time.standard_name,

360 long_name=base_time.long_name,

361 var_name=base_time.var_name,

362 units=base_time.units,

363 bounds=np.array([t_start, t_end]),

364 attributes=base_time.attributes.copy(),

365 )

366 )

367 except iris.exceptions.CoordinateNotFoundError:

368 pass

369

370 scores_cube.rename(f"MAE_of_{base.name()}")

371 return scores_cube

372

373

374def scores_additive_bias(

375 cubes: CubeList, preserved_coordinates: list[str] | str | None = None

376):

377 r"""Calculate the Additive Bias (Mean Error) using scores.

378

379 Acts as a wrapper around the ME calculation from ``scores`` ([scoresa]_, [scoresb]_).

380

381 Parameters

382 ----------

383 cubes: iris.cube.CubeList

384 A CubeList containing exactly two cubes: a base and an "other" model,

385 this can be an analysis and the model.

386 preserved_coordinates: list[str] | str | None, default is None.

387 The coordinates that you wish to preserve in the calculaiton of the

388 ME. For example if you want a map of each time you can preserve

389 ["time","grid_latitude", "grid_longitude"] or if you want a time series

390 you can preserve ["time"], if you want to collapse to a single value

391 use `None`. The default is `None`.

392

393 Returns

394 -------

395 scores_cube: iris.cube.Cube

396 A cube containing the ME between the base and other cube.

397

398 References

399 ----------

400 .. [scoresa] Leeuwenburg, T., Loveday, N., Ebert, E. E., Cook, H.,

401 Khanarmuei, M., Taggart, R. J., Ramanathan, N., Carroll, M., Chong, S.,

402 Griffiths, A., & Sharples, J. (2024) "scores: A Python package for

403 verifying and evaluating models and predictions with xarray". Journal

404 of Open Source Software, vol. 9, 6889. doi: 10.21105/joss.06889

405

406 .. [scoresb] Leeuwenburg, T., Loveday, N., Ramanathan, N., Chong, S.,

407 Taggart, R. J., Shrestha, D., Khanarmuei, M., Cook, H., Bluett, L., Ebert,

408 E. E., Carroll, M., Trotta, B., Bishop, S., Squire, D. T., Griffiths, A.,

409 Pagano, T. C., Fisher, A. J., Mandelbaum, T., Jinghan, F., … Smallwood, J.

410 (2026) "scores: Metrics for the verification, evaluation and optimisation of

411 forecasts, predictions or models (2.5.0)". Zenodo. doi: 10.5281/zenodo.18638494

412 """

413 base, other = _sort_cubes_for_verification(cubes)

414

415 # Copy the coordinates of the input cubes.

416 other_xr = xr.DataArray.from_iris(other)

417 base_xr = xr.DataArray.from_iris(base)

418 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)

419

420 # Scores operates on xarray data arrays, so we transform the iris cube into an array,

421 # apply scores, and then transform it back.

422 scores_cube = xr.DataArray.to_iris(

423 scores.continuous.additive_bias(

424 other_xr,

425 base_xr,

426 preserve_dims=preserve_dims,

427 )

428 )

429

430 # If time is aggregated out, attach a scalar time coordinate with bounds

431 # so plotting can display the aggregated period in the title.

432 try:

433 if not scores_cube.coords("time"): 433 ↛ 457line 433 didn't jump to line 457 because the condition on line 433 was always true

434 base_time = base.coord("time")

435 time_vals = (

436 base_time.bounds.flatten()

437 if base_time.has_bounds()

438 else base_time.points

439 )

440 t_start = float(time_vals[0])

441 t_end = float(time_vals[-1])

442 t_mid = 0.5 * (t_start + t_end)

443

444 scores_cube.add_aux_coord(

445 iris.coords.AuxCoord(

446 t_mid,

447 standard_name=base_time.standard_name,

448 long_name=base_time.long_name,

449 var_name=base_time.var_name,

450 units=base_time.units,

451 bounds=np.array([t_start, t_end]),

452 attributes=base_time.attributes.copy(),

453 )

454 )

455 except iris.exceptions.CoordinateNotFoundError:

456 pass

457 scores_cube.rename(f"Additive_Bias_of_{base.name()}")

458 return scores_cube

459

460

461def scores_correlation_pearsonr(

462 cubes: CubeList, preserved_coordinates: list[str] | str | None = None

463):

464 r"""Calculate the Pearson's Correlation (PC) coefficient using scores.

465

466 Acts as a wrapper around the PC calculation from ``scores`` ([scoresa]_, [scoresb]_).

467

468 Parameters

469 ----------

470 cubes: iris.cube.CubeList

471 A CubeList containing exactly two cubes: a base and an "other" model,

472 this can be an analysis and the model.

473 preserved_coordinates: list[str] | str | None, default is None.

474 The coordinates that you wish to preserve in the calculation of the

475 PC. For example if you want a map of each time you can preserve

476 ["time","grid_latitude", "grid_longitude"] or if you want a time series

477 you can preserve ["time"], if you want to collapse to a single value

478 use `None`. The default is `None`.

479

480 Returns

481 -------

482 scores_cube: iris.cube.Cube

483 A cube containing the PC between the base and other cube.

484

485 References

486 ----------

487 .. [scoresa] Leeuwenburg, T., Loveday, N., Ebert, E. E., Cook, H.,

488 Khanarmuei, M., Taggart, R. J., Ramanathan, N., Carroll, M., Chong, S.,

489 Griffiths, A., & Sharples, J. (2024) "scores: A Python package for

490 verifying and evaluating models and predictions with xarray". Journal

491 of Open Source Software, vol. 9, 6889. doi: 10.21105/joss.06889

492

493 .. [scoresb] Leeuwenburg, T., Loveday, N., Ramanathan, N., Chong, S.,

494 Taggart, R. J., Shrestha, D., Khanarmuei, M., Cook, H., Bluett, L., Ebert,

495 E. E., Carroll, M., Trotta, B., Bishop, S., Squire, D. T., Griffiths, A.,

496 Pagano, T. C., Fisher, A. J., Mandelbaum, T., Jinghan, F., … Smallwood, J.

497 (2026) "scores: Metrics for the verification, evaluation and optimisation of

498 forecasts, predictions or models (2.5.0)". Zenodo. doi: 10.5281/zenodo.18638494

499 """

500 base, other = _sort_cubes_for_verification(cubes)

501

502 # Copy the coordinates of the input cubes.

503 other_xr = xr.DataArray.from_iris(other)

504 base_xr = xr.DataArray.from_iris(base)

505 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)

506

507 # Scores operates on xarray data arrays, so we transform the iris cube into an array,

508 # apply scores, and then transform it back.

509 scores_cube = xr.DataArray.to_iris(

510 scores.continuous.correlation.pearsonr(

511 other_xr,

512 base_xr,

513 preserve_dims=preserve_dims,

514 )

515 )

516

517 # If time is aggregated out, attach a scalar time coordinate with bounds

518 # so plotting can display the aggregated period in the title.

519 try:

520 if not scores_cube.coords("time"): 520 ↛ 545line 520 didn't jump to line 545 because the condition on line 520 was always true

521 base_time = base.coord("time")

522 time_vals = (

523 base_time.bounds.flatten()

524 if base_time.has_bounds()

525 else base_time.points

526 )

527 t_start = float(time_vals[0])

528 t_end = float(time_vals[-1])

529 t_mid = 0.5 * (t_start + t_end)

530

531 scores_cube.add_aux_coord(

532 iris.coords.AuxCoord(

533 t_mid,

534 standard_name=base_time.standard_name,

535 long_name=base_time.long_name,

536 var_name=base_time.var_name,

537 units=base_time.units,

538 bounds=np.array([t_start, t_end]),

539 attributes=base_time.attributes.copy(),

540 )

541 )

542 except iris.exceptions.CoordinateNotFoundError:

543 pass

544

545 scores_cube.rename(f"Pearson_Correlation_of_{base.name()}")

546 return scores_cube

547

548

549def scores_crps_for_ensemble(

550 cubes: Cube | CubeList, method: str = "ecdf", control_member: int = 0

551) -> iris.Constraint:

552 r"""Calculate the CRPS for an ensemble.

553

554 Acts as a wrapper around the crps_for_ensemble from ``scores`` ([scores_a]_, [scores_b]_).

555

556 Lower CRPS values are better (implies experiment distribution is closer to control distribution/observations),

557 larger values are worse (implies distributions are dissimilar).

558 It is applicable across time and spatial scales as the focus is on the distribution of the values.

559 Default method is ecdf. ecdf is exact value from the empirical distributions,

560 whereas fair produces an approximated value based on a random sample of the underlying distribution.

561

562 See [CRPS] for further information.

563

564 Parameters

565 ----------

566 cubes: iris.cube.Cube

567 A Cube containing ensembles data

568

569 Returns

570 -------

571 crps: iris.cube.Cube

572 A cube containing the crps between the ensemble members and the control

573

574 References

575 ----------

576 .. [scores_a] Leeuwenburg, T., Loveday, N., Ebert, E. E., Cook, H.,

577 Khanarmuei, M., Taggart, R. J., Ramanathan, N., Carroll, M., Chong, S.,

578 Griffiths, A., & Sharples, J. (2024) "scores: A Python package for

579 verifying and evaluating models and predictions with xarray". Journal

580 of Open Source Software, vol. 9, 6889. doi: 10.21105/joss.06889

581

582 .. [scores_b] Leeuwenburg, T., Loveday, N., Ramanathan, N., Chong, S.,

583 Taggart, R. J., Shrestha, D., Khanarmuei, M., Cook, H., Bluett, L., Ebert,

584 E. E., Carroll, M., Trotta, B., Bishop, S., Squire, D. T., Griffiths, A.,

585 Pagano, T. C., Fisher, A. J., Mandelbaum, T., Jinghan, F., … Smallwood, J.

586 (2026) "scores: Metrics for the verification, evaluation and optimisation of

587 forecasts, predictions or models (2.5.0)". Zenodo. doi: 10.5281/zenodo.18638494

588

589 .. [CRPS]

590 Hersbach, H., 2000: Decomposition of the Continuous Ranked

591 Probability Score for Ensemble Prediction Systems. Wea.

592 Forecasting, 15, 559–570, https://doi.org/10.1175/1520-0434(2000)015<0559:DOTCRP>2.0.CO;2.

593 """

594 if control_member != 0:

595 logger.warning("control member is usual 0")

596

597 if control_member not in cubes.coords("realization")[0].points:

598 new_control_member = cubes.coords("realization")[0].points[0]

599 logger.warning(

600 f"control member value {control_member} out of bounds, defaulting to control member={new_control_member}"

601 )

602 control_member = new_control_member

603

604 if cubes.coord("time").shape[0] == 1:

605 raise ValueError("Cube has only one time point.")

606

607 if cubes.coord("realization").shape[0] < 3:

608 raise ValueError("Cube should have one control member and at least two members")

609

610 ctrl = cubes.extract(generate_realization_constraint([control_member]))

611 ens_mem = cubes.extract(

612 generate_remove_single_ensemble_member_constraint(control_member)

613 )

614

615 # Realising the data in advance provides a large speedup

616 _ = ctrl.data

617 _ = ens_mem.data

618 del _

619

620 ctrl = xr.DataArray.from_iris(ctrl)

621 ens_mem = xr.DataArray.from_iris(ens_mem)

622

623 crps = xr.DataArray.to_iris(

624 scores.probability.crps_for_ensemble(

625 ens_mem,

626 ctrl,

627 ensemble_member_dim="realization",

628 method=method,

629 preserve_dims="time",

630 )

631 )

632

633 crps.rename(f"CRPS_of_{cubes[0].name()}")

634 _realization_callback(crps)

635 return crps