Coverage for src/CSET/operators/precipitation.py: 99%
204 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-16 13:38 +0000
1# © Crown copyright, Met Office (2022-2026) 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."""
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.
83 The properties used here are based upon [Daviesetal24]_ and [Daviesetal26]_.
85 References
86 ----------
87 .. [Daviesetal24] Davies, P.A., Fowler, H.J, Villalobos-Herrera, R.,
88 Slingo, J., Flack, D.L.A., and Taszarek, M (2024)
89 "A New Conceptual Model for Understanding and Predicting Life-Threatening
90 Rainfall Extremes." Weather and Climate Extremes, vol. 45, 100696,
91 doi: 10.1016/j.wace.2024.100696
92 .. [Daviesetal26] Davies, P. A., Flack, D. L. A., Pirret, J., Fowler, H. J.
93 (2026) "Application of the Davies Four-Stage Conceptual Model for
94 Life-Threatening Rainfall Extremes on the April 2024 United Arab Emirates
95 and Oman Floods." Weather and Climate Extremes, vol. 51, 100846.
96 doi:10.1016/j.wace.2025.100846
98 Examples
99 --------
100 >>> No_MAULs = precipitation.MAUL_properties(maul_mask, u, v, output="number")
101 >>> MAUL_base = precipitation.MAUL_properties(maul_mask, u, v, output="base")
102 >>> MAUL_depth = precipitation.MAUL_properties(maul_mask, u, v, output="depth")
103 >>> Ave_windspeed_below_MAUL = precipitation.MAUL_properties(maul_mask, u, v, output="wind_below")
104 >>> Direction_shear_across_MAUL = precipitation.MAUL_properties(maul_mask, u, v, output="directional_shear")
105 """
106 num_MAULs = iris.cube.CubeList([])
107 maul_d = iris.cube.CubeList([])
108 maul_b = iris.cube.CubeList([])
109 windspeed_below_MAUL = iris.cube.CubeList([])
110 directional_shear_across_MAUL = iris.cube.CubeList([])
111 if output not in ("number", "base", "depth", "wind_below", "directional_shear"):
112 raise ValueError(
113 f"""Unexpected value for output. Expected number, base, depth, wind_below or directional_shear. Got {output}."""
114 )
116 for cube, u, v in zip(
117 iter_maybe(cubes), iter_maybe(u_cubes), iter_maybe(v_cubes), strict=True
118 ):
119 # Check for binary fields.
120 if not np.array_equal(cube.data, cube.data.astype(bool)):
121 raise ValueError(
122 "Data contains values that are not 0 or 1, only masked data should be used."
123 )
124 # Create dummy cubes to store the output. The shape of the dummy cube
125 # depends upon which dimensions are present in the mask cube.
126 number_of_MAULs = next(cube.slices_over("model_level_number")).copy()
127 number_of_MAULs.data[:] = 0.0
128 maul_depth = number_of_MAULs.copy()
129 maul_base = number_of_MAULs.copy()
130 wind_below_maul = number_of_MAULs.copy()
131 directional_shear = number_of_MAULs.copy()
132 # Calculate windspeed and direction.
133 windspeed_and_direction = calculate_vector_wind(u, v)
134 # Select windspeed, hard coded as always in same position from output
135 # of calculate_vector_wind.
136 windspeed = windspeed_and_direction[0]
137 # As above but for wind direction.
138 direction = windspeed_and_direction[1]
139 # Ensure direction in range +/- 180 for difference calculations.
140 direction.data[direction.data > 180.0] -= 360.0
141 # Loop over realization.
142 for mem_number, member in enumerate(cube.slices_over("realization")):
143 # Loop over time.
144 for time_point, time in enumerate(member.slices_over("time")):
145 # Loop over latitude.
146 for lat_point, lat in enumerate(time.slices_over("latitude")):
147 # Loop over longitude.
148 for lon_point, lon in enumerate(lat.slices_over("longitude")):
149 # Label each object in the vertical.
150 labels = label(lon.core_data())
151 # Finds the number of MAULs present based upon the
152 # number of objects identified, if no MAUL is present
153 # the value is set to zero.
154 # The code checks for whether there are multiple
155 # realization and/or time points for correct
156 # indexing of the output data and applies accordingly.
157 if (
158 len(number_of_MAULs.coord("realization").points) != 1
159 and len(number_of_MAULs.coord("time").points) != 1
160 ):
161 number_of_MAULs.data[
162 mem_number, time_point, lat_point, lon_point
163 ] = np.max(labels)
164 elif (
165 len(number_of_MAULs.coord("realization").points) != 1
166 and len(number_of_MAULs.coord("time").points) == 1
167 ):
168 number_of_MAULs.data[mem_number, lat_point, lon_point] = (
169 np.max(labels)
170 )
171 elif (
172 len(number_of_MAULs.coord("time").points) != 1
173 and len(number_of_MAULs.coord("realization").points) == 1
174 ):
175 number_of_MAULs.data[time_point, lat_point, lon_point] = (
176 np.max(labels)
177 )
178 else:
179 number_of_MAULs.data[lat_point, lon_point] = np.max(labels)
180 if output not in ("number", "wind_below", "directional_shear"):
181 # Find the base, top, and depth for each object
182 # using cube metadata.
183 maul_start = []
184 maul_end = []
185 maul_dep = []
186 # Loop over the number of MAULs (plus one to ensure
187 # the case for only one MAUL being present).
188 for maul in range(1, np.max(labels) + 1):
189 # Find all vertical indices belonging to a MAUL.
190 maul_range = np.where(labels == maul)
191 # Find the height at the base of the MAUL
192 # (lowest level).
193 maul_start_point = lon.coord("level_height").points[
194 maul_range[0][0]
195 ]
196 # Find the height at the top of the MAUL
197 # (highest level).
198 maul_end_point = lon.coord("level_height").points[
199 maul_range[0][-1]
200 ]
201 # Calculate the MAUL depth, and store
202 # base and top heights.
203 maul_dep.append(maul_end_point - maul_start_point)
204 maul_start.append(maul_start_point)
205 maul_end.append(maul_end_point)
206 try:
207 # Idendtify where the deepest MAUL is.
208 index = int(
209 np.where(maul_dep == np.max(maul_dep))[0][0]
210 )
211 # As with number the code checks for whether
212 # there are multiple realization and/or time
213 # points for correct indexing of the output data
214 # and applies accordingly.
215 if (
216 len(number_of_MAULs.coord("realization").points)
217 != 1
218 and len(number_of_MAULs.coord("time").points) != 1
219 ):
220 # Store the deepest MAUL.
221 maul_depth.data[
222 mem_number, time_point, lat_point, lon_point
223 ] = np.max(maul_dep)
224 # Store the base height of the deepest MAUL.
225 maul_base.data[
226 mem_number, time_point, lat_point, lon_point
227 ] = maul_start[index]
228 elif (
229 len(number_of_MAULs.coord("realization").points)
230 != 1
231 and len(number_of_MAULs.coord("time").points) == 1
232 ):
233 maul_depth.data[
234 mem_number, lat_point, lon_point
235 ] = np.max(maul_dep)
236 maul_base.data[mem_number, lat_point, lon_point] = (
237 maul_start[index]
238 )
239 elif (
240 len(number_of_MAULs.coord("time").points) != 1
241 and len(number_of_MAULs.coord("realization").points)
242 == 1
243 ):
244 maul_depth.data[
245 time_point, lat_point, lon_point
246 ] = np.max(maul_dep)
247 maul_base.data[time_point, lat_point, lon_point] = (
248 maul_start[index]
249 )
250 else:
251 maul_depth.data[lat_point, lon_point] = np.max(
252 maul_dep
253 )
254 maul_base.data[lat_point, lon_point] = maul_start[
255 index
256 ]
257 # Here a ValueError is raised if a MAUL is not found, however
258 # this is a valid answer, and so output data is set to NaN.
259 # The dimensionality logic for output data is identical
260 # to that used previously.
261 except ValueError:
262 if (
263 len(number_of_MAULs.coord("realization").points)
264 != 1
265 and len(number_of_MAULs.coord("time").points) != 1
266 ):
267 maul_depth.data[
268 mem_number, time_point, lat_point, lon_point
269 ] = np.nan
270 maul_base.data[
271 mem_number, time_point, lat_point, lon_point
272 ] = np.nan
273 elif (
274 len(number_of_MAULs.coord("realization").points)
275 != 1
276 and len(number_of_MAULs.coord("time").points) == 1
277 ):
278 maul_depth.data[
279 mem_number, lat_point, lon_point
280 ] = np.nan
281 maul_base.data[mem_number, lat_point, lon_point] = (
282 np.nan
283 )
284 elif (
285 len(number_of_MAULs.coord("time").points) != 1
286 and len(number_of_MAULs.coord("realization").points)
287 == 1
288 ):
289 maul_depth.data[
290 time_point, lat_point, lon_point
291 ] = np.nan
292 maul_base.data[time_point, lat_point, lon_point] = (
293 np.nan
294 )
295 else:
296 maul_depth.data[lat_point, lon_point] = np.nan
297 maul_base.data[lat_point, lon_point] = np.nan
298 # Separate loop for calculating wind properties.
299 elif output not in ("number"):
300 # Find the base, top, and depth for each object
301 # using cube metadata.
302 maul_start = []
303 maul_end = []
304 maul_dep = []
305 # Loop over the number of MAULs. The loop starts
306 # at one as a value of zero implies there is not
307 # a MAUL present, so the first MAUL is one.
308 # Given this labelling convention plus one is required
309 # to ensure that the correct number of MAULs are
310 # looped over.
311 for maul in range(1, np.max(labels) + 1):
312 # Find all vertical indices belonging to a MAUL.
313 maul_range = np.where(labels == maul)
314 # Find the height at the base of the MAUL
315 # (lowest level).
316 maul_start_point = lon.coord("level_height").points[
317 maul_range[0][0]
318 ]
319 # Find the height at the top of the MAUL
320 # (highest level).
321 maul_end_point = lon.coord("level_height").points[
322 maul_range[0][-1]
323 ]
324 # Calculate the MAUL depth, and store
325 # base and top heights.
326 maul_dep.append(maul_end_point - maul_start_point)
327 maul_start.append(maul_start_point)
328 maul_end.append(maul_end_point)
329 try:
330 # Identify where the deepest MAUL is.
331 index = np.argmax(maul_dep)
332 maul_base_value = maul_start[index]
333 maul_top_value = maul_end[index]
334 height_index = np.abs(
335 lon.coord("level_height").points - maul_base_value
336 ).argmin()
337 top_index = np.abs(
338 lon.coord("level_height").points - maul_top_value
339 ).argmin()
340 # As with number the code checks for whether
341 # there are multiple realization and/or time
342 # points for correct indexing of the output data
343 # and applies accordingly.
344 if (
345 len(number_of_MAULs.coord("realization").points)
346 != 1
347 and len(number_of_MAULs.coord("time").points) != 1
348 ):
349 # Store and calculate the windspeed below the
350 # deepest MAUL.
351 wind_below_maul.data[
352 mem_number, time_point, lat_point, lon_point
353 ] = np.mean(
354 windspeed[
355 mem_number,
356 time_point,
357 0:height_index,
358 lat_point,
359 lon_point,
360 ].data
361 )
362 # Store and calculate the directional wind shear (difference)
363 # across the deepest MAUL from the top to the bottom.
364 directional_shear.data[
365 mem_number, time_point, lat_point, lon_point
366 ] = (
367 direction[
368 mem_number,
369 time_point,
370 top_index,
371 lat_point,
372 lon_point,
373 ].data
374 - direction[
375 mem_number,
376 time_point,
377 height_index,
378 lat_point,
379 lon_point,
380 ].data
381 )
382 elif (
383 len(number_of_MAULs.coord("realization").points)
384 != 1
385 and len(number_of_MAULs.coord("time").points) == 1
386 ):
387 wind_below_maul.data[
388 mem_number, lat_point, lon_point
389 ] = np.mean(
390 windspeed[
391 mem_number,
392 0:height_index,
393 lat_point,
394 lon_point,
395 ].data
396 )
397 directional_shear.data[
398 mem_number, lat_point, lon_point
399 ] = (
400 direction[
401 mem_number, top_index, lat_point, lon_point
402 ].data
403 - direction[
404 mem_number,
405 height_index,
406 lat_point,
407 lon_point,
408 ].data
409 )
410 elif (
411 len(number_of_MAULs.coord("time").points) != 1
412 and len(number_of_MAULs.coord("realization").points)
413 == 1
414 ):
415 wind_below_maul.data[
416 time_point, lat_point, lon_point
417 ] = np.mean(
418 windspeed[
419 time_point,
420 0:height_index,
421 lat_point,
422 lon_point,
423 ].data
424 )
425 directional_shear.data[
426 time_point, lat_point, lon_point
427 ] = (
428 direction[
429 time_point, top_index, lat_point, lon_point
430 ].data
431 - direction[
432 time_point,
433 height_index,
434 lat_point,
435 lon_point,
436 ].data
437 )
438 else:
439 wind_below_maul.data[lat_point, lon_point] = (
440 np.mean(
441 windspeed[
442 0:height_index, lat_point, lon_point
443 ].data
444 )
445 )
446 directional_shear.data[lat_point, lon_point] = (
447 direction[top_index, lat_point, lon_point].data
448 - direction[
449 height_index, lat_point, lon_point
450 ].data
451 )
453 # Here a ValueError is raised if a MAUL is not found, or an
454 # IndexError if the MAUL starts at the surface and so there
455 # is no wind below the MAUL however these are a valid answers,
456 # and so output data is set to NaN.
457 # The dimensionality logic for output data is identical
458 # to that used previously.
459 except (ValueError, IndexError):
460 if (
461 len(number_of_MAULs.coord("realization").points)
462 != 1
463 and len(number_of_MAULs.coord("time").points) != 1
464 ):
465 wind_below_maul.data[
466 mem_number, time_point, lat_point, lon_point
467 ] = np.nan
468 directional_shear.data[
469 mem_number, time_point, lat_point, lon_point
470 ] = np.nan
471 elif (
472 len(number_of_MAULs.coord("realization").points)
473 != 1
474 and len(number_of_MAULs.coord("time").points) == 1
475 ):
476 wind_below_maul.data[
477 mem_number, lat_point, lon_point
478 ] = np.nan
479 directional_shear.data[
480 mem_number, lat_point, lon_point
481 ] = np.nan
482 elif (
483 len(number_of_MAULs.coord("time").points) != 1
484 and len(number_of_MAULs.coord("realization").points)
485 == 1
486 ):
487 wind_below_maul.data[
488 time_point, lat_point, lon_point
489 ] = np.nan
490 directional_shear.data[
491 time_point, lat_point, lon_point
492 ] = np.nan
493 else:
494 wind_below_maul.data[lat_point, lon_point] = np.nan
495 directional_shear.data[lat_point, lon_point] = (
496 np.nan
497 )
499 # Ensure directional shear differences are in range +/- 180 degrees.
500 directional_shear.data[directional_shear.data > 180.0] -= 360.0
501 directional_shear.data[directional_shear.data < -180.0] += 360.0
503 # Units and renaming for number, depth and base (the other case).
504 match output:
505 case "number":
506 number_of_MAULs.units = "1"
507 number_of_MAULs.rename("Number_of_MAULs")
508 num_MAULs.append(number_of_MAULs)
509 case "depth":
510 maul_depth.units = "m"
511 maul_depth.rename("MAUL_depth")
512 maul_d.append(maul_depth)
513 case "base":
514 maul_base.units = "m"
515 maul_base.rename("MAUL_base_height")
516 maul_b.append(maul_base)
517 case "wind_below":
518 wind_below_maul.units = "m s^-1"
519 wind_below_maul.rename("windspeed_below_MAUL")
520 windspeed_below_MAUL.append(wind_below_maul)
521 case _:
522 directional_shear.units = "degrees"
523 directional_shear.rename("directional_shear_across_MAUL")
524 directional_shear_across_MAUL.append(directional_shear)
526 # Output data.
527 match output:
528 case "number" if len(num_MAULs) == 1:
529 return num_MAULs[0]
530 case "number":
531 return num_MAULs
532 case "depth" if len(maul_d) == 1:
533 return maul_d[0]
534 case "depth":
535 return maul_d
536 case "base" if len(maul_b) == 1:
537 return maul_b[0]
538 case "base":
539 return maul_b
540 case "wind_below" if len(windspeed_below_MAUL) == 1:
541 return windspeed_below_MAUL[0]
542 case "wind_below": 542 ↛ 543line 542 didn't jump to line 543 because the pattern on line 542 never matched
543 return windspeed_below_MAUL
544 case "directional_shear" if len(directional_shear_across_MAUL) == 1:
545 return directional_shear_across_MAUL[0]
546 case _:
547 return directional_shear_across_MAUL
550def convert_rainfall_depth_to_rate(cubes, **kwargs):
551 """Convert rainfall depth to rate.
553 Convert rainfall depth (e.g. mm or kg m-2)
554 over a time interval into a rainfall rate (kg m-2 s-1).
556 The conversion uses the duration associated with the time coordinate:
557 - If time bounds are present, the bounds define the accumulation interval
558 - Otherwise, the interval is inferred from differences between time points
560 Arguments
561 ---------
562 cubes: iris.cube.Cube | iris.cube.CubeList
563 Cube(s) containing rainfall accumulation (depth), with units convertible
564 to kg m-2 (equivalent to mm of water).
566 Each cube must include a time coordinate, optionally with bounds.
568 kwargs:
569 Additional keyword arguments (currently unused, present for API compatibility).
571 Returns
572 -------
573 iris.cube.Cube | iris.cube.CubeList
574 Cube(s) with rainfall expressed as a rate in kg m-2 s-1.
576 The returned object matches the type of the input:
577 - single Cube → single Cube
578 - CubeList → CubeList
580 Raises
581 ------
582 ValueError
583 - If no time coordinate is present
584 - If only a single time point is available without bounds
585 - If any inferred duration is non-positive
587 Notes
588 -----
589 - Conversion relies on the equivalence:
590 1 mm of rainfall ≡ 1 kg m-2
591 - iris does not know that mm = kg m-2
593 - Unit handling:
594 * Cubes already in rate units (convertible to kg m-2 s-1) are left unchanged
595 * Cubes not representing accumulation (not convertible to kg m-2) are skipped
597 - Time handling:
598 * If units are time-reference units (e.g. "hours since ..."),
599 only the base unit (e.g. hours) is used for duration conversion
601 - Broadcasting:
602 The duration array is reshaped to match the time dimension of the cube
603 before division.
605 Examples
606 --------
607 >>> rate = precipitation.convert_rainfall_depth_to_rate(cube)
608 >>> rate_list = precipitation.convert_rainfall_depth_to_rate(cube_list)
609 """
610 from cf_units import Unit
612 cubes_list = iris.cube.CubeList(iter_maybe(cubes))
613 output = iris.cube.CubeList()
614 for cube in cubes_list:
615 # Identify input type
616 is_rate = cube.units.is_convertible("kg m-2 s-1") or cube.units.is_convertible(
617 "mm s-1"
618 )
619 is_mass_accum = cube.units.is_convertible("kg m-2")
620 is_depth_accum = cube.units.is_convertible("mm")
622 # Skip rates and unrelated variables
623 if is_rate or not (is_mass_accum or is_depth_accum):
624 output.append(cube)
625 continue
627 # Time coordinate is required for rainfall accumulations
628 try:
629 time_coord = cube.coord("time")
630 except iris.exceptions.CoordinateNotFoundError as exc:
631 raise ValueError("No time coordinate; cannot convert rainfall.") from exc
633 # Get accumulation duration
634 if time_coord.bounds is not None:
635 duration = time_coord.bounds[:, 1] - time_coord.bounds[:, 0]
636 else:
637 t = time_coord.points
639 if t.size < 2:
640 raise ValueError("Cannot infer duration from a single time point")
642 dt = np.diff(t)
643 duration = np.concatenate([dt, dt[-1:]]) # assume last interval repeats
645 # Convert duration to seconds
646 units = time_coord.units
647 if units.is_time_reference(): 647 ↛ 651line 647 didn't jump to line 651 because the condition on line 647 was always true
648 base_unit = str(units).split(" since ")[0].strip()
649 duration = Unit(base_unit).convert(duration, "seconds")
650 else:
651 duration = units.convert(duration, "seconds")
653 if np.any(duration <= 0):
654 raise ValueError("Non-positive rainfall accumulation interval detected.")
656 # Normalise rainfall accumulation units before dividing
657 # e.g. if rainfall amount is in cm, then the conversion will still work
659 data = cube.lazy_data()
661 if is_depth_accum:
662 factor = cube.units.convert(1.0, "mm")
663 data = data * factor
664 else:
665 factor = cube.units.convert(1.0, "kg m-2")
666 data = data * factor
668 # Reshape duration for broadcasting along time dimension
669 reshape = [1] * cube.ndim
670 time_dim = cube.coord_dims("time")[0]
671 reshape[time_dim] = -1
672 duration = duration.reshape(reshape)
674 # Convert depth(amount) to rate
675 # Numerically: mm s-1 == kg m-2 s-1
676 data = data / duration
677 new_cube = cube.copy(data=data)
678 new_cube.units = "kg m-2 s-1"
679 output.append(new_cube)
681 return output[0] if isinstance(cubes, iris.cube.Cube) else output