Coverage for src/CSET/operators/scoreswrappers.py: 90%
40 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-16 13:58 +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 scores
21import xarray as xr
22from iris.cube import Cube
23from iris.util import reverse
25from CSET._common import is_increasing
26from CSET.operators._utils import fully_equalise_attributes, get_cube_yxcoordname
27from CSET.operators.misc import _extract_common_time_points
28from CSET.operators.regrid import regrid_onto_cube
31def _sort_cubes_for_verification(cubes: iris.cube.CubeList):
32 """Prepare cubes ready for verification in scores.
34 Parameters
35 ----------
36 cubes: iris.cube.CubeList
37 A CubeList of exact 2 cubes, one from each model.
39 Returns
40 -------
41 base: iris.cube.Cube
42 The cube from the "analysis" in the same format as the other model.
43 other: iris.cube.Cube
44 The cube from the model in the same format as the base model.
46 Raises
47 ------
48 ValueError: "cubes should contain exactly 2 cubes."
49 If any other number of cubes are present.
51 Notes
52 -----
53 This operator is used for sorting the data into the correct format. It
54 is likely going to need to be refactored out of CSET and perhaps moved into
55 `CSET._utils` given common code between here and `misc.difference`.
56 """
57 # Set cubes into correct format using code from difference operator
58 if len(cubes) != 2:
59 raise ValueError("cubes should contain exactly 2 cubes.")
60 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
61 other: Cube = cubes.extract_cube(
62 iris.Constraint(
63 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
64 )
65 )
67 # If cubes contain a pressure coordinate, ensure it is increasing.
68 for cube in cubes:
69 try:
70 if len(cube.coord("pressure").points) > 2: 70 ↛ 68line 70 didn't jump to line 68 because the condition on line 70 was always true
71 if not is_increasing(cube.coord("pressure").points): 71 ↛ 72line 71 didn't jump to line 72 because the condition on line 71 was never true
72 reverse(cube, "pressure")
74 except iris.exceptions.CoordinateNotFoundError:
75 pass
77 # Get spatial coord names.
78 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
79 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
81 # Ensure cubes to compare are on common differencing grid.
82 # This is triggered if either
83 # i) latitude and longitude shapes are not the same. Note grid points
84 # are not compared directly as these can differ through rounding
85 # errors.
86 # ii) or variables are known to often sit on different grid staggering
87 # in different models (e.g. cell center vs cell edge), as is the case
88 # for UM and LFRic comparisons.
89 # In future greater choice of regridding method might be applied depending
90 # on variable type. Linear regridding can in general be appropriate for smooth
91 # variables. Care should be taken with interpretation of differences
92 # given this dependency on regridding.
93 if (
94 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
95 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
96 ) or (
97 base.long_name
98 in [
99 "eastward_wind_at_10m",
100 "northward_wind_at_10m",
101 "northward_wind_at_cell_centres",
102 "eastward_wind_at_cell_centres",
103 "zonal_wind_at_pressure_levels",
104 "meridional_wind_at_pressure_levels",
105 "potential_vorticity_at_pressure_levels",
106 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
107 ]
108 ):
109 logging.debug(
110 "Linear regridding base cube to other grid to compute differences"
111 )
112 base = regrid_onto_cube(base, other, method="Linear")
114 # Figure out if we are comparing between UM and LFRic; flip array if so.
115 base_lat_direction = is_increasing(base.coord(base_lat_name).points)
116 other_lat_direction = is_increasing(other.coord(other_lat_name).points)
117 if base_lat_direction != other_lat_direction: 117 ↛ 118line 117 didn't jump to line 118 because the condition on line 117 was never true
118 reverse(other, other_lat_name)
120 # Extract just common time points.
121 base, other = _extract_common_time_points(base, other)
123 # Equalise attributes so we can merge.
124 fully_equalise_attributes([base, other])
125 logging.debug("Base: %s\nOther: %s", base, other)
127 return base, other
130def scores_rmse(
131 cubes: iris.cube.CubeList, preserved_coordinates: list[str] | str | None = None
132):
133 r"""Calculate the Root Mean Square Error (RMSE) using scores.
135 Acts as a wrapper around the RMSE calculation from `scores` ([scoresa]_, [scoresb]_).
136 It is calculated as
138 .. math:: RMSE = \sqrt{\frac{1}{N} \Sigma(forecast - observations)^2}
140 Parameters
141 ----------
142 cubes: iris.cube.CubeList
143 A CubeList containing exactly two cubes: a base and an "other" model,
144 this can be an analysis and the model.
145 preserved_coordinates: list[str] | str | None, default is None.
146 The coordinates that you wish to preserve in the calculaiton of the
147 RMSE. For example if you want a map of each time you can preserve
148 ["time","grid_latitude", "grid_longitude"] or if you want a time series
149 you can preserve ["time"], if you want to collapse to a single value
150 use `None`. The default is `None`.
152 Returns
153 -------
154 RMSE: iris.cube.Cube
155 A cube containing the RMSE between the base and other cube.
157 References
158 ----------
159 .. [scoresa] Leeuwenburg, T., Loveday, N., Ebert, E. E., Cook, H.,
160 Khanarmuei, M., Taggart, R. J., Ramanathan, N., Carroll, M., Chong, S.,
161 Griffiths, A., & Sharples, J. (2024) "scores: A Python package for
162 verifying and evaluating models and predictions with xarray". Journal
163 of Open Source Software, vol. 9, 6889. doi: 10.21105/joss.06889
165 .. [scoresb] Leeuwenburg, T., Loveday, N., Ramanathan, N., Chong, S.,
166 Taggart, R. J., Shrestha, D., Khanarmuei, M., Cook, H., Bluett, L., Ebert,
167 E. E., Carroll, M., Trotta, B., Bishop, S., Squire, D. T., Griffiths, A.,
168 Pagano, T. C., Fisher, A. J., Mandelbaum, T., Jinghan, F., … Smallwood, J.
169 (2026) "scores: Metrics for the verification, evaluation and optimisation of
170 forecasts, predictions or models (2.5.0)". Zenodo. doi: 10.5281/zenodo.18638494
172 """
173 base, other = _sort_cubes_for_verification(cubes)
174 # Scores operators on xarray data arrays, so we transform the iris cube into an array,
175 # apply scores, and then transform it back.
176 RMSE = xr.DataArray.to_iris(
177 scores.continuous.rmse(
178 xr.DataArray.from_iris(other),
179 xr.DataArray.from_iris(base),
180 preserve_dims=preserved_coordinates,
181 )
182 )
183 RMSE.rename(f"RMSE_of_{base.name()}")
184 return RMSE