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