Coverage for src/CSET/operators/convection.py: 100%
37 statements
« prev ^ index » next coverage.py v7.10.6, created at 2025-09-05 21:08 +0000
1# © Crown copyright, Met Office (2022-2024) 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"""A module containing different diagnostics for convection.
17The diagnostics are calculated from output from the Unified Model, although
18precalculated values in the required input form may also be used.
20"""
22import copy
23import logging
24import warnings
26import numpy as np
29def cape_ratio(SBCAPE, MUCAPE, MUCIN, MUCIN_thresh=-75.0):
30 r"""Ratio of two fields, one filtered to allow physical values to be output.
32 Parameters
33 ----------
34 SBCAPE: Cube
35 Surface-based convective available potential energy as calculated by the
36 model. If using the UM please use STASH ``m01s20i114``
37 MUCAPE: Cube
38 Most-unstable convective available potential energy as calculated by the
39 model. If using the UM please use STASH ``m01s20i112``
40 MUCIN: Cube
41 Most-unstable convective inhibition associated with the most-unstable
42 ascent as calculated by the model. If using the UM please use STASH
43 ``m01s20i113``
44 MUCIN_thresh: float, optional, default is -75. J/kg.
45 Threshold to filter the MUCAPE by values are realistically realisable.
47 Returns
48 -------
49 Cube
51 Notes
52 -----
53 This diagnostic is based on Clark et al. (2012) [Clarketal2012]_. It is
54 based around the idea that for elevated convection the convective
55 instability is not based at the surface. This utilises two flavours of CAPE:
56 the surface-based CAPE (SBCAPE) and the most-unstable CAPE (MUCAPE). The
57 MUCAPE is filtered by the MUCIN associated with that parcel's ascent to
58 ensure that any CAPE can at least theoretically be released. The default
59 value is set at -75 J/kg but it can be changes depending on location and
60 users requirements.
62 .. math:: 1 - (\frac{SBCAPE}{MUCAPE})
64 The ratio is defined in this way such that if SBCAPE=MUCAPE the ratio will
65 equal 1. If the ratio was reversed when MUCAPE exists and SBCAPE is zero the
66 ratio would be undefined.
68 The diagnostic varies smoothly between zero and unity. A value of 0 implies
69 an environment is suitable for surface-based convection. A value of 1
70 implies an environment is suitable for elevated convection. Values between
71 imply transition convection with values closer to one imply elevated
72 convection is more likely and values closer to zero implying that
73 surface-based convection is more likely.
75 Further details about this diagnostic for elevated convection identification
76 can be found in Flack et al. (2023) [FlackCAPE2023]_.
78 Expected applicability ranges: Convective-scale models will be noisier than
79 parametrized models as they are more responsive to the convection, and thus
80 it may be more sensible to view as a larger spatial average rather than on
81 the native resolution.
83 Interpretation notes: UM stash for CAPE and CIN are calculated at the end of
84 the timestep. Therefore this diagnostic is applicable after precipitation
85 has occurred, not before as is the usual interpretation of CAPE related
86 diagnostics.
88 References
89 ----------
90 .. [Clarketal2012] Clark, A. J., Kain J. S., Marsh P. T., Correia J., Xue
91 M., and Kong F., (2012) "Forecasting tornado pathlengths using a
92 three-dimensional object identification algorithm applied to
93 convection-allowing forecasts." Weather and Forecasting, vol. 27,
94 1090–1113, doi: 10.1175/WAF-D-11-00147.1
95 .. [FlackCAPE2023] Flack, D.L.A., Lehnert, M., Lean, H.W., and Willington,
96 S. (2023) "Characteristics of Diagnostics for Identifying Elevated
97 Convection over the British Isles in a Convection-Allowing Model."
98 Weather and Forecasting, vol. 30, 1079-1094, doi: 10.1175/WAF-D-22-0219.1
100 Examples
101 --------
102 >>> CAPE_ratios=convection.cape_ratio(
103 SBCAPE,MUCAPE,MUCIN)
104 >>> iplt.pcolormesh(CAPE_ratios[0,:,:],cmap=mpl.cm.RdBu)
105 >>> plt.gca().coastlines('10m')
106 >>> plt.colorbar()
107 >>> plt.clim(0,1)
108 >>> plt.show()
110 >>> CAPE_ratios=convection.cape_ratio(
111 SBCAPE,MUCAPE,MUCIN,MUCIN_thresh=-1.5)
112 >>> iplt.pcolormesh(CAPE_ratios[0,:,:],cmap=mpl.cm.RdBu)
113 >>> plt.gca().coastlines('10m')
114 >>> plt.clim(0,1)
115 >>> plt.colorbar()
116 >>> plt.show()
117 """
118 # Load in the data into the new arrays.
119 MUCAPE_data = copy.deepcopy(MUCAPE.data)
120 if isinstance(MUCAPE_data, np.ma.MaskedArray):
121 MUCAPE_data = MUCAPE_data.filled(np.nan)
122 # Remove all MUCAPE below MUCIN threshold.
123 MUCAPE_data[MUCIN.data <= MUCIN_thresh] = np.nan
124 with warnings.catch_warnings():
125 # Ignore possible divide by zero warnings, as they are replaced by NaNs.
126 warnings.filterwarnings("ignore", category=RuntimeWarning)
127 # Now calculate the main diagnostic.
128 EC_Flagb = 1 - (SBCAPE.data / MUCAPE_data)
129 if isinstance(EC_Flagb, np.ma.MaskedArray):
130 EC_Flagb = EC_Flagb.filled(np.nan)
131 # Filter to reduce NaN values and -inf values for plotting ease.
132 # There are multiple types of NaN values so need to convert them all to same type.
133 EC_Flagb[np.isnan(EC_Flagb)] = np.nan
134 EC_Flagb[np.isinf(EC_Flagb)] = np.nan
135 # Take the coordinates from an existing cube and replace the data.
136 cape_ratio_cube = SBCAPE.copy()
137 cape_ratio_cube.data = EC_Flagb
138 # Rename and remove STASH code.
139 cape_ratio_cube.var_name = "cape_ratio"
140 cape_ratio_cube.attributes.pop("STASH", None)
141 return cape_ratio_cube
144def inflow_layer_properties(EIB, BLheight, Orography):
145 r"""Filter to create a binary mask identifying elevated convection.
147 Parameters
148 ----------
149 EIB: Cube
150 Effective inflow layer base (precalculated or as identified by the
151 model). If using the UM please use STASH ``m01s20i119``.
152 BLheight: Cube
153 Boundary layer height (precalculated or as identified by the model). If
154 using the UM please use STASH ``m01s00i025``.
155 Orography: Cube
156 Model or actual orography, expected to be 2 dimensional. If 3 or 4
157 dimensional cube given converts to 2 dimensions assuming static
158 orography field in ensemble realization and time. If using the UM please
159 use STASH ``m01s00i033``.
161 Returns
162 -------
163 Cube
165 Notes
166 -----
167 This diagnostic is based on the concept of an effective inflow layer. This
168 concept was first introduced by Thompson et al. (2007) [Thompsonetal2007]_.
169 The inflow layer defined the region of air that is most likely to be
170 ingested into the convective event. It is defined by thresholding the CAPE
171 and CIN values: CAPE > 100 J/kg and \|CIN\| < 250 J/kg.
173 To turn this into a diagnostic for elevated convection the inflow layer base
174 is filtered against the boundary layer height. The model orography is added
175 to the boundary layer height to ensure reference height consistency as the
176 BL height is defined above ground level and the inflow layer base is defined
177 above sea level in the model output.
179 .. math:: EIB > BLheight + Orography
181 This is a binary diagnostic. It has a value of 0 to imply the environment is
182 suitable for surface-based convection. It has a value of 1 to indicate the
183 environment is suitable to produce elevated convection.
185 Further details about this diagnostic for elevated convection identification
186 can be found in Flack et al. (2023) [Flackinf2023]_.
188 Expected applicability ranges: Convective-scale models will be noisier than
189 parametrized models as they are more responsive to the convection, and thus
190 it may be more sensible to view as a larger spatial average rather than at
191 native resolution.
193 Interpretation notes: The effective inflow layer base diagnostic from UM
194 STASH is dependent upon the UM CAPE and CIN diagnostics. These diagnostics
195 are calculated at the end of the timestep. Therefore this diagnostic is
196 applicable after precipitation has occurred, not before as is the usual
197 interpretation of CAPE related diagnostics.
199 You might encounter warnings with the following text ``Orography assumed not
200 to vary with ensemble member.`` or ``Orography assumed not to vary with time
201 and ensemble member.`` these warnings are expected when the orography files
202 are not 2-dimensional, and do not cause any problems unless ordering is not
203 as expected.
205 References
206 ----------
207 .. [Thompsonetal2007] Thompson, R. L. Mead, C. M., and Edwards, R., (2007)
208 "Effective Storm-Relative Helicity and Bulk Shear in Supercell
209 Thunderstorm Environments." Weather and Forecasting, vol. 22, 102-115,
210 doi: 10.1175/WAF969.1
211 .. [Flackinf2023] Flack, D.L.A., Lehnert, M., Lean, H.W., and Willington, S.
212 (2023) "Characteristics of Diagnostics for Identifying Elevated
213 Convection over the British Isles in a Convection-Allowing Model."
214 Weather and Forecasting, vol. 30, 1079-1094, doi: 10.1175/WAF-D-22-0219.1
216 Examples
217 --------
218 >>> Inflow_properties=convection.inflow_layer_properties(EIB,BLheight,Orography)
219 >>> iplt.pcolormesh(Inflow_properties[0,:,:],cmap=mpl.cm.Purples)
220 >>> plt.gca().coastlines('10m')
221 >>> plt.colorbar()
222 >>> plt.clim(0,1)
223 >>> plt.show()
225 """
226 # Setup new array for output of the diagnostic.
227 EC_Flagd = np.zeros(EIB.shape)
228 # Check dimensions for Orography cube and replace with 2D array if not 2D.
229 if Orography.ndim == 3:
230 Orography = Orography.slices_over("realization").next()
231 logging.warning("Orography assumed not to vary with ensemble member")
232 elif Orography.ndim == 4:
233 Orography = Orography.slices_over(("time", "realization")).next()
234 logging.warning("Orography assumed not to vary with time or ensemble member. ")
235 # Masked arrays are not respected, so convert masked values into NaNs.
236 if isinstance(EIB.data, np.ma.MaskedArray):
237 EIB.data = EIB.data.filled(np.nan)
238 # Change points where Effective inflow layer base is larger than boundary
239 # layer height to 1 implying elevated convection.
240 EC_Flagd[EIB.data > (BLheight.data + Orography.data)] = 1.0
241 # Take the coordinates from an existing cube and replace the data.
242 inflow_properties_cube = EIB.copy()
243 inflow_properties_cube.data = EC_Flagd
244 # Rename and remove STASH code.
245 inflow_properties_cube.var_name = "inflow_layer_properties"
246 inflow_properties_cube.attributes.pop("STASH", None)
247 return inflow_properties_cube