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

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

125

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

141

142 # Equalise attributes so we can merge.

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

144 logging.debug("Base: %s\nOther: %s", base, other)

145

146 return base, other

147

148

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.

155

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 and helps to convert from iris to xarray coordinate dimension names.

161 """

162 if preserved_coordinates is None:

163 return None

164

165 coord_names = (

166 [preserved_coordinates]

167 if isinstance(preserved_coordinates, str)

168 else preserved_coordinates

169 )

170 preserve_dims: list[str] = []

171

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

178

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

187

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)

192

193 return preserve_dims

194

195

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

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

198

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

200 It is calculated as

201

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

203

204 Parameters

205 ----------

206 cubes: iris.cube.CubeList

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

208 this can be an analysis and the model.

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

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

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

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

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

214 use `None`. The default is `None`.

215

216 Returns

217 -------

218 scores_cube: iris.cube.Cube

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

220

221 References

222 ----------

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

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

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

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

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

228

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

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

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

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

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

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

235 """

236 base, other = _sort_cubes_for_verification(cubes)

237

238 # Copy the coordinates of the input cubes.

239 other_xr = xr.DataArray.from_iris(other)

240 base_xr = xr.DataArray.from_iris(base)

241 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)

242

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

244 # apply scores, and then transform it back.

245 scores_cube = xr.DataArray.to_iris(

246 scores.continuous.rmse(

247 other_xr,

248 base_xr,

249 preserve_dims=preserve_dims,

250 )

251 )

252

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

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

255 try:

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

257 base_time = base.coord("time")

258 time_vals = (

259 base_time.bounds.flatten()

260 if base_time.has_bounds()

261 else base_time.points

262 )

263 t_start = float(time_vals[0])

264 t_end = float(time_vals[-1])

265 t_mid = 0.5 * (t_start + t_end)

266

267 scores_cube.add_aux_coord(

268 iris.coords.AuxCoord(

269 t_mid,

270 standard_name=base_time.standard_name,

271 long_name=base_time.long_name,

272 var_name=base_time.var_name,

273 units=base_time.units,

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

275 attributes=base_time.attributes.copy(),

276 )

277 )

278 except iris.exceptions.CoordinateNotFoundError:

279 pass

280

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

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

283 # scores_cube.add_aux_coord(time_coord)

284 return scores_cube

285

286

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

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

289

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

291

292 Parameters

293 ----------

294 cubes: iris.cube.CubeList

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

296 this can be an analysis and the model.

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

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

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

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

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

302 use `None`. The default is `None`.

303

304 Returns

305 -------

306 scores_cube: iris.cube.Cube

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

308

309 References

310 ----------

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

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

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

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

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

316

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

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

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

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

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

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

323 """

324 base, other = _sort_cubes_for_verification(cubes)

325

326 # Copy the coordinates of the input cubes.

327 other_xr = xr.DataArray.from_iris(other)

328 base_xr = xr.DataArray.from_iris(base)

329 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)

330

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

332 # apply scores, and then transform it back.

333 scores_cube = xr.DataArray.to_iris(

334 scores.continuous.mae(

335 other_xr,

336 base_xr,

337 preserve_dims=preserve_dims,

338 )

339 )

340

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

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

343 try:

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

345 base_time = base.coord("time")

346 time_vals = (

347 base_time.bounds.flatten()

348 if base_time.has_bounds()

349 else base_time.points

350 )

351 t_start = float(time_vals[0])

352 t_end = float(time_vals[-1])

353 t_mid = 0.5 * (t_start + t_end)

354

355 scores_cube.add_aux_coord(

356 iris.coords.AuxCoord(

357 t_mid,

358 standard_name=base_time.standard_name,

359 long_name=base_time.long_name,

360 var_name=base_time.var_name,

361 units=base_time.units,

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

363 attributes=base_time.attributes.copy(),

364 )

365 )

366 except iris.exceptions.CoordinateNotFoundError:

367 pass

368

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

370 return scores_cube

371

372

373def scores_additive_bias(

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

375):

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

377

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

379

380 Parameters

381 ----------

382 cubes: iris.cube.CubeList

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

384 this can be an analysis and the model.

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

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

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

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

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

390 use `None`. The default is `None`.

391

392 Returns

393 -------

394 scores_cube: iris.cube.Cube

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

396

397 References

398 ----------

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

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

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

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

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

404

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

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

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

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

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

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

411 """

412 base, other = _sort_cubes_for_verification(cubes)

413

414 # Copy the coordinates of the input cubes.

415 other_xr = xr.DataArray.from_iris(other)

416 base_xr = xr.DataArray.from_iris(base)

417 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)

418

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

420 # apply scores, and then transform it back.

421 scores_cube = xr.DataArray.to_iris(

422 scores.continuous.additive_bias(

423 other_xr,

424 base_xr,

425 preserve_dims=preserve_dims,

426 )

427 )

428

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

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

431 try:

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

433 base_time = base.coord("time")

434 time_vals = (

435 base_time.bounds.flatten()

436 if base_time.has_bounds()

437 else base_time.points

438 )

439 t_start = float(time_vals[0])

440 t_end = float(time_vals[-1])

441 t_mid = 0.5 * (t_start + t_end)

442

443 scores_cube.add_aux_coord(

444 iris.coords.AuxCoord(

445 t_mid,

446 standard_name=base_time.standard_name,

447 long_name=base_time.long_name,

448 var_name=base_time.var_name,

449 units=base_time.units,

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

451 attributes=base_time.attributes.copy(),

452 )

453 )

454 except iris.exceptions.CoordinateNotFoundError:

455 pass

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

457 return scores_cube

458

459

460def scores_correlation_pearsonr(

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

462):

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

464

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

466

467 Parameters

468 ----------

469 cubes: iris.cube.CubeList

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

471 this can be an analysis and the model.

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

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

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

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

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

477 use `None`. The default is `None`.

478

479 Returns

480 -------

481 scores_cube: iris.cube.Cube

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

483

484 References

485 ----------

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

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

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

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

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

491

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

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

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

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

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

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

498 """

499 base, other = _sort_cubes_for_verification(cubes)

500

501 # Copy the coordinates of the input cubes.

502 other_xr = xr.DataArray.from_iris(other)

503 base_xr = xr.DataArray.from_iris(base)

504 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)

505

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

507 # apply scores, and then transform it back.

508 scores_cube = xr.DataArray.to_iris(

509 scores.continuous.correlation.pearsonr(

510 other_xr,

511 base_xr,

512 preserve_dims=preserve_dims,

513 )

514 )

515

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

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

518 try:

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

520 base_time = base.coord("time")

521 time_vals = (

522 base_time.bounds.flatten()

523 if base_time.has_bounds()

524 else base_time.points

525 )

526 t_start = float(time_vals[0])

527 t_end = float(time_vals[-1])

528 t_mid = 0.5 * (t_start + t_end)

529

530 scores_cube.add_aux_coord(

531 iris.coords.AuxCoord(

532 t_mid,

533 standard_name=base_time.standard_name,

534 long_name=base_time.long_name,

535 var_name=base_time.var_name,

536 units=base_time.units,

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

538 attributes=base_time.attributes.copy(),

539 )

540 )

541 except iris.exceptions.CoordinateNotFoundError:

542 pass

543

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

545 return scores_cube

546

547

548def scores_crps_for_ensemble(

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

550) -> iris.Constraint:

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

552

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

554

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

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

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

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

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

560

561 See [CRPS] for further information.

562

563 Parameters

564 ----------

565 cubes: iris.cube.Cube

566 A Cube containing ensembles data

567

568 Returns

569 -------

570 crps: iris.cube.Cube

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

572

573 References

574 ----------

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

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

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

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

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

580

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

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

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

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

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

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

587

588 .. [CRPS]

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

590 Probability Score for Ensemble Prediction Systems. Wea.

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

592 """

593 if control_member != 0:

594 logging.warning("control member is usual 0")

595

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

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

598 logging.warning(

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

600 )

601 control_member = new_control_member

602

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

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

605

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

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

608

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

610 ens_mem = cubes.extract(

611 generate_remove_single_ensemble_member_constraint(control_member)

612 )

613

614 # Realising the data in advance provides a large speedup

615 _ = ctrl.data

616 _ = ens_mem.data

617 del _

618

619 ctrl = xr.DataArray.from_iris(ctrl)

620 ens_mem = xr.DataArray.from_iris(ens_mem)

621

622 crps = xr.DataArray.to_iris(

623 scores.probability.crps_for_ensemble(

624 ens_mem,

625 ctrl,

626 ensemble_member_dim="realization",

627 method=method,

628 preserve_dims="time",

629 )

630 )

631

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

633 _realization_callback(crps)

634 return crps