Coverage for test_ecape.py: 100%
79 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-06-04 23:10 -0600
« prev ^ index » next coverage.py v7.2.7, created at 2023-06-04 23:10 -0600
1from pathlib import Path
3import numpy as np
4from metpy.calc import dewpoint_from_specific_humidity, most_unstable_parcel
5from metpy.units import units
6from pytest import approx
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)
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.
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
30MetPy's CAPE calculations were chosen for ease of readability and implementation.
31 """
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 """
39 sounding_loc = Path("./sounding.txt")
40 data = np.genfromtxt(sounding_loc, delimiter=",")
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
49 dew_point_temperature = dewpoint_from_specific_humidity(pressure, temperature, specific_humidity)
51 return height, pressure, temperature, specific_humidity, u_wind, v_wind, dew_point_temperature
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 """
60 height, pressure, temperature, specific_humidity, u_wind, v_wind, dew_point_temperature = sample_data()
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}")
66 assert ecape
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 """
75 el_z = 11750.0 * units("m")
77 psi = calc_psi(el_z)
78 assert psi.magnitude == approx(0.0034, rel=0.001)
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 """
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')
92 ecape_a = calc_ecape_a(sr_wind, psi, ncape, cape)
93 assert ecape_a.magnitude == approx(3.343908138651551e03 , rel=0.0001)
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 """
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]
109 integral_arg = calc_integral_arg(mseo_bar, mseo_star, t0)
111 for test, verify in zip(integral_arg, int_arg):
112 assert test.magnitude == approx(verify, rel=0.0001)
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 """
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
128 ncape = calc_ncape(integral_arg, height, lfc_idx, el_idx)
129 assert ncape.magnitude == 7.604878130037112e02
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 """
138 height, pressure, temperature, specific_humidity, u_wind, v_wind, dew_point_temperature = sample_data()
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')
145 mseo_bar_sounding, mseo_star_sounding = calc_mse(pressure, height, temperature, specific_humidity)
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)
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
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
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
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 """
171 height, pressure, temperature, specific_humidity, u_wind, v_wind, dew_point_temperature = sample_data()
173 sr_wind = calc_sr_wind(pressure, u_wind, v_wind, height)
175 assert sr_wind.magnitude == 15.929554223814558
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.
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
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 """
190 height, pressure, temperature, specific_humidity, u_wind, v_wind, dew_point_temperature = sample_data()
192 el_idx, el_z = calc_el_height(pressure, height, temperature, dew_point_temperature, most_unstable_parcel)
194 assert el_idx == 115
195 assert el_z.magnitude == 11500
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.
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
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 """
210 height, pressure, temperature, specific_humidity, u_wind, v_wind, dew_point_temperature = sample_data()
212 lfc_idx, lfc_z = calc_lfc_height(pressure, height, temperature, dew_point_temperature, most_unstable_parcel)
214 assert lfc_idx == 18
215 assert lfc_z.magnitude == 1800