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", "directional_shear"],

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", "directional_shear"]

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, the

49 average windspeed below the MAUL, or the difference in

50 wind direction across the MAUL, respectively.

53 Returns

54 -------

55 cube: iris.cube.Cube | iris.cube.CubeList

56 A Cube or CubeList depending upon the output specified.

58 Raises

59 ------

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

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

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

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

65 Notes

66 -----

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

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

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

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

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

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

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

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

75 change in wind direction across the MAUL (top - base) is also calculated.

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

78 When a MAUL is not present the output will be set to NaN for depth, base, wind below

79 and directional shear. Should the MAUL start at the surface the wind below will also

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

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

84 Examples

85 --------

86 >>> No_MAULs = precipitation.MAUL_properties(maul_mask, u, v, output="number")

87 >>> MAUL_base = precipitation.MAUL_properties(maul_mask, u, v, output="base")

88 >>> MAUL_depth = precipitation.MAUL_properties(maul_mask, u, v, output="depth")

89 >>> Ave_windspeed_below_MAUL = precipitation.MAUL_properties(maul_mask, u, v, output="wind_below")

90 >>> Direction_shear_across_MAUL = precipitation.MAUL_properties(maul_mask, u, v, output="directional_shear")

91 """

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

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

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

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

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

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

98 raise ValueError(

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

100 )

101

102 for cube, u, v in zip(

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

104 ):

105 # Check for binary fields.

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

107 raise ValueError(

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

109 )

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

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

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

113 number_of_MAULs.data[:] = 0.0

114 maul_depth = number_of_MAULs.copy()

115 maul_base = number_of_MAULs.copy()

116 wind_below_maul = number_of_MAULs.copy()

117 directional_shear = number_of_MAULs.copy()

118 # Calculate windspeed and direction.

119 windspeed_and_direction = calculate_vector_wind(u, v)

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

121 # of calculate_vector_wind.

122 windspeed = windspeed_and_direction[0]

123 # As above but for wind direction.

124 direction = windspeed_and_direction[1]

125 # Ensure direction in range +/- 180 for difference calculations.

126 direction.data[direction.data > 180.0] -= 360.0

127 # Loop over realization.

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

129 # Loop over time.

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

131 # Loop over latitude.

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

133 # Loop over longitude.

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

135 # Label each object in the vertical.

136 labels = label(lon.core_data())

137 # Finds the number of MAULs present based upon the

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

139 # the value is set to zero.

140 # The code checks for whether there are multiple

141 # realization and/or time points for correct

142 # indexing of the output data and applies accordingly.

143 if (

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

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

146 ):

147 number_of_MAULs.data[

148 mem_number, time_point, lat_point, lon_point

149 ] = np.max(labels)

150 elif (

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

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

153 ):

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

155 np.max(labels)

156 )

157 elif (

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

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

160 ):

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

162 np.max(labels)

163 )

164 else:

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

166 if output not in ("number", "wind_below", "directional_shear"):

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

168 # using cube metadata.

169 maul_start = []

170 maul_end = []

171 maul_dep = []

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

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

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

175 # Find all vertical indices belonging to a MAUL.

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

177 # Find the height at the base of the MAUL

178 # (lowest level).

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

180 maul_range[0][0]

181 ]

182 # Find the height at the top of the MAUL

183 # (highest level).

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

185 maul_range[0][-1]

186 ]

187 # Calculate the MAUL depth, and store

188 # base and top heights.

189 maul_dep.append(maul_end_point - maul_start_point)

190 maul_start.append(maul_start_point)

191 maul_end.append(maul_end_point)

192 try:

193 # Idendtify where the deepest MAUL is.

194 index = int(

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

196 )

197 # As with number the code checks for whether

198 # there are multiple realization and/or time

199 # points for correct indexing of the output data

200 # and applies accordingly.

201 if (

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

203 != 1

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

205 ):

206 # Store the deepest MAUL.

207 maul_depth.data[

208 mem_number, time_point, lat_point, lon_point

209 ] = np.max(maul_dep)

210 # Store the base height of the deepest MAUL.

211 maul_base.data[

212 mem_number, time_point, lat_point, lon_point

213 ] = maul_start[index]

214 elif (

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

216 != 1

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

218 ):

219 maul_depth.data[

220 mem_number, lat_point, lon_point

221 ] = np.max(maul_dep)

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

223 maul_start[index]

224 )

225 elif (

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

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

228 == 1

229 ):

230 maul_depth.data[

231 time_point, lat_point, lon_point

232 ] = np.max(maul_dep)

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

234 maul_start[index]

235 )

236 else:

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

238 maul_dep

239 )

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

241 index

242 ]

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

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

245 # The dimensionality logic for output data is identical

246 # to that used previously.

247 except ValueError:

248 if (

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

250 != 1

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

252 ):

253 maul_depth.data[

254 mem_number, time_point, lat_point, lon_point

255 ] = np.nan

256 maul_base.data[

257 mem_number, time_point, lat_point, lon_point

258 ] = np.nan

259 elif (

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

261 != 1

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

263 ):

264 maul_depth.data[

265 mem_number, lat_point, lon_point

266 ] = np.nan

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

268 np.nan

269 )

270 elif (

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

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

273 == 1

274 ):

275 maul_depth.data[

276 time_point, lat_point, lon_point

277 ] = np.nan

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

279 np.nan

280 )

281 else:

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

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

284 # Separate loop for calculating wind properties.

285 elif output not in ("number"):

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

287 # using cube metadata.

288 maul_start = []

289 maul_end = []

290 maul_dep = []

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

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

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

294 # Find all vertical indices belonging to a MAUL.

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

296 # Find the height at the base of the MAUL

297 # (lowest level).

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

299 maul_range[0][0]

300 ]

301 # Find the height at the top of the MAUL

302 # (highest level).

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

304 maul_range[0][-1]

305 ]

306 # Calculate the MAUL depth, and store

307 # base and top heights.

308 maul_dep.append(maul_end_point - maul_start_point)

309 maul_start.append(maul_start_point)

310 maul_end.append(maul_end_point)

311 try:

312 # Identify where the deepest MAUL is.

313 index = int(

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

315 )

316 maul_base_value = maul_start[index]

317 maul_top_value = maul_end[index]

318 height_index = np.abs(

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

320 ).argmin()

321 top_index = np.abs(

322 lon.coord("level_height").points - maul_top_value

323 ).argmin()

324 # As with number the code checks for whether

325 # there are multiple realization and/or time

326 # points for correct indexing of the output data

327 # and applies accordingly.

328 if (

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

330 != 1

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

332 ):

333 # Store and calculate the windspeed below the

334 # deepest MAUL.

335 wind_below_maul.data[

336 mem_number, time_point, lat_point, lon_point

337 ] = np.mean(

338 windspeed[

339 mem_number,

340 time_point,

341 0:height_index,

342 lat_point,

343 lon_point,

344 ].data

345 )

346 # Store and calculate the directional wind shear (difference)

347 # across the deepest MAUL.

348 directional_shear.data[

349 mem_number, time_point, lat_point, lon_point

350 ] = (

351 direction[

352 mem_number,

353 time_point,

354 top_index,

355 lat_point,

356 lon_point,

357 ].data

358 - direction[

359 mem_number,

360 time_point,

361 height_index,

362 lat_point,

363 lon_point,

364 ].data

365 )

366 elif (

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

368 != 1

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

370 ):

371 wind_below_maul.data[

372 mem_number, lat_point, lon_point

373 ] = np.mean(

374 windspeed[

375 mem_number,

376 0:height_index,

377 lat_point,

378 lon_point,

379 ].data

380 )

381 directional_shear.data[

382 mem_number, lat_point, lon_point

383 ] = (

384 direction[

385 mem_number, top_index, lat_point, lon_point

386 ].data

387 - direction[

388 mem_number,

389 height_index,

390 lat_point,

391 lon_point,

392 ].data

393 )

394 elif (

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

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

397 == 1

398 ):

399 wind_below_maul.data[

400 time_point, lat_point, lon_point

401 ] = np.mean(

402 windspeed[

403 time_point,

404 0:height_index,

405 lat_point,

406 lon_point,

407 ].data

408 )

409 directional_shear.data[

410 time_point, lat_point, lon_point

411 ] = (

412 direction[

413 time_point, top_index, lat_point, lon_point

414 ].data

415 - direction[

416 time_point,

417 height_index,

418 lat_point,

419 lon_point,

420 ].data

421 )

422 else:

423 wind_below_maul.data[lat_point, lon_point] = (

424 np.mean(

425 windspeed[

426 0:height_index, lat_point, lon_point

427 ].data

428 )

429 )

430 directional_shear.data[lat_point, lon_point] = (

431 direction[top_index, lat_point, lon_point].data

432 - direction[

433 height_index, lat_point, lon_point

434 ].data

435 )

436

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

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

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

440 # and so output data is set to NaN.

441 # The dimensionality logic for output data is identical

442 # to that used previously.

443 except (ValueError, IndexError):

444 if (

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

446 != 1

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

448 ):

449 wind_below_maul.data[

450 mem_number, time_point, lat_point, lon_point

451 ] = np.nan

452 directional_shear.data[

453 mem_number, time_point, lat_point, lon_point

454 ] = np.nan

455 elif (

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

457 != 1

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

459 ):

460 wind_below_maul.data[

461 mem_number, lat_point, lon_point

462 ] = np.nan

463 directional_shear.data[

464 mem_number, lat_point, lon_point

465 ] = np.nan

466 elif (

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

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

469 == 1

470 ):

471 wind_below_maul.data[

472 time_point, lat_point, lon_point

473 ] = np.nan

474 directional_shear.data[

475 time_point, lat_point, lon_point

476 ] = np.nan

477 else:

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

479 directional_shear.data[lat_point, lon_point] = (

480 np.nan

481 )

482

483 # Ensure directional shear differences are in range +/- 180 degrees.

484 directional_shear.data[directional_shear.data > 180.0] -= 360.0

485 directional_shear.data[directional_shear.data < -180.0] += 360.0

486

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

488 match output:

489 case "number":

490 number_of_MAULs.units = "1"

491 number_of_MAULs.rename("Number_of_MAULs")

492 num_MAULs.append(number_of_MAULs)

493 case "depth":

494 maul_depth.units = "m"

495 maul_depth.rename("MAUL_depth")

496 maul_d.append(maul_depth)

497 case "base":

498 maul_base.units = "m"

499 maul_base.rename("MAUL_base_height")

500 maul_b.append(maul_base)

501 case "wind_below":

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

503 wind_below_maul.rename("windspeed_below_MAUL")

504 windspeed_below_MAUL.append(wind_below_maul)

505 case _:

506 directional_shear.units = "degrees"

507 directional_shear.rename("directional_shear_across_MAUL")

508 directional_shear_across_MAUL.append(directional_shear)

509

510 # Output data.

511 match output:

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

513 return num_MAULs[0]

514 case "number":

515 return num_MAULs

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

517 return maul_d[0]

518 case "depth":

519 return maul_d

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

521 return maul_b[0]

522 case "base":

523 return maul_b

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

525 return windspeed_below_MAUL[0]

526 case "wind_below":

527 return windspeed_below_MAUL

528 case "directional_shear" if len(directional_shear_across_MAUL) == 1:

529 return directional_shear_across_MAUL[0]

530 case _:

531 return directional_shear_across_MAUL