Coverage for src/CSET/operators/scoreswrappers.py: 85%
47 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-17 15:44 +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 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")
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
136 # Equalise attributes so we can merge.
137 fully_equalise_attributes(CubeList([base, other]))
138 logging.debug("Base: %s\nOther: %s", base, other)
140 return base, other
143def scores_rmse(cubes: CubeList, preserved_coordinates: list[str] | str | None = None):
144 r"""Calculate the Root Mean Square Error (RMSE) using scores.
146 Acts as a wrapper around the RMSE calculation from ``scores`` ([scoresa]_, [scoresb]_).
147 It is calculated as
149 .. math:: RMSE = \sqrt{\frac{1}{N} \Sigma(forecast - observations)^2}
151 Parameters
152 ----------
153 cubes: iris.cube.CubeList
154 A CubeList containing exactly two cubes: a base and an "other" model,
155 this can be an analysis and the model.
156 preserved_coordinates: list[str] | str | None, default is None.
157 The coordinates that you wish to preserve in the calculaiton of the
158 RMSE. For example if you want a map of each time you can preserve
159 ["time","grid_latitude", "grid_longitude"] or if you want a time series
160 you can preserve ["time"], if you want to collapse to a single value
161 use `None`. The default is `None`.
163 Returns
164 -------
165 RMSE: iris.cube.Cube
166 A cube containing the RMSE between the base and other cube.
168 References
169 ----------
170 .. [scoresa] Leeuwenburg, T., Loveday, N., Ebert, E. E., Cook, H.,
171 Khanarmuei, M., Taggart, R. J., Ramanathan, N., Carroll, M., Chong, S.,
172 Griffiths, A., & Sharples, J. (2024) "scores: A Python package for
173 verifying and evaluating models and predictions with xarray". Journal
174 of Open Source Software, vol. 9, 6889. doi: 10.21105/joss.06889
176 .. [scoresb] Leeuwenburg, T., Loveday, N., Ramanathan, N., Chong, S.,
177 Taggart, R. J., Shrestha, D., Khanarmuei, M., Cook, H., Bluett, L., Ebert,
178 E. E., Carroll, M., Trotta, B., Bishop, S., Squire, D. T., Griffiths, A.,
179 Pagano, T. C., Fisher, A. J., Mandelbaum, T., Jinghan, F., … Smallwood, J.
180 (2026) "scores: Metrics for the verification, evaluation and optimisation of
181 forecasts, predictions or models (2.5.0)". Zenodo. doi: 10.5281/zenodo.18638494
182 """
183 base, other = _sort_cubes_for_verification(cubes)
184 # Scores operators on xarray data arrays, so we transform the iris cube into an array,
185 # apply scores, and then transform it back.
186 RMSE = xr.DataArray.to_iris(
187 scores.continuous.rmse(
188 xr.DataArray.from_iris(other),
189 xr.DataArray.from_iris(base),
190 preserve_dims=preserved_coordinates,
191 )
192 )
193 RMSE.rename(f"RMSE_of_{base.name()}")
194 return RMSE