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

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 """

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

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

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

76 if output not in ("number", "base", "depth"):

77 raise ValueError(

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

79 )

81 for cube in iter_maybe(cubes):

82 # Check for binary fields.

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

84 raise ValueError(

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

86 )

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

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

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

90 number_of_MAULs.data[:] = 0.0

91 maul_depth = number_of_MAULs.copy()

92 maul_base = number_of_MAULs.copy()

93 # Loop over realization.

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

95 # Loop over time.

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

97 # Loop over latitude.

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

99 # Loop over longitude.

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

101 # Label each object in the vertical.

102 labels = label(lon.core_data())

103 # Finds the number of MAULs present based upon the

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

105 # the value is set to zero.

106 # The code checks for whether there are multiple

107 # realization and/or time points for correct

108 # indexing of the output data and applies accordingly.

109 if (

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

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

112 ):

113 number_of_MAULs.data[

114 mem_number, time_point, lat_point, lon_point

115 ] = np.max(labels)

116 elif (

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

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

119 ):

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

121 np.max(labels)

122 )

123 elif (

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

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

126 ):

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

128 np.max(labels)

129 )

130 else:

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

132 if output != "number":

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

134 # using cube metadata.

135 maul_start = []

136 maul_end = []

137 maul_dep = []

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

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

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

141 # Find all vertical indices belonging to a MAUL.

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

143 # Find the height at the base of the MAUL

144 # (lowest level).

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

146 maul_range[0][0]

147 ]

148 # Find the height at the top of the MAUL

149 # (highest level).

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

151 maul_range[0][-1]

152 ]

153 # Calculate the MAUL depth, and store

154 # base and top heights.

155 maul_dep.append(maul_end_point - maul_start_point)

156 maul_start.append(maul_start_point)

157 maul_end.append(maul_end_point)

158 try:

159 # Idendtify where the deepest MAUL is.

160 index = int(

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

162 )

163 # As with number the code checks for whether

164 # there are multiple realization and/or time

165 # points for correct indexing of the output data

166 # and applies accordingly.

167 if (

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

169 != 1

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

171 ):

172 # Store the deepest MAUL.

173 maul_depth.data[

174 mem_number, time_point, lat_point, lon_point

175 ] = np.max(maul_dep)

176 # Store the base height of the deepest MAUL.

177 maul_base.data[

178 mem_number, time_point, lat_point, lon_point

179 ] = maul_start[index]

180 elif (

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

182 != 1

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

184 ):

185 maul_depth.data[

186 mem_number, lat_point, lon_point

187 ] = np.max(maul_dep)

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

189 maul_start[index]

190 )

191 elif (

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

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

194 == 1

195 ):

196 maul_depth.data[

197 time_point, lat_point, lon_point

198 ] = np.max(maul_dep)

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

200 maul_start[index]

201 )

202 else:

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

204 maul_dep

205 )

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

207 index

208 ]

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

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

211 # The dimensionality logic for output data is identical

212 # to that used previously.

213 except ValueError:

214 if (

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, time_point, lat_point, lon_point

221 ] = np.nan

222 maul_base.data[

223 mem_number, time_point, lat_point, lon_point

224 ] = np.nan

225 elif (

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

227 != 1

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

229 ):

230 maul_depth.data[

231 mem_number, lat_point, lon_point

232 ] = np.nan

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

234 np.nan

235 )

236 elif (

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

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

239 == 1

240 ):

241 maul_depth.data[

242 time_point, lat_point, lon_point

243 ] = np.nan

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

245 np.nan

246 )

247 else:

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

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

250

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

252 match output:

253 case "number":

254 number_of_MAULs.units = "1"

255 number_of_MAULs.rename("Number_of_MAULs")

256 num_MAULs.append(number_of_MAULs)

257 case "depth":

258 maul_depth.units = "m"

259 maul_depth.rename("MAUL_depth")

260 maul_d.append(maul_depth)

261 case _:

262 maul_base.units = "m"

263 maul_base.rename("MAUL_base_height")

264 maul_b.append(maul_base)

265

266 # Output data.

267 match output:

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

269 return num_MAULs[0]

270 case "number":

271 return num_MAULs

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

273 return maul_d[0]

274 case "depth":

275 return maul_d

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

277 return maul_b[0]

278 case _:

279 return maul_b

Coverage for src / CSET / operators / precipitation.py: 100%

93 statements