Coverage for src/CSET/operators/precipitation.py: 100%
93 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-09 07:49 +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
27def MAUL_properties(
28 cubes: iris.cube.Cube | iris.cube.CubeList,
29 output: Literal["number", "base", "depth"],
30) -> iris.cube.Cube | iris.cube.CubeList:
31 """Identify properties of Moist Absolutely Unstable Layers.
33 Parameters
34 ----------
35 cubes: iris.cube.Cube | iris.cube.CubeList
36 A cube or cubelist of a mask(s) as to whether a MAUL exists.
37 This input must be a binary field.
38 output: Literal["number", "base", "depth"]
39 The output is the desired property required. It can be
40 number, base, depth for the number of MAULs, base height
41 of the deepest MAUL, or the depth of the deepest MAUL,
42 respectively.
44 Returns
45 -------
46 cube: iris.cube.Cube | iris.cube.CubeList
47 A Cube or CubeList depending upon the output specified.
49 Raises
50 ------
51 ValueError: Data contains values that are not 0 or 1, only masked data should be used.
52 This error is raised when a mask field is not provided to the operator.
53 ValueError: Unexpected value for output. Expected number, depth or base. Got {output}.
54 This error is raised when the wrong output string is specified.
56 Notes
57 -----
58 Having been provided with a mask field for identifying whether Moist
59 Absolutely Unstable Layers (MAULs) are present, based on criteria
60 set out in a recipe. The operator applies image processing to the mask
61 to each point in the latitude/longitude coordinates. It uses the image
62 processing to identify continuous layers (1s), and labels them.
63 It identifies the number of layesr by identifying the maximum label number,
64 and then finds the top and base of each layer. Depending on the output
65 desired it will output information for the deepest MAUL.
67 When a MAUL is not present the output will be set to NaN for depth and base.
68 If number of MAULs is the desired output it will be set to zero.
70 The MAUL diagnostic is applicable anywhere in the globe and across all scales.
71 The properties used here are based upon [Daviesetal24]_ and [Daviesetal26]_.
73 References
74 ----------
75 .. [Daviesetal24] Davies, P.A., Fowler, H.J, Villalobos-Herrera, R.,
76 Slingo, J., Flack, D.L.A., and Taszarek, M (2024)
77 "A New Conceptual Model for Understanding and Predicting Life-Threatening
78 Rainfall Extremes." Weather and Climate Extremes, vol. 45, 100696,
79 doi: 10.1016/j.wace.2024.100696
80 .. [Daviesetal26] Davies, P. A., Flack, D. L. A., Pirret, J., Fowler, H. J.
81 (2026) "Application of the Davies Four-Stage Conceptual Model for
82 Life-Threatening Rainfall Extremes on the April 2024 United Arab Emirates
83 and Oman Floods." Weather and Climate Extremes, vol. 51, 100846.
84 doi:10.1016/j.wace.2025.100846
85 """
86 num_MAULs = iris.cube.CubeList([])
87 maul_d = iris.cube.CubeList([])
88 maul_b = iris.cube.CubeList([])
90 if output not in ("number", "base", "depth"):
91 raise ValueError(
92 f"""Unexpected value for output. Expected number, depth or base. Got {output}."""
93 )
95 for cube in iter_maybe(cubes):
96 # Check for binary fields.
97 if not np.array_equal(cube.data, cube.data.astype(bool)):
98 raise ValueError(
99 "Data contains values that are not 0 or 1, only masked data should be used."
100 )
101 # Create dummy cubes to store the output. The shape of the dummy cube
102 # depends upon which dimensions are present in the mask cube.
103 number_of_MAULs = next(cube.slices_over("model_level_number")).copy()
104 number_of_MAULs.data[:] = 0.0
105 maul_depth = number_of_MAULs.copy()
106 maul_base = number_of_MAULs.copy()
107 # Loop over realization.
108 for mem_number, member in enumerate(cube.slices_over("realization")):
109 # Loop over time.
110 for time_point, time in enumerate(member.slices_over("time")):
111 # Loop over latitude.
112 for lat_point, lat in enumerate(time.slices_over("latitude")):
113 # Loop over longitude.
114 for lon_point, lon in enumerate(lat.slices_over("longitude")):
115 # Label each object in the vertical.
116 labels = label(lon.core_data())
117 # Finds the number of MAULs present based upon the
118 # number of objects identified, if no MAUL is present
119 # the value is set to zero.
120 # The code checks for whether there are multiple
121 # realization and/or time points for correct
122 # indexing of the output data and applies accordingly.
123 if (
124 len(number_of_MAULs.coord("realization").points) != 1
125 and len(number_of_MAULs.coord("time").points) != 1
126 ):
127 number_of_MAULs.data[
128 mem_number, time_point, lat_point, lon_point
129 ] = np.max(labels)
130 elif (
131 len(number_of_MAULs.coord("realization").points) != 1
132 and len(number_of_MAULs.coord("time").points) == 1
133 ):
134 number_of_MAULs.data[mem_number, lat_point, lon_point] = (
135 np.max(labels)
136 )
137 elif (
138 len(number_of_MAULs.coord("time").points) != 1
139 and len(number_of_MAULs.coord("realization").points) == 1
140 ):
141 number_of_MAULs.data[time_point, lat_point, lon_point] = (
142 np.max(labels)
143 )
144 else:
145 number_of_MAULs.data[lat_point, lon_point] = np.max(labels)
146 if output != "number":
147 # Find the base, top, and depth for each object
148 # using cube metadata.
149 maul_start = []
150 maul_end = []
151 maul_dep = []
152 # Loop over the number of MAULs (plus one to ensure
153 # the case for only one MAUL being present).
154 for maul in range(1, np.max(labels) + 1):
155 # Find all vertical indices belonging to a MAUL.
156 maul_range = np.where(labels == maul)
157 # Find the height at the base of the MAUL
158 # (lowest level).
159 maul_start_point = lon.coord("level_height").points[
160 maul_range[0][0]
161 ]
162 # Find the height at the top of the MAUL
163 # (highest level).
164 maul_end_point = lon.coord("level_height").points[
165 maul_range[0][-1]
166 ]
167 # Calculate the MAUL depth, and store
168 # base and top heights.
169 maul_dep.append(maul_end_point - maul_start_point)
170 maul_start.append(maul_start_point)
171 maul_end.append(maul_end_point)
172 try:
173 # Idendtify where the deepest MAUL is.
174 index = int(
175 np.where(maul_dep == np.max(maul_dep))[0][0]
176 )
177 # As with number the code checks for whether
178 # there are multiple realization and/or time
179 # points for correct indexing of the output data
180 # and applies accordingly.
181 if (
182 len(number_of_MAULs.coord("realization").points)
183 != 1
184 and len(number_of_MAULs.coord("time").points) != 1
185 ):
186 # Store the deepest MAUL.
187 maul_depth.data[
188 mem_number, time_point, lat_point, lon_point
189 ] = np.max(maul_dep)
190 # Store the base height of the deepest MAUL.
191 maul_base.data[
192 mem_number, time_point, lat_point, lon_point
193 ] = maul_start[index]
194 elif (
195 len(number_of_MAULs.coord("realization").points)
196 != 1
197 and len(number_of_MAULs.coord("time").points) == 1
198 ):
199 maul_depth.data[
200 mem_number, lat_point, lon_point
201 ] = np.max(maul_dep)
202 maul_base.data[mem_number, lat_point, lon_point] = (
203 maul_start[index]
204 )
205 elif (
206 len(number_of_MAULs.coord("time").points) != 1
207 and len(number_of_MAULs.coord("realization").points)
208 == 1
209 ):
210 maul_depth.data[
211 time_point, lat_point, lon_point
212 ] = np.max(maul_dep)
213 maul_base.data[time_point, lat_point, lon_point] = (
214 maul_start[index]
215 )
216 else:
217 maul_depth.data[lat_point, lon_point] = np.max(
218 maul_dep
219 )
220 maul_base.data[lat_point, lon_point] = maul_start[
221 index
222 ]
223 # Here a ValueError is raised if a MAUL is not found, however
224 # this is a valid answer, and so output data is set to NaN.
225 # The dimensionality logic for output data is identical
226 # to that used previously.
227 except ValueError:
228 if (
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, time_point, lat_point, lon_point
235 ] = np.nan
236 maul_base.data[
237 mem_number, time_point, lat_point, lon_point
238 ] = np.nan
239 elif (
240 len(number_of_MAULs.coord("realization").points)
241 != 1
242 and len(number_of_MAULs.coord("time").points) == 1
243 ):
244 maul_depth.data[
245 mem_number, lat_point, lon_point
246 ] = np.nan
247 maul_base.data[mem_number, lat_point, lon_point] = (
248 np.nan
249 )
250 elif (
251 len(number_of_MAULs.coord("time").points) != 1
252 and len(number_of_MAULs.coord("realization").points)
253 == 1
254 ):
255 maul_depth.data[
256 time_point, lat_point, lon_point
257 ] = np.nan
258 maul_base.data[time_point, lat_point, lon_point] = (
259 np.nan
260 )
261 else:
262 maul_depth.data[lat_point, lon_point] = np.nan
263 maul_base.data[lat_point, lon_point] = np.nan
265 # Units and renaming for number, depth and base (the other case).
266 match output:
267 case "number":
268 number_of_MAULs.units = "1"
269 number_of_MAULs.rename("Number_of_MAULs")
270 num_MAULs.append(number_of_MAULs)
271 case "depth":
272 maul_depth.units = "m"
273 maul_depth.rename("MAUL_depth")
274 maul_d.append(maul_depth)
275 case _:
276 maul_base.units = "m"
277 maul_base.rename("MAUL_base_height")
278 maul_b.append(maul_base)
280 # Output data.
281 match output:
282 case "number" if len(num_MAULs) == 1:
283 return num_MAULs[0]
284 case "number":
285 return num_MAULs
286 case "depth" if len(maul_d) == 1:
287 return maul_d[0]
288 case "depth":
289 return maul_d
290 case "base" if len(maul_b) == 1:
291 return maul_b[0]
292 case _:
293 return maul_b