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

1# © Crown copyright, Met Office (2022-2026) 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"""Operators to perform various kinds of image processing."""

17from typing import Literal

19import iris

20import iris.cube

21import numpy as np

22from skimage.measure import label

24from CSET._common import iter_maybe

25from CSET.operators.wind import calculate_vector_wind

28def MAUL_properties(

29 cubes: iris.cube.Cube | iris.cube.CubeList,

30 u_cubes: iris.cube.Cube | iris.cube.CubeList,

31 v_cubes: iris.cube.Cube | iris.cube.CubeList,

32 output: Literal["number", "base", "depth", "wind_below"],

33) -> iris.cube.Cube | iris.cube.CubeList:

34 """Identify properties of Moist Absolutely Unstable Layers.

36 Parameters

37 ----------

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

39 A cube or cubelist of a mask(s) as to whether a MAUL exists.

40 This input must be a binary field.

41 u_cubes: iris.cube.Cube | iris.cube.CubeList

42 A cube or cubelist of the wind in the u direction.

43 v_cubes: iris.cube.Cube | iris.cube.CubeList

44 A cube or cubelist of the wind in the v direction.

45 output: Literal["number", "base", "depth", "wind_below"]

46 The output is the desired property required. It can be

47 number, base, depth for the number of MAULs, base height

48 of the deepest MAUL, the depth of the deepest MAUL, or the

49 average windspeed below the MAUL, respectively.

52 Returns

53 -------

54 cube: iris.cube.Cube | iris.cube.CubeList

55 A Cube or CubeList depending upon the output specified.

57 Raises

58 ------

59 ValueError: Data contains values that are not 0 or 1, only masked data should be used.

60 This error is raised when a mask field is not provided to the operator.

61 ValueError: Unexpected value for output. Expected number, base, depth or wind_below. Got {output}.

62 This error is raised when the wrong output string is specified.

64 Notes

65 -----

66 Having been provided with a mask field for identifying whether Moist

67 Absolutely Unstable Layers (MAULs) are present, based on criteria

68 set out in a recipe. The operator applies image processing to the mask

69 to each point in the latitude/longitude coordinates. It uses the image

70 processing to identify continuous layers (1s), and labels them.

71 It identifies the number of layers by identifying the maximum label number,

72 and then finds the top and base of each layer. It will also find the average

73 windspeed below the MAUL for indications of presence of low-level jets.

74 Depending on the output desired it will output information for the deepest MAUL.

76 When a MAUL is not present the output will be set to NaN for depth and base.

77 If number of MAULs is the desired output it will be set to zero.

79 The MAUL diagnostic is applicable anywhere in the globe and across all scales.

80 The properties used here are based upon [Daviesetal24]_ and [Daviesetal26]_.

82 References

83 ----------

84 .. [Daviesetal24] Davies, P.A., Fowler, H.J, Villalobos-Herrera, R.,

85 Slingo, J., Flack, D.L.A., and Taszarek, M (2024)

86 "A New Conceptual Model for Understanding and Predicting Life-Threatening

87 Rainfall Extremes." Weather and Climate Extremes, vol. 45, 100696,

88 doi: 10.1016/j.wace.2024.100696

89 .. [Daviesetal26] Davies, P. A., Flack, D. L. A., Pirret, J., Fowler, H. J.

90 (2026) "Application of the Davies Four-Stage Conceptual Model for

91 Life-Threatening Rainfall Extremes on the April 2024 United Arab Emirates

92 and Oman Floods." Weather and Climate Extremes, vol. 51, 100846.

93 doi:10.1016/j.wace.2025.100846

94 """

95 num_MAULs = iris.cube.CubeList([])

96 maul_d = iris.cube.CubeList([])

97 maul_b = iris.cube.CubeList([])

98 windspeed_below_MAUL = iris.cube.CubeList([])

100 if output not in ("number", "base", "depth", "wind_below"):

101 raise ValueError(

102 f"""Unexpected value for output. Expected number, base, depth or wind_below. Got {output}."""

103 )

104

105 for cube, u, v in zip(

106 iter_maybe(cubes), iter_maybe(u_cubes), iter_maybe(v_cubes), strict=True

107 ):

108 # Check for binary fields.

109 if not np.array_equal(cube.data, cube.data.astype(bool)):

110 raise ValueError(

111 "Data contains values that are not 0 or 1, only masked data should be used."

112 )

113 # Create dummy cubes to store the output. The shape of the dummy cube

114 # depends upon which dimensions are present in the mask cube.

115 number_of_MAULs = next(cube.slices_over("model_level_number")).copy()

116 number_of_MAULs.data[:] = 0.0

117 maul_depth = number_of_MAULs.copy()

118 maul_base = number_of_MAULs.copy()

119 wind_below_maul = number_of_MAULs.copy()

120 # Calculate windspeed and direction.

121 windspeed_and_direction = calculate_vector_wind(u, v)

122 # Select windspeed, hard coded as always in same position from output

123 # of calculate_vector_wind.

124 windspeed = windspeed_and_direction[0]

125 # Loop over realization.

126 for mem_number, member in enumerate(cube.slices_over("realization")):

127 # Loop over time.

128 for time_point, time in enumerate(member.slices_over("time")):

129 # Loop over latitude.

130 for lat_point, lat in enumerate(time.slices_over("latitude")):

131 # Loop over longitude.

132 for lon_point, lon in enumerate(lat.slices_over("longitude")):

133 # Label each object in the vertical.

134 labels = label(lon.core_data())

135 # Finds the number of MAULs present based upon the

136 # number of objects identified, if no MAUL is present

137 # the value is set to zero.

138 # The code checks for whether there are multiple

139 # realization and/or time points for correct

140 # indexing of the output data and applies accordingly.

141 if (

142 len(number_of_MAULs.coord("realization").points) != 1

143 and len(number_of_MAULs.coord("time").points) != 1

144 ):

145 number_of_MAULs.data[

146 mem_number, time_point, lat_point, lon_point

147 ] = np.max(labels)

148 elif (

149 len(number_of_MAULs.coord("realization").points) != 1

150 and len(number_of_MAULs.coord("time").points) == 1

151 ):

152 number_of_MAULs.data[mem_number, lat_point, lon_point] = (

153 np.max(labels)

154 )

155 elif (

156 len(number_of_MAULs.coord("time").points) != 1

157 and len(number_of_MAULs.coord("realization").points) == 1

158 ):

159 number_of_MAULs.data[time_point, lat_point, lon_point] = (

160 np.max(labels)

161 )

162 else:

163 number_of_MAULs.data[lat_point, lon_point] = np.max(labels)

164 if output not in ("number", "wind_below"):

165 # Find the base, top, and depth for each object

166 # using cube metadata.

167 maul_start = []

168 maul_end = []

169 maul_dep = []

170 # Loop over the number of MAULs (plus one to ensure

171 # the case for only one MAUL being present).

172 for maul in range(1, np.max(labels) + 1):

173 # Find all vertical indices belonging to a MAUL.

174 maul_range = np.where(labels == maul)

175 # Find the height at the base of the MAUL

176 # (lowest level).

177 maul_start_point = lon.coord("level_height").points[

178 maul_range[0][0]

179 ]

180 # Find the height at the top of the MAUL

181 # (highest level).

182 maul_end_point = lon.coord("level_height").points[

183 maul_range[0][-1]

184 ]

185 # Calculate the MAUL depth, and store

186 # base and top heights.

187 maul_dep.append(maul_end_point - maul_start_point)

188 maul_start.append(maul_start_point)

189 maul_end.append(maul_end_point)

190 try:

191 # Idendtify where the deepest MAUL is.

192 index = int(

193 np.where(maul_dep == np.max(maul_dep))[0][0]

194 )

195 # As with number the code checks for whether

196 # there are multiple realization and/or time

197 # points for correct indexing of the output data

198 # and applies accordingly.

199 if (

200 len(number_of_MAULs.coord("realization").points)

201 != 1

202 and len(number_of_MAULs.coord("time").points) != 1

203 ):

204 # Store the deepest MAUL.

205 maul_depth.data[

206 mem_number, time_point, lat_point, lon_point

207 ] = np.max(maul_dep)

208 # Store the base height of the deepest MAUL.

209 maul_base.data[

210 mem_number, time_point, lat_point, lon_point

211 ] = maul_start[index]

212 elif (

213 len(number_of_MAULs.coord("realization").points)

214 != 1

215 and len(number_of_MAULs.coord("time").points) == 1

216 ):

217 maul_depth.data[

218 mem_number, lat_point, lon_point

219 ] = np.max(maul_dep)

220 maul_base.data[mem_number, lat_point, lon_point] = (

221 maul_start[index]

222 )

223 elif (

224 len(number_of_MAULs.coord("time").points) != 1

225 and len(number_of_MAULs.coord("realization").points)

226 == 1

227 ):

228 maul_depth.data[

229 time_point, lat_point, lon_point

230 ] = np.max(maul_dep)

231 maul_base.data[time_point, lat_point, lon_point] = (

232 maul_start[index]

233 )

234 else:

235 maul_depth.data[lat_point, lon_point] = np.max(

236 maul_dep

237 )

238 maul_base.data[lat_point, lon_point] = maul_start[

239 index

240 ]

241 # Here a ValueError is raised if a MAUL is not found, however

242 # this is a valid answer, and so output data is set to NaN.

243 # The dimensionality logic for output data is identical

244 # to that used previously.

245 except ValueError:

246 if (

247 len(number_of_MAULs.coord("realization").points)

248 != 1

249 and len(number_of_MAULs.coord("time").points) != 1

250 ):

251 maul_depth.data[

252 mem_number, time_point, lat_point, lon_point

253 ] = np.nan

254 maul_base.data[

255 mem_number, time_point, lat_point, lon_point

256 ] = np.nan

257 elif (

258 len(number_of_MAULs.coord("realization").points)

259 != 1

260 and len(number_of_MAULs.coord("time").points) == 1

261 ):

262 maul_depth.data[

263 mem_number, lat_point, lon_point

264 ] = np.nan

265 maul_base.data[mem_number, lat_point, lon_point] = (

266 np.nan

267 )

268 elif (

269 len(number_of_MAULs.coord("time").points) != 1

270 and len(number_of_MAULs.coord("realization").points)

271 == 1

272 ):

273 maul_depth.data[

274 time_point, lat_point, lon_point

275 ] = np.nan

276 maul_base.data[time_point, lat_point, lon_point] = (

277 np.nan

278 )

279 else:

280 maul_depth.data[lat_point, lon_point] = np.nan

281 maul_base.data[lat_point, lon_point] = np.nan

282 # Separate loop for calculating wind properties.

283 elif output not in ("number"):

284 # Find the base, top, and depth for each object

285 # using cube metadata.

286 maul_start = []

287 maul_end = []

288 maul_dep = []

289 # Loop over the number of MAULs. The loop starts

290 # at one as a value of zero implies there is not

291 # a MAUL present, so the first MAUL is one.

292 # Given this labelling convention plus one is required

293 # to ensure that the correct number of MAULs are

294 # looped over.

295 for maul in range(1, np.max(labels) + 1):

296 # Find all vertical indices belonging to a MAUL.

297 maul_range = np.where(labels == maul)

298 # Find the height at the base of the MAUL

299 # (lowest level).

300 maul_start_point = lon.coord("level_height").points[

301 maul_range[0][0]

302 ]

303 # Find the height at the top of the MAUL

304 # (highest level).

305 maul_end_point = lon.coord("level_height").points[

306 maul_range[0][-1]

307 ]

308 # Calculate the MAUL depth, and store

309 # base and top heights.

310 maul_dep.append(maul_end_point - maul_start_point)

311 maul_start.append(maul_start_point)

312 maul_end.append(maul_end_point)

313 try:

314 # Identify where the deepest MAUL is.

315 index = np.argmax(maul_dep)

316 maul_base_value = maul_start[index]

317 height_index = np.abs(

318 lon.coord("level_height").points - maul_base_value

319 ).argmin()

320 # As with number the code checks for whether

321 # there are multiple realization and/or time

322 # points for correct indexing of the output data

323 # and applies accordingly.

324 if (

325 len(number_of_MAULs.coord("realization").points)

326 != 1

327 and len(number_of_MAULs.coord("time").points) != 1

328 ):

329 # Store and calculate the windspeed below the

330 # deepest MAUL.

331 wind_below_maul.data[

332 mem_number, time_point, lat_point, lon_point

333 ] = np.mean(

334 windspeed[

335 mem_number,

336 time_point,

337 0:height_index,

338 lat_point,

339 lon_point,

340 ].data

341 )

342 elif (

343 len(number_of_MAULs.coord("realization").points)

344 != 1

345 and len(number_of_MAULs.coord("time").points) == 1

346 ):

347 wind_below_maul.data[

348 mem_number, lat_point, lon_point

349 ] = np.mean(

350 windspeed[

351 mem_number,

352 0:height_index,

353 lat_point,

354 lon_point,

355 ].data

356 )

357 elif (

358 len(number_of_MAULs.coord("time").points) != 1

359 and len(number_of_MAULs.coord("realization").points)

360 == 1

361 ):

362 wind_below_maul.data[

363 time_point, lat_point, lon_point

364 ] = np.mean(

365 windspeed[

366 time_point,

367 0:height_index,

368 lat_point,

369 lon_point,

370 ].data

371 )

372 else:

373 wind_below_maul.data[lat_point, lon_point] = (

374 np.mean(

375 windspeed[

376 0:height_index, lat_point, lon_point

377 ].data

378 )

379 )

380 # Here a ValueError is raised if a MAUL is not found, or an

381 # IndexError if the MAUL starts at the surface and so there

382 # is no wind below the MAUL however these are a valid answers,

383 # and so output data is set to NaN.

384 # The dimensionality logic for output data is identical

385 # to that used previously.

386 except (ValueError, IndexError):

387 if (

388 len(number_of_MAULs.coord("realization").points)

389 != 1

390 and len(number_of_MAULs.coord("time").points) != 1

391 ):

392 wind_below_maul.data[

393 mem_number, time_point, lat_point, lon_point

394 ] = np.nan

395 elif (

396 len(number_of_MAULs.coord("realization").points)

397 != 1

398 and len(number_of_MAULs.coord("time").points) == 1

399 ):

400 wind_below_maul.data[

401 mem_number, lat_point, lon_point

402 ] = np.nan

403 elif (

404 len(number_of_MAULs.coord("time").points) != 1

405 and len(number_of_MAULs.coord("realization").points)

406 == 1

407 ):

408 wind_below_maul.data[

409 time_point, lat_point, lon_point

410 ] = np.nan

411 else:

412 wind_below_maul.data[lat_point, lon_point] = np.nan

413

414 # Units and renaming for number, depth and base (the other case).

415 match output:

416 case "number":

417 number_of_MAULs.units = "1"

418 number_of_MAULs.rename("Number_of_MAULs")

419 num_MAULs.append(number_of_MAULs)

420 case "depth":

421 maul_depth.units = "m"

422 maul_depth.rename("MAUL_depth")

423 maul_d.append(maul_depth)

424 case "base":

425 maul_base.units = "m"

426 maul_base.rename("MAUL_base_height")

427 maul_b.append(maul_base)

428 case _:

429 wind_below_maul.units = "m s^-1"

430 wind_below_maul.rename("windspeed_below_MAUL")

431 windspeed_below_MAUL.append(wind_below_maul)

432

433 # Output data.

434 match output:

435 case "number" if len(num_MAULs) == 1:

436 return num_MAULs[0]

437 case "number":

438 return num_MAULs

439 case "depth" if len(maul_d) == 1:

440 return maul_d[0]

441 case "depth":

442 return maul_d

443 case "base" if len(maul_b) == 1:

444 return maul_b[0]

445 case "base":

446 return maul_b

447 case "wind_below" if len(windspeed_below_MAUL) == 1:

448 return windspeed_below_MAUL[0]

449 case _:

450 return windspeed_below_MAUL