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

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 xarray as xr

25from iris.cube import Cube, CubeList

26from iris.util import reverse

28from CSET._common import is_increasing

29from CSET.operators._utils import fully_equalise_attributes, get_cube_yxcoordname

30from CSET.operators.misc import _extract_common_time_points

31from CSET.operators.regrid import regrid_onto_cube

34def _sort_cubes_for_verification(cubes: CubeList):

35 """Prepare cubes ready for verification in scores.

37 Parameters

38 ----------

39 cubes: iris.cube.CubeList

40 A CubeList of exact 2 cubes, one from each model.

42 Returns

43 -------

44 base: iris.cube.Cube

45 The cube from the "analysis" in the same format as the other model.

46 other: iris.cube.Cube

47 The cube from the model in the same format as the base model.

49 Raises

50 ------

51 ValueError: "cubes should contain exactly 2 cubes."

52 If any other number of cubes are present.

54 Notes

55 -----

56 This operator is used for sorting the data into the correct format. It

57 is likely going to need to be refactored out of CSET and perhaps moved into

58 `CSET._utils` given common code between here and `misc.difference`.

59 """

60 # Set cubes into correct format using code from difference operator

61 if len(cubes) != 2:

62 raise ValueError("cubes should contain exactly 2 cubes.")

63 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))

64 other: Cube = cubes.extract_cube(

65 iris.Constraint(

66 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes

67 )

68 )

70 # If cubes contain a pressure coordinate, ensure it is increasing.

71 for cube in cubes:

72 try:

73 if len(cube.coord("pressure").points) > 2: 73 ↛ 71line 73 didn't jump to line 71 because the condition on line 73 was always true

74 if not is_increasing(cube.coord("pressure").points): 74 ↛ 75line 74 didn't jump to line 75 because the condition on line 74 was never true

75 reverse(cube, "pressure")

77 except iris.exceptions.CoordinateNotFoundError:

78 pass

80 # Extract just common time points.

81 base, other = _extract_common_time_points(base, other)

83 # Get spatial coord names.

84 base_lat_name, base_lon_name = get_cube_yxcoordname(base)

85 other_lat_name, other_lon_name = get_cube_yxcoordname(other)

87 # Ensure cubes to compare are on common differencing grid.

88 # This is triggered if either

89 # i) latitude and longitude shapes are not the same. Note grid points

90 # are not compared directly as these can differ through rounding

91 # errors.

92 # ii) or variables are known to often sit on different grid staggering

93 # in different models (e.g. cell center vs cell edge), as is the case

94 # for UM and LFRic comparisons.

95 # In future greater choice of regridding method might be applied depending

96 # on variable type. Linear regridding can in general be appropriate for smooth

97 # variables. Care should be taken with interpretation of differences

98 # given this dependency on regridding.

99 if (

100 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape

101 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape

102 ) or (

103 base.long_name

104 in [

105 "eastward_wind_at_10m",

106 "northward_wind_at_10m",

107 "northward_wind_at_cell_centres",

108 "eastward_wind_at_cell_centres",

109 "zonal_wind_at_pressure_levels",

110 "meridional_wind_at_pressure_levels",

111 "potential_vorticity_at_pressure_levels",

112 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",

113 ]

114 ):

115 logging.debug(

116 "Linear regridding base cube to other grid to compute differences"

117 )

118 base = regrid_onto_cube(base, other, method="Linear")

119

120 # Figure out if we are comparing between UM and LFRic; flip array if so.

121 base_lat_direction = is_increasing(base.coord(base_lat_name).points)

122 other_lat_direction = is_increasing(other.coord(other_lat_name).points)

123 if base_lat_direction != other_lat_direction: 123 ↛ 125line 123 didn't jump to line 125 because the condition on line 123 was never true

124 # Copy base cube for correct coordinate information.

125 other_tmp = base.copy()

126 # Flip the data and place in the copied cube.

127 other_tmp.data = np.flip(

128 other.data, other.coord(other_lat_name).cube_dims(other)

129 )

130 # Use original name and units from the other cube.

131 other_tmp.rename(other.name())

132 other_tmp.units = other.units

133 # Replace the cube.

134 other = other_tmp

135

136 # Equalise attributes so we can merge.

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

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

139

140 return base, other

141

142

143def _resolve_preserve_dims(

144 cube: Cube,

145 data_array: xr.DataArray,

146 preserved_coordinates: list[str] | str | None,

147) -> list[str] | None:

148 """Resolve preserve coordinates to xarray dimension names.

149

150 The ``scores`` package expects preserve dimensions to match xarray

151 dimension names. In Iris data, commonly used coordinates such as ``time``

152 may be auxiliary coordinates attached to a differently named dimension

153 (e.g. ``dim0``). This helper maps coordinate names to their underlying

154 dimension names and helps to convert from iris to xarray coordinate dimension names.

155 """

156 if preserved_coordinates is None:

157 return None

158

159 coord_names = (

160 [preserved_coordinates]

161 if isinstance(preserved_coordinates, str)

162 else preserved_coordinates

163 )

164 preserve_dims: list[str] = []

165

166 for coord_name in coord_names:

167 # Already an xarray dimension name.

168 if coord_name in data_array.dims: 168 ↛ 174line 168 didn't jump to line 174 because the condition on line 168 was always true

169 if coord_name not in preserve_dims: 169 ↛ 171line 169 didn't jump to line 171 because the condition on line 169 was always true

170 preserve_dims.append(coord_name)

171 continue

172

173 # Otherwise, map coordinate name to dimension index/indices.

174 try:

175 dim_indices = cube.coord_dims(coord_name)

176 except iris.exceptions.CoordinateNotFoundError:

177 # Keep original name so scores raises a clear error for unknown keys.

178 if coord_name not in preserve_dims:

179 preserve_dims.append(coord_name)

180 continue

181

182 for dim_index in dim_indices:

183 dim_name = data_array.dims[dim_index]

184 if dim_name not in preserve_dims:

185 preserve_dims.append(dim_name)

186

187 return preserve_dims

188

189

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

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

192

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

194 It is calculated as

195

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

197

198 Parameters

199 ----------

200 cubes: iris.cube.CubeList

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

202 this can be an analysis and the model.

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

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

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

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

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

208 use `None`. The default is `None`.

209

210 Returns

211 -------

212 scores_cube: iris.cube.Cube

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

214

215 References

216 ----------

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

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

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

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

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

222

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

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

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

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

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

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

229 """

230 base, other = _sort_cubes_for_verification(cubes)

231

232 # Copy the coordinates of the input cubes.

233 other_xr = xr.DataArray.from_iris(other)

234 base_xr = xr.DataArray.from_iris(base)

235 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)

236

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

238 # apply scores, and then transform it back.

239 scores_cube = xr.DataArray.to_iris(

240 scores.continuous.rmse(

241 other_xr,

242 base_xr,

243 preserve_dims=preserve_dims,

244 )

245 )

246

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

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

249 try:

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

251 base_time = base.coord("time")

252 time_vals = (

253 base_time.bounds.flatten()

254 if base_time.has_bounds()

255 else base_time.points

256 )

257 t_start = float(time_vals[0])

258 t_end = float(time_vals[-1])

259 t_mid = 0.5 * (t_start + t_end)

260

261 scores_cube.add_aux_coord(

262 iris.coords.AuxCoord(

263 t_mid,

264 standard_name=base_time.standard_name,

265 long_name=base_time.long_name,

266 var_name=base_time.var_name,

267 units=base_time.units,

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

269 attributes=base_time.attributes.copy(),

270 )

271 )

272 except iris.exceptions.CoordinateNotFoundError:

273 pass

274

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

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

277 # scores_cube.add_aux_coord(time_coord)

278 return scores_cube

279

280

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

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

283 base, other = _sort_cubes_for_verification(cubes)

284

285 # Copy the coordinates of the input cubes.

286 other_xr = xr.DataArray.from_iris(other)

287 base_xr = xr.DataArray.from_iris(base)

288 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)

289

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

291 # apply scores, and then transform it back.

292 scores_cube = xr.DataArray.to_iris(

293 scores.continuous.mae(

294 other_xr,

295 base_xr,

296 preserve_dims=preserve_dims,

297 )

298 )

299

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

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

302 try:

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

304 base_time = base.coord("time")

305 time_vals = (

306 base_time.bounds.flatten()

307 if base_time.has_bounds()

308 else base_time.points

309 )

310 t_start = float(time_vals[0])

311 t_end = float(time_vals[-1])

312 t_mid = 0.5 * (t_start + t_end)

313

314 scores_cube.add_aux_coord(

315 iris.coords.AuxCoord(

316 t_mid,

317 standard_name=base_time.standard_name,

318 long_name=base_time.long_name,

319 var_name=base_time.var_name,

320 units=base_time.units,

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

322 attributes=base_time.attributes.copy(),

323 )

324 )

325 except iris.exceptions.CoordinateNotFoundError:

326 pass

327

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

329 return scores_cube

330

331

332def scores_additive_bias(

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

334):

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

336 base, other = _sort_cubes_for_verification(cubes)

337

338 # Copy the coordinates of the input cubes.

339 other_xr = xr.DataArray.from_iris(other)

340 base_xr = xr.DataArray.from_iris(base)

341 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)

342

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

344 # apply scores, and then transform it back.

345 scores_cube = xr.DataArray.to_iris(

346 scores.continuous.additive_bias(

347 other_xr,

348 base_xr,

349 preserve_dims=preserve_dims,

350 )

351 )

352

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

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

355 try:

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

357 base_time = base.coord("time")

358 time_vals = (

359 base_time.bounds.flatten()

360 if base_time.has_bounds()

361 else base_time.points

362 )

363 t_start = float(time_vals[0])

364 t_end = float(time_vals[-1])

365 t_mid = 0.5 * (t_start + t_end)

366

367 scores_cube.add_aux_coord(

368 iris.coords.AuxCoord(

369 t_mid,

370 standard_name=base_time.standard_name,

371 long_name=base_time.long_name,

372 var_name=base_time.var_name,

373 units=base_time.units,

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

375 attributes=base_time.attributes.copy(),

376 )

377 )

378 except iris.exceptions.CoordinateNotFoundError:

379 pass

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

381 return scores_cube

382

383

384def scores_correlation_pearsonr(

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

386):

387 """Calculate the Pearson's Correlation using scores."""

388 base, other = _sort_cubes_for_verification(cubes)

389

390 # Copy the coordinates of the input cubes.

391 other_xr = xr.DataArray.from_iris(other)

392 base_xr = xr.DataArray.from_iris(base)

393 preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)

394

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

396 # apply scores, and then transform it back.

397 scores_cube = xr.DataArray.to_iris(

398 scores.continuous.correlation.pearsonr(

399 other_xr,

400 base_xr,

401 preserve_dims=preserve_dims,

402 )

403 )

404

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

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

407 try:

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

409 base_time = base.coord("time")

410 time_vals = (

411 base_time.bounds.flatten()

412 if base_time.has_bounds()

413 else base_time.points

414 )

415 t_start = float(time_vals[0])

416 t_end = float(time_vals[-1])

417 t_mid = 0.5 * (t_start + t_end)

418

419 scores_cube.add_aux_coord(

420 iris.coords.AuxCoord(

421 t_mid,

422 standard_name=base_time.standard_name,

423 long_name=base_time.long_name,

424 var_name=base_time.var_name,

425 units=base_time.units,

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

427 attributes=base_time.attributes.copy(),

428 )

429 )

430 except iris.exceptions.CoordinateNotFoundError:

431 pass

432

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

434 return scores_cube