Coverage for src/CSET/operators/misc.py: 100%

1# © Crown copyright, Met Office (2022-2025) and CSET contributors.

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

7# http://www.apache.org/licenses/LICENSE-2.0

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"""Miscellaneous operators."""

17import itertools

18import logging

19from collections.abc import Iterable

21import iris

22import numpy as np

23from iris.cube import Cube, CubeList

25from CSET._common import iter_maybe

26from CSET.operators._utils import fully_equalise_attributes, get_cube_yxcoordname

27from CSET.operators.regrid import regrid_onto_cube

30def noop(x, **kwargs):

31 """Return its input without doing anything to it.

33 Useful for constructing diagnostic chains.

35 Arguments

36 ---------

37 x: Any

38 Input to return.

40 Returns

41 -------

42 x: Any

43 The input that was given.

44 """

45 return x

48def remove_attribute(

49 cubes: Cube | CubeList, attribute: str | Iterable, **kwargs

50) -> CubeList:

51 """Remove a cube attribute.

53 If the attribute is not on the cube, the cube is passed through unchanged.

55 Arguments

56 ---------

57 cubes: Cube | CubeList

58 One or more cubes to remove the attribute from.

59 attribute: str | Iterable

60 Name of attribute (or Iterable of names) to remove.

62 Returns

63 -------

64 cubes: CubeList

65 CubeList of cube(s) with the attribute removed.

66 """

67 # Ensure cubes is a CubeList.

68 if not isinstance(cubes, CubeList):

69 cubes = CubeList(iter_maybe(cubes))

71 for cube in cubes:

72 for attr in iter_maybe(attribute):

73 cube.attributes.pop(attr, None)

74 return cubes

77def addition(addend_1, addend_2):

78 """Addition of two fields.

80 Parameters

81 ----------

82 addend_1: Cube

83 Any field to have another field added to it.

84 addend_2: Cube

85 Any field to be added to another field.

87 Returns

88 -------

89 Cube

91 Raises

92 ------

93 ValueError, iris.exceptions.NotYetImplementedError

94 When the cubes are not compatible.

96 Notes

97 -----

98 This is a simple operator designed for combination of diagnostics or

99 creating new diagnostics by using recipes.

100

101 Examples

102 --------

103 >>> field_addition = misc.addition(kinetic_energy_u, kinetic_energy_v)

104

105 """

106 return addend_1 + addend_2

107

108

109def subtraction(minuend, subtrahend):

110 """Subtraction of two fields.

111

112 Parameters

113 ----------

114 minuend: Cube

115 Any field to have another field subtracted from it.

116 subtrahend: Cube

117 Any field to be subtracted from to another field.

118

119 Returns

120 -------

121 Cube

122

123 Raises

124 ------

125 ValueError, iris.exceptions.NotYetImplementedError

126 When the cubes are not compatible.

127

128 Notes

129 -----

130 This is a simple operator designed for combination of diagnostics or

131 creating new diagnostics by using recipes. It can be used for model

132 differences to allow for comparisons between the same field in different

133 models or model configurations.

134

135 Examples

136 --------

137 >>> model_diff = misc.subtraction(temperature_model_A, temperature_model_B)

138

139 """

140 return minuend - subtrahend

141

142

143def division(numerator, denominator):

144 """Division of two fields.

145

146 Parameters

147 ----------

148 numerator: Cube

149 Any field to have the ratio taken with respect to another field.

150 denominator: Cube

151 Any field used to divide another field or provide the reference

152 value in a ratio.

153

154 Returns

155 -------

156 Cube

157

158 Raises

159 ------

160 ValueError

161 When the cubes are not compatible.

162

163 Notes

164 -----

165 This is a simple operator designed for combination of diagnostics or

166 creating new diagnostics by using recipes.

167

168 Examples

169 --------

170 >>> bowen_ratio = misc.division(sensible_heat_flux, latent_heat_flux)

171

172 """

173 return numerator / denominator

174

175

176def multiplication(multiplicand, multiplier):

177 """Multiplication of two fields.

178

179 Parameters

180 ----------

181 multiplicand: Cube

182 Any field to be multiplied by another field.

183 multiplier: Cube

184 Any field to be multiplied to another field.

185

186 Returns

187 -------

188 Cube

189

190 Raises

191 ------

192 ValueError

193 When the cubes are not compatible.

194

195 Notes

196 -----

197 This is a simple operator designed for combination of diagnostics or

198 creating new diagnostics by using recipes.

199

200 Examples

201 --------

202 >>> filtered_CAPE_ratio = misc.multiplication(CAPE_ratio, inflow_layer_properties)

203

204 """

205 return multiplicand * multiplier

206

207

208def combine_cubes_into_cubelist(first: Cube | CubeList, **kwargs) -> CubeList:

209 """Operator that combines multiple cubes or CubeLists into one.

210

211 Arguments

212 ---------

213 first: Cube | CubeList

214 First cube or CubeList to merge into CubeList.

215 second: Cube | CubeList

216 Second cube or CubeList to merge into CubeList. This must be a named

217 argument.

218 third: Cube | CubeList

219 There can be any number of additional arguments, they just need unique

220 names.

221 ...

222

223 Returns

224 -------

225 combined_cubelist: CubeList

226 Combined CubeList containing all cubes/CubeLists.

227

228 Raises

229 ------

230 TypeError:

231 If the provided arguments are not either a Cube or CubeList.

232 """

233 # Create empty CubeList to store cubes/CubeList.

234 all_cubes = CubeList()

235 # Combine all CubeLists into a single flat iterable.

236 for item in itertools.chain(iter_maybe(first), *map(iter_maybe, kwargs.values())):

237 # Check each item is a Cube, erroring if not.

238 if isinstance(item, Cube):

239 # Add cube to CubeList.

240 all_cubes.append(item)

241 else:

242 raise TypeError("Not a Cube or CubeList!", item)

243 return all_cubes

244

245

246def difference(cubes: CubeList):

247 """Difference of two fields.

248

249 Parameters

250 ----------

251 cubes: CubeList

252 A list of exactly two cubes. One must have the cset_comparison_base

253 attribute set to 1, and will be used as the base of the comparison.

254

255 Returns

256 -------

257 Cube

258

259 Raises

260 ------

261 ValueError

262 When the cubes are not compatible.

263

264 Notes

265 -----

266 This is a simple operator designed for combination of diagnostics or

267 creating new diagnostics by using recipes. It can be used for model

268 differences to allow for comparisons between the same field in different

269 models or model configurations.

270

271 Examples

272 --------

273 >>> model_diff = misc.difference(temperature_model_A, temperature_model_B)

274

275 """

276 if len(cubes) != 2:

277 raise ValueError("cubes should contain exactly 2 cubes.")

278 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))

279 other: Cube = cubes.extract_cube(

280 iris.Constraint(

281 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes

282 )

283 )

284

285 # Get spatial coord names.

286 base_lat_name, base_lon_name = get_cube_yxcoordname(base)

287 other_lat_name, other_lon_name = get_cube_yxcoordname(other)

288

289 # Ensure cubes to compare are on common differencing grid.

290 # This is triggered if either

291 # i) latitude and longitude shapes are not the same. Note grid points

292 # are not compared directly as these can differ through rounding

293 # errors.

294 # ii) or variables are known to often sit on different grid staggering

295 # in different models (e.g. cell center vs cell edge), as is the case

296 # for UM and LFRic comparisons.

297 # In future greater choice of regridding method might be applied depending

298 # on variable type. Linear regridding can in general be appropriate for smooth

299 # variables. Care should be taken with interpretation of differences

300 # given this dependency on regridding.

301 if (

302 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape

303 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape

304 ) or (

305 base.long_name

306 in [

307 "eastward_wind_at_10m",

308 "northward_wind_at_10m",

309 "northward_wind_at_cell_centres",

310 "eastward_wind_at_cell_centres",

311 "zonal_wind_at_pressure_levels",

312 "meridional_wind_at_pressure_levels",

313 "potential_vorticity_at_pressure_levels",

314 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",

315 ]

316 ):

317 logging.debug(

318 "Linear regridding base cube to other grid to compute differences"

319 )

320 base = regrid_onto_cube(base, other, method="Linear")

321

322 def is_increasing(sequence: list) -> bool:

323 """Determine the direction of an ordered sequence.

324

325 Returns "increasing" or "decreasing" depending on whether the sequence

326 is going up or down. The sequence should already be monotonic, with no

327 duplicate values. An iris DimCoord's points fulfills this criteria.

328 """

329 return sequence[0] < sequence[1]

330

331 # Figure out if we are comparing between UM and LFRic; flip array if so.

332 base_lat_direction = is_increasing(base.coord(base_lat_name).points)

333 other_lat_direction = is_increasing(other.coord(other_lat_name).points)

334 if base_lat_direction != other_lat_direction:

335 other.data = np.flip(other.data, other.coord(other_lat_name).cube_dims(other))

336

337 # Extract just common time points.

338 base, other = _extract_common_time_points(base, other)

339

340 # Equalise attributes so we can merge.

341 fully_equalise_attributes([base, other])

342 logging.debug("Base: %s\nOther: %s", base, other)

343

344 # This currently relies on the cubes having the same underlying data layout.

345 difference = base.copy()

346

347 # Differences don't have a standard name; long name gets a suffix. We are

348 # assuming we can rely on cubes having a long name, so we don't check for

349 # its presents.

350 difference.standard_name = None

351 if base.standard_name:

352 difference.long_name = base.standard_name + "_difference"

353 else:

354 difference.long_name = base.long_name + "_difference"

355 if base.var_name:

356 difference.var_name = base.var_name + "_difference"

357

358 difference.data = base.data - other.data

359 return difference

360

361

362def _extract_common_time_points(base: Cube, other: Cube) -> tuple[Cube, Cube]:

363 """Extract common time points from cubes to allow comparison."""

364 # Get the name of the first non-scalar time coordinate.

365 time_coord = next(

366 map(

367 lambda coord: coord.name(),

368 filter(

369 lambda coord: coord.shape > (1,) and coord.name() in ["time", "hour"],

370 base.coords(),

371 ),

372 ),

373 None,

374 )

375 if not time_coord:

376 logging.debug("No time coord, skipping equalisation.")

377 return (base, other)

378 base_time_coord = base.coord(time_coord)

379 other_time_coord = other.coord(time_coord)

380 logging.debug("Base: %s\nOther: %s", base_time_coord, other_time_coord)

381 if time_coord == "hour":

382 # We directly compare points when comparing coordinates with

383 # non-absolute units, such as hour. We can't just check the units are

384 # equal as iris automatically converts to datetime objects in the

385 # comparison for certain coordinate names.

386 base_times = base_time_coord.points

387 other_times = other_time_coord.points

388 shared_times = set.intersection(set(base_times), set(other_times))

389 else:

390 # Units don't match, so converting to datetimes for comparison.

391 base_times = base_time_coord.units.num2date(base_time_coord.points)

392 other_times = other_time_coord.units.num2date(other_time_coord.points)

393 shared_times = set.intersection(set(base_times), set(other_times))

394 logging.debug("Shared times: %s", shared_times)

395 time_constraint = iris.Constraint(

396 coord_values={time_coord: lambda cell: cell.point in shared_times}

397 )

398 # Extract points matching the shared times.

399 base = base.extract(time_constraint)

400 other = other.extract(time_constraint)

401 if base is None or other is None:

402 raise ValueError("No common time points found!")

403 return (base, other)

404

405

406def convert_units(cubes: iris.cube.Cube | iris.cube.CubeList, units: str):

407 """Convert the units of a cube.

408

409 Arguments

410 ---------

411 cubes: iris.cube.Cube | iris.cube.CubeList

412 A Cube or CubeList of a field for its units to be converted.

413

414 units: str

415 The unit that the original field is to be converted to. It takes

416 CF compliant units.

417

418 Returns

419 -------

420 iris.cube.Cube | iris.cube.CubeList

421 The field converted into the specified units.

422

423 Examples

424 --------

425 >>> T_in_F = misc.convert_units(temperature_in_K, "Fahrenheit")

426

427 """

428 new_cubelist = iris.cube.CubeList([])

429 for cube in iter_maybe(cubes):

430 # Copy cube to keep original data.

431 cube_a = cube.copy()

432 # Convert cube units.

433 cube_a.convert_units(units)

434 # Rename new cube so units in title to help with colorbar selection.

435 cube_a.rename(f"{cube.name()}_{units.replace(' ', '_')}")

436 new_cubelist.append(cube_a)

437 if len(new_cubelist) == 1:

438 return new_cubelist[0]

439 else:

440 return new_cubelist