Coverage for src/CSET/operators/precipitation.py: 80%
176 statements
« prev ^ index » next coverage.py v7.15.1, created at 2026-07-15 14:27 +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"],
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 )
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"): 283 ↛ 286line 283 didn't jump to line 286 because the condition on line 283 was never true
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
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": 424 ↛ 428line 424 didn't jump to line 428 because the pattern on line 424 always matched
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)
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": 445 ↛ exitline 445 didn't return from function 'MAUL_properties' because the pattern on line 445 always matched
446 return maul_b
449def convert_rainfall_depth_to_rate(cubes, **kwargs):
450 """Convert rainfall depth to rate.
452 Convert rainfall depth (e.g. mm or kg m-2)
453 over a time interval into a rainfall rate (kg m-2 s-1).
455 The conversion uses the duration associated with the time coordinate:
456 - If time bounds are present, the bounds define the accumulation interval
457 - Otherwise, the interval is inferred from differences between time points
459 Arguments
460 ---------
461 cubes: iris.cube.Cube | iris.cube.CubeList
462 Cube(s) containing rainfall accumulation (depth), with units convertible
463 to kg m-2 (equivalent to mm of water).
465 Each cube must include a time coordinate, optionally with bounds.
467 kwargs:
468 Additional keyword arguments (currently unused, present for API compatibility).
470 Returns
471 -------
472 iris.cube.Cube | iris.cube.CubeList
473 Cube(s) with rainfall expressed as a rate in kg m-2 s-1.
475 The returned object matches the type of the input:
476 - single Cube → single Cube
477 - CubeList → CubeList
479 Raises
480 ------
481 ValueError
482 - If no time coordinate is present
483 - If only a single time point is available without bounds
484 - If any inferred duration is non-positive
486 Notes
487 -----
488 - Conversion relies on the equivalence:
489 1 mm of rainfall ≡ 1 kg m-2
490 - iris does not know that mm = kg m-2
492 - Unit handling:
493 * Cubes already in rate units (convertible to kg m-2 s-1) are left unchanged
494 * Cubes not representing accumulation (not convertible to kg m-2) are skipped
496 - Time handling:
497 * If units are time-reference units (e.g. "hours since ..."),
498 only the base unit (e.g. hours) is used for duration conversion
500 - Broadcasting:
501 The duration array is reshaped to match the time dimension of the cube
502 before division.
504 Examples
505 --------
506 >>> rate = precipitation.convert_rainfall_depth_to_rate(cube)
507 >>> rate_list = precipitation.convert_rainfall_depth_to_rate(cube_list)
508 """
509 from cf_units import Unit
511 cubes_list = iris.cube.CubeList(iter_maybe(cubes))
512 output = iris.cube.CubeList()
513 for cube in cubes_list:
514 # Identify input type
515 is_rate = cube.units.is_convertible("kg m-2 s-1") or cube.units.is_convertible(
516 "mm s-1"
517 )
518 is_mass_accum = cube.units.is_convertible("kg m-2")
519 is_depth_accum = cube.units.is_convertible("mm")
521 # Skip rates and unrelated variables
522 if is_rate or not (is_mass_accum or is_depth_accum):
523 output.append(cube)
524 continue
526 # Time coordinate is required for rainfall accumulations
527 try:
528 time_coord = cube.coord("time")
529 except iris.exceptions.CoordinateNotFoundError as exc:
530 raise ValueError("No time coordinate; cannot convert rainfall.") from exc
532 # Get accumulation duration
533 if time_coord.bounds is not None:
534 duration = time_coord.bounds[:, 1] - time_coord.bounds[:, 0]
535 else:
536 t = time_coord.points
538 if t.size < 2:
539 raise ValueError("Cannot infer duration from a single time point")
541 dt = np.diff(t)
542 duration = np.concatenate([dt, dt[-1:]]) # assume last interval repeats
544 # Convert duration to seconds
545 units = time_coord.units
546 if units.is_time_reference(): 546 ↛ 550line 546 didn't jump to line 550 because the condition on line 546 was always true
547 base_unit = str(units).split(" since ")[0].strip()
548 duration = Unit(base_unit).convert(duration, "seconds")
549 else:
550 duration = units.convert(duration, "seconds")
552 if np.any(duration <= 0):
553 raise ValueError("Non-positive rainfall accumulation interval detected.")
555 # Normalise rainfall accumulation units before dividing
556 # e.g. if rainfall amount is in cm, then the conversion will still work
558 data = cube.lazy_data()
560 if is_depth_accum:
561 factor = cube.units.convert(1.0, "mm")
562 data = data * factor
563 else:
564 factor = cube.units.convert(1.0, "kg m-2")
565 data = data * factor
567 # Reshape duration for broadcasting along time dimension
568 reshape = [1] * cube.ndim
569 time_dim = cube.coord_dims("time")[0]
570 reshape[time_dim] = -1
571 duration = duration.reshape(reshape)
573 # Convert depth(amount) to rate
574 # Numerically: mm s-1 == kg m-2 s-1
575 data = data / duration
576 new_cube = cube.copy(data=data)
577 new_cube.units = "kg m-2 s-1"
578 output.append(new_cube)
580 return output[0] if isinstance(cubes, iris.cube.Cube) else output