Coverage for src/CSET/operators/scoreswrappers.py: 90%
72 statements
« prev ^ index » next coverage.py v7.15.0, created at 2026-07-07 09: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 scores_rmse(cubes: CubeList, preserved_coordinates: list[str] | str | None = None):
150 r"""Calculate the Root Mean Square Error (RMSE) using scores.
152 Acts as a wrapper around the RMSE calculation from ``scores`` ([scoresa]_, [scoresb]_).
153 It is calculated as
155 .. math:: RMSE = \sqrt{\frac{1}{N} \Sigma(forecast - observations)^2}
157 Parameters
158 ----------
159 cubes: iris.cube.CubeList
160 A CubeList containing exactly two cubes: a base and an "other" model,
161 this can be an analysis and the model.
162 preserved_coordinates: list[str] | str | None, default is None.
163 The coordinates that you wish to preserve in the calculaiton of the
164 RMSE. For example if you want a map of each time you can preserve
165 ["time","grid_latitude", "grid_longitude"] or if you want a time series
166 you can preserve ["time"], if you want to collapse to a single value
167 use `None`. The default is `None`.
169 Returns
170 -------
171 RMSE: iris.cube.Cube
172 A cube containing the RMSE between the base and other cube.
174 References
175 ----------
176 .. [scoresa] Leeuwenburg, T., Loveday, N., Ebert, E. E., Cook, H.,
177 Khanarmuei, M., Taggart, R. J., Ramanathan, N., Carroll, M., Chong, S.,
178 Griffiths, A., & Sharples, J. (2024) "scores: A Python package for
179 verifying and evaluating models and predictions with xarray". Journal
180 of Open Source Software, vol. 9, 6889. doi: 10.21105/joss.06889
182 .. [scoresb] Leeuwenburg, T., Loveday, N., Ramanathan, N., Chong, S.,
183 Taggart, R. J., Shrestha, D., Khanarmuei, M., Cook, H., Bluett, L., Ebert,
184 E. E., Carroll, M., Trotta, B., Bishop, S., Squire, D. T., Griffiths, A.,
185 Pagano, T. C., Fisher, A. J., Mandelbaum, T., Jinghan, F., … Smallwood, J.
186 (2026) "scores: Metrics for the verification, evaluation and optimisation of
187 forecasts, predictions or models (2.5.0)". Zenodo. doi: 10.5281/zenodo.18638494
188 """
189 base, other = _sort_cubes_for_verification(cubes)
190 # Scores operators on xarray data arrays, so we transform the iris cube into an array,
191 # apply scores, and then transform it back.
192 RMSE = xr.DataArray.to_iris(
193 scores.continuous.rmse(
194 xr.DataArray.from_iris(other),
195 xr.DataArray.from_iris(base),
196 preserve_dims=preserved_coordinates,
197 )
198 )
199 RMSE.rename(f"RMSE_of_{base.name()}")
200 return RMSE
203def scores_crps_for_ensemble(
204 cubes: Cube | CubeList, method: str = "ecdf", control_member: int = 0
205) -> iris.Constraint:
206 r"""Calculate the CRPS for an ensemble.
208 Acts as a wrapper around the crps_for_ensemble from ``scores`` ([scores_a]_, [scores_b]_).
210 Lower CRPS values are better (implies experiment distribution is closer to control distribution/observations),
211 larger values are worse (implies distributions are dissimilar).
212 It is applicable across time and spatial scales as the focus is on the distribution of the values.
213 Default method is ecdf. ecdf is exact value from the empirical distributions,
214 whereas fair produces an approximated value based on a random sample of the underlying distribution.
216 See [CRPS] for further information.
218 Parameters
219 ----------
220 cubes: iris.cube.Cube
221 A Cube containing ensembles data
223 Returns
224 -------
225 crps: iris.cube.Cube
226 A cube containing the crps between the ensemble members and the control
228 References
229 ----------
230 .. [scores_a] Leeuwenburg, T., Loveday, N., Ebert, E. E., Cook, H.,
231 Khanarmuei, M., Taggart, R. J., Ramanathan, N., Carroll, M., Chong, S.,
232 Griffiths, A., & Sharples, J. (2024) "scores: A Python package for
233 verifying and evaluating models and predictions with xarray". Journal
234 of Open Source Software, vol. 9, 6889. doi: 10.21105/joss.06889
236 .. [scores_b] Leeuwenburg, T., Loveday, N., Ramanathan, N., Chong, S.,
237 Taggart, R. J., Shrestha, D., Khanarmuei, M., Cook, H., Bluett, L., Ebert,
238 E. E., Carroll, M., Trotta, B., Bishop, S., Squire, D. T., Griffiths, A.,
239 Pagano, T. C., Fisher, A. J., Mandelbaum, T., Jinghan, F., … Smallwood, J.
240 (2026) "scores: Metrics for the verification, evaluation and optimisation of
241 forecasts, predictions or models (2.5.0)". Zenodo. doi: 10.5281/zenodo.18638494
243 .. [CRPS]
244 Hersbach, H., 2000: Decomposition of the Continuous Ranked
245 Probability Score for Ensemble Prediction Systems. Wea.
246 Forecasting, 15, 559–570, https://doi.org/10.1175/1520-0434(2000)015<0559:DOTCRP>2.0.CO;2.
247 """
248 if control_member != 0:
249 logging.warning("control member is usual 0")
251 if control_member not in cubes.coords("realization")[0].points:
252 new_control_member = cubes.coords("realization")[0].points[0]
253 logging.warning(
254 f"control member value {control_member} out of bounds, defaulting to control member={new_control_member}"
255 )
256 control_member = new_control_member
258 if cubes.coord("time").shape[0] == 1:
259 raise ValueError("Cube has only one time coordinate.")
261 if cubes.coord("realization").shape[0] < 3:
262 raise ValueError("Cube should have one control member and at least two members")
264 ctrl = cubes.extract(generate_realization_constraint([control_member]))
265 ens_mem = cubes.extract(
266 generate_remove_single_ensemble_member_constraint(control_member)
267 )
269 # Realising the data in advance provides a large speedup
270 _ = ctrl.data
271 _ = ens_mem.data
272 del _
274 ctrl = xr.DataArray.from_iris(ctrl)
275 ens_mem = xr.DataArray.from_iris(ens_mem)
277 crps = xr.DataArray.to_iris(
278 scores.probability.crps_for_ensemble(
279 ens_mem,
280 ctrl,
281 ensemble_member_dim="realization",
282 method=method,
283 preserve_dims="time",
284 )
285 )
287 crps.rename(f"CRPS_of_{cubes[0].name()}")
288 _realization_callback(crps)
289 return crps