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