Coverage for test_ecape.py: 100%

79 statements  

« prev     ^ index     » next       coverage.py v7.2.7, created at 2023-06-04 23:10 -0600

1from pathlib import Path 

2 

3import numpy as np 

4from metpy.calc import dewpoint_from_specific_humidity, most_unstable_parcel 

5from metpy.units import units 

6from pytest import approx 

7 

8from src.ecape.calc import ( 

9 calc_ecape, 

10 calc_ecape_a, 

11 calc_el_height, 

12 calc_integral_arg, 

13 calc_lfc_height, 

14 calc_mse, 

15 calc_ncape, 

16 calc_psi, 

17 calc_sr_wind, 

18) 

19 

20""" 

21Note: 

22There is a ~10% difference in ECAPE between calc_ecape and Peters' published matlab scripts. 

23This is primarily due to a difference in calculated CAPE. The tests below describe other sources of error. 

24 

25Since: 

26 - the methods here are within ~1% of Peters' calculations when CAPE is equivalent in the sample data 

27 - Peters et. al. specifically mention MetPy for determining CAPE 

28 - MetPy is a reliable, open-source, and frequently used meteorological calculation package 

29  

30MetPy's CAPE calculations were chosen for ease of readability and implementation. 

31 """ 

32 

33def sample_data(): 

34 """ 

35 values via author's published matlab scripts (sigma = 1.6) 

36 https://figshare.com/articles/software/ECAPE_scripts/21859818 

37 """ 

38 

39 sounding_loc = Path("./sounding.txt") 

40 data = np.genfromtxt(sounding_loc, delimiter=",") 

41 

42 height = data[:, 0] * units("m") # height, m 

43 pressure = data[:, 1] * units("Pa") # pressure, Pa 

44 temperature = data[:, 2] * units("K") # temp, K 

45 specific_humidity = data[:, 3] * units("kg/kg") # water vapor mass fraction (specific humidity), kg / kg 

46 u_wind = data[:, 4] * units("m/s") # u wind, m / s 

47 v_wind = data[:, 5] * units("m/s") # v wind, m / s 

48 

49 dew_point_temperature = dewpoint_from_specific_humidity(pressure, temperature, specific_humidity) 

50 

51 return height, pressure, temperature, specific_humidity, u_wind, v_wind, dew_point_temperature 

52 

53 

54def test_end_to_end_ecape(): 

55 """ 

56 values via author's published matlab scripts (sigma = 1.6) 

57 https://figshare.com/articles/software/ECAPE_scripts/21859818 

58 """ 

59 

60 height, pressure, temperature, specific_humidity, u_wind, v_wind, dew_point_temperature = sample_data() 

61 

62 for cape_type in ['most_unstable', 'surface_based', 'mixed_layer']: 

63 ecape = calc_ecape(height, pressure, temperature, specific_humidity, u_wind, v_wind, cape_type) 

64 print(f"{cape_type}: {ecape}") 

65 

66 assert ecape 

67 

68 

69def test_calc_psi(): 

70 """ 

71 values via author's published matlab scripts (sigma = 1.6) 

72 https://figshare.com/articles/software/ECAPE_scripts/21859818 

73 """ 

74 

75 el_z = 11750.0 * units("m") 

76 

77 psi = calc_psi(el_z) 

78 assert psi.magnitude == approx(0.0034, rel=0.001) 

79 

80 

81def test_calc_ecape_a(): 

82 """ 

83 values via author's published matlab scripts (sigma = 1.6) 

84 https://figshare.com/articles/software/ECAPE_scripts/21859818 

85 """ 

86 

87 sr_wind = 16.662798431352986 * units('m/s') 

88 psi = 0.003401863644631 * units('dimensionless') 

89 ncape = 7.604878130037112e02 * units('m**2/s**2') 

90 cape = 3.530029673046427e03 * units('J/kg') 

91 

92 ecape_a = calc_ecape_a(sr_wind, psi, ncape, cape) 

93 assert ecape_a.magnitude == approx(3.343908138651551e03 , rel=0.0001) 

94 

95 

96def test_calc_integral_arg(): 

97 """ 

98 values via author's published matlab scripts (sigma = 1.6) 

99 https://figshare.com/articles/software/ECAPE_scripts/21859818 

100 """ 

101 

102 intarg_loc = Path("./intarg.txt") 

103 data = np.genfromtxt(intarg_loc, delimiter=",") 

104 mseo_bar = data[:, 0] * units("J/kg") 

105 mseo_star = data[:, 1] * units("J/kg") 

106 t0 = data[:, 2] * units('K') 

107 int_arg = data[:, 3] 

108 

109 integral_arg = calc_integral_arg(mseo_bar, mseo_star, t0) 

110 

111 for test, verify in zip(integral_arg, int_arg): 

112 assert test.magnitude == approx(verify, rel=0.0001) 

113 

114 

115def test_calc_ncape(): 

116 """ 

117 values via author's published matlab scripts (sigma = 1.6) 

118 https://figshare.com/articles/software/ECAPE_scripts/21859818 

119 """ 

120 

121 ncape_loc = Path("./ncape.txt") 

122 data = np.genfromtxt(ncape_loc, delimiter=",") 

123 integral_arg = data[:, 0] * units('m/s**2') 

124 height = data[:, 1] * units('m') 

125 lfc_idx = 16 

126 el_idx = 117 

127 

128 ncape = calc_ncape(integral_arg, height, lfc_idx, el_idx) 

129 assert ncape.magnitude == 7.604878130037112e02 

130 

131 

132def test_calc_mse(): 

133 """ 

134 values via author's published matlab scripts (sigma = 1.6) 

135 https://figshare.com/articles/software/ECAPE_scripts/21859818 

136 """ 

137 

138 height, pressure, temperature, specific_humidity, u_wind, v_wind, dew_point_temperature = sample_data() 

139 

140 mseo_loc = Path("./mseo.txt") 

141 data = np.genfromtxt(mseo_loc, delimiter=",") 

142 mse0_bar = data[:, 1] * units('J/kg') 

143 mse0_star = data[:, 2] * units('J/kg') 

144 

145 mseo_bar_sounding, mseo_star_sounding = calc_mse(pressure, height, temperature, specific_humidity) 

146 

147 for test, verify in zip(mseo_bar_sounding, mse0_bar): 

148 assert test.magnitude == approx(verify, rel=.005) 

149 for test, verify in zip(mseo_star_sounding, mse0_star): 

150 assert test.magnitude == approx(verify, rel=.005) 

151 

152 

153def test_calc_sr_wind(): 

154 """ 

155 values via author's published matlab scripts (sigma = 1.6) 

156 https://figshare.com/articles/software/ECAPE_scripts/21859818 

157 

158 Note: there is a difference in Bunkers right mover compenents 

159 between the author's compute_VSR.m and MetPy's metpy.calc.bunkers_storm_motion 

160 

161 For the example sounding: 

162 compute_VSR.m: C_x: 15.634223483535706, C_y: 4.74162399790177, V_sr: 16.662798431352986 

163 metpy.calc.bunkers_storm_motion: C_x: 14.635206387655197, C_y: 4.668955528029753, V_sr: 15.929554223814558 

164 ~4.6 % lower V_sr -> ~1.1 % lower ECAPE result in this case 

165 Given modest effect, no further interrogation conducted, although I expect use of WCA vs. mean to be the cause 

166 

167 Given equivalent bunkers-right components, the remainder of the function produces equivalent results 

168 MetPy's bunkers_storm_motion was used for code readability and consistency. 

169 """ 

170 

171 height, pressure, temperature, specific_humidity, u_wind, v_wind, dew_point_temperature = sample_data() 

172 

173 sr_wind = calc_sr_wind(pressure, u_wind, v_wind, height) 

174 

175 assert sr_wind.magnitude == 15.929554223814558 

176 

177 

178def test_calc_el_height(): 

179 """ 

180 Note: there is a difference in el & lfc idx/heights between the authors published scripts and 

181 the following calculations. This may be due to the selection of indexes or the el/lfc calculations. 

182 

183 COMPUTE_ECAPE.m: MU_EL_idx: 117, MU_EL: 11750 -> ECAPE 3087 

184 calc_el_height.py: el_idx: 115, el_z: 11500 -> ECAPE: 3089 

185 

186 Given a modest (~.1 %) impact on the resultant ECAPE (via the ncape & psi calculations) 

187 a process leveraging metpy.calc.el was used for code readability and consistency. 

188 """ 

189 

190 height, pressure, temperature, specific_humidity, u_wind, v_wind, dew_point_temperature = sample_data() 

191 

192 el_idx, el_z = calc_el_height(pressure, height, temperature, dew_point_temperature, most_unstable_parcel) 

193 

194 assert el_idx == 115 

195 assert el_z.magnitude == 11500 

196 

197 

198def test_calc_lfc_height(): 

199 """ 

200 Note: there is a difference in el & lfc idx/heights between the authors published scripts and 

201 the following calculations. This may be due to the selection of indexes or the el/lfc calculations. 

202 

203 COMPUTE_ECAPE.m: MU_LFC_idx: 16, MU_EL: 1650 -> ECAPE 3086 

204 calc_el_height.py: lfc_idx: 18, lfc_z: 1800 -> ECAPE: 3089 

205 

206 Given a modest (~.1 %) impact on the resultant ECAPE (via the ncape calculation) 

207 a process leveraging metpy.calc.lfc was used for code readability and consistency. 

208 """ 

209 

210 height, pressure, temperature, specific_humidity, u_wind, v_wind, dew_point_temperature = sample_data() 

211 

212 lfc_idx, lfc_z = calc_lfc_height(pressure, height, temperature, dew_point_temperature, most_unstable_parcel) 

213 

214 assert lfc_idx == 18 

215 assert lfc_z.magnitude == 1800