Coverage for src/CSET/operators/imageprocessing.py: 100%
57 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-07 15:12 +0000
1# © Crown copyright, Met Office (2022-2025) 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"""Operators to perform various kinds of image processing."""
17import logging
19import iris
20import iris.cube
21import numpy as np
22from skimage.metrics import structural_similarity
24from CSET._common import is_increasing
25from CSET.operators._utils import fully_equalise_attributes, get_cube_yxcoordname
26from CSET.operators.misc import _extract_common_time_points
27from CSET.operators.regrid import regrid_onto_cube
29logger = logging.getLogger(__name__)
32def _SSIM_cube_preparation(
33 cubes: iris.cube.CubeList,
34) -> (iris.cube.Cube, iris.cube.Cube, str):
35 """Prepare the cubes for the SSIM calculation.
37 Parameters
38 ----------
39 cubes: iris.cube.CubeList
40 A list of exactly two cubes. One must have the cset_comparison_base
41 attribute set to 1, and will be used as the base of the comparison.
42 The cubes must contain a time coordinate.
44 Returns
45 -------
46 iris.cube.Cube, iris.cube.Cube, str
48 Raises
49 ------
50 ValueError
51 When the cubes are not compatible, or no time coordinate.
52 """
53 if len(cubes) != 2:
54 raise ValueError("cubes should contain exactly 2 cubes.")
55 base: iris.cube.Cube = cubes.extract_cube(
56 iris.AttributeConstraint(cset_comparison_base=1)
57 )
58 other: iris.cube.Cube = cubes.extract_cube(
59 iris.Constraint(
60 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
61 )
62 )
64 # Get spatial coord names.
65 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
66 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
68 # Ensure cubes to compare are on common differencing grid.
69 # This is triggered if either
70 # i) latitude and longitude shapes are not the same. Note grid points
71 # are not compared directly as these can differ through rounding
72 # errors.
73 # ii) or variables are known to often sit on different grid staggering
74 # in different models (e.g. cell center vs cell edge), as is the case
75 # for UM and LFRic comparisons.
76 # In future greater choice of regridding method might be applied depending
77 # on variable type. Linear regridding can in general be appropriate for smooth
78 # variables. Care should be taken with interpretation of differences
79 # given this dependency on regridding.
80 if (
81 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
82 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
83 ) or (
84 base.long_name
85 in [
86 "eastward_wind_at_10m",
87 "northward_wind_at_10m",
88 "northward_wind_at_cell_centres",
89 "eastward_wind_at_cell_centres",
90 "zonal_wind_at_pressure_levels",
91 "meridional_wind_at_pressure_levels",
92 "potential_vorticity_at_pressure_levels",
93 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
94 ]
95 ):
96 logger.debug("Linear regridding base cube to other grid to compute differences")
97 base = regrid_onto_cube(base, other, method="Linear")
99 # Figure out if we are comparing between UM and LFRic; flip array if so.
100 base_lat_direction = is_increasing(base.coord(base_lat_name).points)
101 other_lat_direction = is_increasing(other.coord(other_lat_name).points)
102 if base_lat_direction != other_lat_direction:
103 other.data = np.flip(other.data, other.coord(other_lat_name).cube_dims(other))
105 # Extract just common time points.
106 base, other = _extract_common_time_points(base, other)
108 # Equalise attributes so we can merge.
109 fully_equalise_attributes([base, other])
110 logger.debug("Base: %s\nOther: %s", base, other)
112 # Get the name of the first non-scalar time coordinate.
113 time_coord = next(
114 (
115 coord.name()
116 for coord in filter(
117 lambda coord: coord.shape > (1,) and coord.name() in ["time", "hour"],
118 base.coords(),
119 )
120 ),
121 None,
122 )
124 if time_coord is None:
125 raise ValueError("Cubes should contain a time coordinate.")
126 # Create and empty CubeList for storing the time or realization data.
127 return base, other, time_coord
130def spatial_structural_similarity_model_comparisons(
131 cubes: iris.cube.CubeList, sigma: float = 1.5
132) -> iris.cube.Cube:
133 r"""Calculate the structural similarity and produces a spatial plot.
135 Parameters
136 ----------
137 cubes: iris.cube.CubeList
138 A list of exactly two cubes. One must have the cset_comparison_base
139 attribute set to 1, and will be used as the base of the comparison.
140 The cubes must contain a time coordinate.
141 sigma: float, optional
142 The standard deviation of the Gaussian kernel to be used. The default
143 is set to 1.5 to mimic the human eye following [Wangetal2004]_.
145 Returns
146 -------
147 iris.cube.Cube
149 Raises
150 ------
151 ValueError
152 When the cubes are not compatible, or no time coordinate.
154 Notes
155 -----
156 This diagnostic was introduced by Wang et al. (2004) [Wangetal2004]_. It is
157 an image processing diagnostic that takes into account three factors asscoiated
158 with an image: i) luminace, ii) contrast, iii) structure. In calculation terms
159 it is a combination of the intensity, variance, and co-variance of an image. It is
160 calculated as follows:
162 .. math:: SSIM(x,y) = \frac{(2\mu_{x}\mu_{y} + C_{1})(2\sigma_{xy} + C_{2})}{(\mu^{2}_{x}\mu^{2}_{y} + C_{1})(\sigma^{2}_{x}\sigma^{2}_{y} + C_{2})}
164 for images, x and y, and small constancts C1 and C2, with the other symbols having
165 their usual statistical meaning.
167 The diagnostic varies between positive and negative one, on the most part.
168 However, should the data being compared lie outside of the specified data
169 range values larger or smaller can occur. Values close to or exactly 1 imply
170 a perceptably similar image; values close to or exactly -1 imply an anticorrelated
171 image; small values imply the fields are perceptably different.
173 The diagnostic has been setup with default values that are designed to mimic the
174 human eye and so are universally applicable irrespective of model resolution.
175 However, it should be noted that the mean structural similarity is not
176 identical to the domain mean of the structural similarity. This difference
177 occurs because the former is calculated over the mean of all the windows
178 (Gaussian kernels) rather than the mean of the grid boxes.
180 Further details, including caveats, can be found in Wang et al. (2004)
181 [Wangetal2004]_.
183 Examples
184 --------
185 >>> SSIM = imageprocessing.spatial_structural_similarity_model_comparisons(
186 cubes, sigma=1.5)
187 >>> iplt.pcolormesh(SSIM[0,:], cmap=mpl.cm.bwr)
188 >>> plt.gca().coastlines('10m')
189 >>> plt.clim(-1, 1)
190 >>> plt.colorbar()
191 >>> plt.show()
192 """
193 base, other, time_coord = _SSIM_cube_preparation(cubes)
194 ssim = iris.cube.CubeList()
196 # Loop over realization and time coordinates.
197 for base_r, other_r in zip(
198 base.slices_over("realization"),
199 other.slices_over("realization"),
200 strict=True,
201 ):
202 for base_t, other_t in zip(
203 base_r.slices_over(time_coord), other_r.slices_over(time_coord), strict=True
204 ):
205 # Use the full array as output will be as a 2D map.
206 ssim_map = base_t.copy()
207 _, ssim_map.data = structural_similarity(
208 other_t.data,
209 base_t.data,
210 data_range=base_t.data.max() - base_t.data.min(),
211 gaussian_weights=True,
212 sigma=sigma,
213 full=True,
214 )
215 ssim.append(ssim_map)
216 # Merge the cube slices into one cube, rename, and change units.
217 ssim = ssim.merge_cube()
218 ssim.standard_name = None
219 ssim.long_name = "structural_similarity"
220 ssim.units = "1"
221 return ssim
224def mean_structural_similarity_model_comparisons(
225 cubes: iris.cube.CubeList, sigma: float = 1.5
226) -> iris.cube.Cube:
227 r"""Calculate the mean structural similarity and produces a timeseries.
229 Parameters
230 ----------
231 cubes: iris.cube.CubeList
232 A list of exactly two cubes. One must have the cset_comparison_base
233 attribute set to 1, and will be used as the base of the comparison.
234 The cubes must contain a time coordinate.
235 sigma: float, optional
236 The standard deviation of the Gaussian kernel to be used. The default
237 is set to 1.5 to mimic the human eye following [Wangetal2004a]_.
239 Returns
240 -------
241 iris.cube.Cube
243 Raises
244 ------
245 ValueError
246 When the cubes are not compatible, or no time coordinate.
248 Notes
249 -----
250 This diagnostic was introduced by Wang et al. (2004) [Wangetal2004a]_. It is
251 an image processing diagnostic that takes into account three factors asscoiated
252 with an image: i) luminace, ii) contrast, iii) structure. In calculation terms
253 it is a combination of the intensity, variance, and co-variance of an image. It is
254 calculated as follows:
256 .. math:: SSIM(x,y) = \frac{(2\mu_{x}\mu_{y} + C_{1})(2\sigma_{xy} + C_{2})}{(\mu^{2}_{x}\mu^{2}_{y} + C_{1})(\sigma^{2}_{x}\sigma^{2}_{y} + C_{2})}
258 for images, x and y, and small constancts C1 and C2, with the other symbols having
259 their usual statistical meaning.
261 The diagnostic varies between positive and negative one, on the most part.
262 However, should the data being compared lie outside of the specified data
263 range values larger or smaller can occur. Values close to or exactly 1 imply
264 a perceptably similar image; values close to or exactly -1 imply an anticorrelated
265 image; small values imply the fields are perceptably different.
267 The diagnostic has been setup with default values that are designed to mimic the
268 human eye and so are universally applicable irrespective of model resolution.
269 However, it should be noted that the mean structural similarity is not
270 identical to the domain mean of the structural similarity. This difference
271 occurs because the former is calculated over the mean of all the windows
272 (Gaussian kernels) rather than the mean of the grid boxes.
274 Further details, including caveats, can be found in Wang et al. (2004)
275 [Wangetal2004a]_.
277 Examples
278 --------
279 >>> MSSIM = imageprocessing.structural_similarity_model_comparisons(
280 cubes,sigma=1.5, spatial_plot=False)
281 >>> iplt.plot(MSSIM)
282 """
283 base, other, time_coord = _SSIM_cube_preparation(cubes)
284 ssim = iris.cube.CubeList()
286 # Loop over realization and time coordinates.
287 for base_r, other_r in zip(
288 base.slices_over("realization"),
289 other.slices_over("realization"),
290 strict=True,
291 ):
292 for base_t, other_t in zip(
293 base_r.slices_over(time_coord), other_r.slices_over(time_coord), strict=True
294 ):
295 # The MSSIM (Mean structural similarity) is compression to
296 # a single point. Therefore, copying cube data for one
297 # point in the domain to keep cube consistency.
298 mssim = base_t[0, 0].copy()
299 mssim.data = structural_similarity(
300 other_t.data,
301 base_t.data,
302 data_range=base_t.data.max() - base_t.data.min(),
303 gaussian_weights=True,
304 sigma=sigma,
305 )
306 ssim.append(mssim)
307 # Merge the cube slices into one cube, rename, and change units.
308 ssim = ssim.merge_cube()
309 ssim.standard_name = None
310 ssim.long_name = "structural_similarity"
311 ssim.units = "1"
312 return ssim