Coverage for src/CSET/operators/convection.py: 88%
38 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-07 15:12 +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
28logger = logging.getLogger(__name__)
31def cape_ratio(SBCAPE, MUCAPE, MUCIN, MUCIN_thresh=-75.0):
32 r"""Ratio of two fields, one filtered to allow physical values to be output.
34 Parameters
35 ----------
36 SBCAPE: Cube
37 Surface-based convective available potential energy as calculated by the
38 model. If using the UM please use STASH ``m01s20i114``
39 MUCAPE: Cube
40 Most-unstable convective available potential energy as calculated by the
41 model. If using the UM please use STASH ``m01s20i112``
42 MUCIN: Cube
43 Most-unstable convective inhibition associated with the most-unstable
44 ascent as calculated by the model. If using the UM please use STASH
45 ``m01s20i113``
46 MUCIN_thresh: float, optional, default is -75. J/kg.
47 Threshold to filter the MUCAPE by values are realistically realisable.
49 Returns
50 -------
51 Cube
53 Notes
54 -----
55 This diagnostic is based on Clark et al. (2012) [Clarketal2012]_. It is
56 based around the idea that for elevated convection the convective
57 instability is not based at the surface. This utilises two flavours of CAPE:
58 the surface-based CAPE (SBCAPE) and the most-unstable CAPE (MUCAPE). The
59 MUCAPE is filtered by the MUCIN associated with that parcel's ascent to
60 ensure that any CAPE can at least theoretically be released. The default
61 value is set at -75 J/kg but it can be changes depending on location and
62 users requirements.
64 .. math:: 1 - (\frac{SBCAPE}{MUCAPE})
66 The ratio is defined in this way such that if SBCAPE=MUCAPE the ratio will
67 equal 1. If the ratio was reversed when MUCAPE exists and SBCAPE is zero the
68 ratio would be undefined.
70 The diagnostic varies smoothly between zero and unity. A value of 0 implies
71 an environment is suitable for surface-based convection. A value of 1
72 implies an environment is suitable for elevated convection. Values between
73 imply transition convection with values closer to one imply elevated
74 convection is more likely and values closer to zero implying that
75 surface-based convection is more likely.
77 Further details about this diagnostic for elevated convection identification
78 can be found in Flack et al. (2023) [FlackCAPE2023]_.
80 Expected applicability ranges: Convective-scale models will be noisier than
81 parametrized models as they are more responsive to the convection, and thus
82 it may be more sensible to view as a larger spatial average rather than on
83 the native resolution.
85 Interpretation notes: UM stash for CAPE and CIN are calculated at the end of
86 the timestep. Therefore this diagnostic is applicable after precipitation
87 has occurred, not before as is the usual interpretation of CAPE related
88 diagnostics.
90 Examples
91 --------
92 >>> CAPE_ratios=convection.cape_ratio(
93 SBCAPE,MUCAPE,MUCIN)
94 >>> iplt.pcolormesh(CAPE_ratios[0,:,:],cmap=mpl.cm.RdBu)
95 >>> plt.gca().coastlines('10m')
96 >>> plt.colorbar()
97 >>> plt.clim(0,1)
98 >>> plt.show()
100 >>> CAPE_ratios=convection.cape_ratio(
101 SBCAPE,MUCAPE,MUCIN,MUCIN_thresh=-1.5)
102 >>> iplt.pcolormesh(CAPE_ratios[0,:,:],cmap=mpl.cm.RdBu)
103 >>> plt.gca().coastlines('10m')
104 >>> plt.clim(0,1)
105 >>> plt.colorbar()
106 >>> plt.show()
107 """
108 # Load in the data into the new arrays.
109 MUCAPE_data = copy.deepcopy(MUCAPE.data)
110 if isinstance(MUCAPE_data, np.ma.MaskedArray):
111 MUCAPE_data = MUCAPE_data.filled(np.nan)
112 # Remove all MUCAPE below MUCIN threshold.
113 MUCAPE_data[MUCIN.data <= MUCIN_thresh] = np.nan
114 with warnings.catch_warnings():
115 # Ignore possible divide by zero warnings, as they are replaced by NaNs.
116 warnings.filterwarnings("ignore", category=RuntimeWarning)
117 # Now calculate the main diagnostic.
118 EC_Flagb = 1 - (SBCAPE.data / MUCAPE_data)
119 if isinstance(EC_Flagb, np.ma.MaskedArray):
120 EC_Flagb = EC_Flagb.filled(np.nan)
121 # Filter to reduce NaN values and -inf values for plotting ease.
122 # There are multiple types of NaN values so need to convert them all to same type.
123 EC_Flagb[np.isnan(EC_Flagb)] = np.nan
124 EC_Flagb[np.isinf(EC_Flagb)] = np.nan
125 # Take the coordinates from an existing cube and replace the data.
126 cape_ratio_cube = SBCAPE.copy()
127 cape_ratio_cube.data = EC_Flagb
128 # Rename and remove STASH code.
129 cape_ratio_cube.var_name = "cape_ratio"
130 cape_ratio_cube.attributes.pop("STASH", None)
131 return cape_ratio_cube
134def inflow_layer_properties(EIB, BLheight, Orography):
135 r"""Filter to create a binary mask identifying elevated convection.
137 Parameters
138 ----------
139 EIB: Cube
140 Effective inflow layer base (precalculated or as identified by the
141 model). If using the UM please use STASH ``m01s20i119``.
142 BLheight: Cube
143 Boundary layer height (precalculated or as identified by the model). If
144 using the UM please use STASH ``m01s00i025``.
145 Orography: Cube
146 Model or actual orography, expected to be 2 dimensional. If 3 or 4
147 dimensional cube given converts to 2 dimensions assuming static
148 orography field in ensemble realization and time. If using the UM please
149 use STASH ``m01s00i033``.
151 Returns
152 -------
153 Cube
155 Notes
156 -----
157 This diagnostic is based on the concept of an effective inflow layer. This
158 concept was first introduced by Thompson et al. (2007) [Thompsonetal2007]_.
159 The inflow layer defined the region of air that is most likely to be
160 ingested into the convective event. It is defined by thresholding the CAPE
161 and CIN values: CAPE > 100 J/kg and \|CIN\| < 250 J/kg.
163 To turn this into a diagnostic for elevated convection the inflow layer base
164 is filtered against the boundary layer height. The model orography is added
165 to the boundary layer height to ensure reference height consistency as the
166 BL height is defined above ground level and the inflow layer base is defined
167 above sea level in the model output.
169 .. math:: EIB > BLheight + Orography
171 This is a binary diagnostic. It has a value of 0 to imply the environment is
172 suitable for surface-based convection. It has a value of 1 to indicate the
173 environment is suitable to produce elevated convection.
175 Further details about this diagnostic for elevated convection identification
176 can be found in Flack et al. (2023) [Flackinf2023]_.
178 Expected applicability ranges: Convective-scale models will be noisier than
179 parametrized models as they are more responsive to the convection, and thus
180 it may be more sensible to view as a larger spatial average rather than at
181 native resolution.
183 Interpretation notes: The effective inflow layer base diagnostic from UM
184 STASH is dependent upon the UM CAPE and CIN diagnostics. These diagnostics
185 are calculated at the end of the timestep. Therefore this diagnostic is
186 applicable after precipitation has occurred, not before as is the usual
187 interpretation of CAPE related diagnostics.
189 You might encounter warnings with the following text ``Orography assumed not
190 to vary with ensemble member.`` or ``Orography assumed not to vary with time
191 and ensemble member.`` these warnings are expected when the orography files
192 are not 2-dimensional, and do not cause any problems unless ordering is not
193 as expected.
195 Examples
196 --------
197 >>> Inflow_properties=convection.inflow_layer_properties(EIB,BLheight,Orography)
198 >>> iplt.pcolormesh(Inflow_properties[0,:,:],cmap=mpl.cm.Purples)
199 >>> plt.gca().coastlines('10m')
200 >>> plt.colorbar()
201 >>> plt.clim(0,1)
202 >>> plt.show()
204 """
205 # Setup new array for output of the diagnostic.
206 EC_Flagd = np.zeros(EIB.shape)
207 # Check dimensions for Orography cube and replace with 2D array if not 2D.
208 if Orography.ndim == 3: 208 ↛ 209line 208 didn't jump to line 209 because the condition on line 208 was never true
209 Orography = Orography.slices_over("realization").next()
210 logger.warning("Orography assumed not to vary with ensemble member")
211 elif Orography.ndim == 4: 211 ↛ 212line 211 didn't jump to line 212 because the condition on line 211 was never true
212 Orography = Orography.slices_over(("time", "realization")).next()
213 logger.warning("Orography assumed not to vary with time or ensemble member. ")
214 # Masked arrays are not respected, so convert masked values into NaNs.
215 if isinstance(EIB.data, np.ma.MaskedArray):
216 EIB.data = EIB.data.filled(np.nan)
217 # Change points where Effective inflow layer base is larger than boundary
218 # layer height to 1 implying elevated convection.
219 EC_Flagd[EIB.data > (BLheight.data + Orography.data)] = 1.0
220 # Take the coordinates from an existing cube and replace the data.
221 inflow_properties_cube = EIB.copy()
222 inflow_properties_cube.data = EC_Flagd
223 # Rename and remove STASH code.
224 inflow_properties_cube.var_name = "inflow_layer_properties"
225 inflow_properties_cube.attributes.pop("STASH", None)
226 return inflow_properties_cube