Sounding Calculation Examples
Use functions from metpy.calc
to perform a number of calculations using sounding data.
The code below uses example data to perform many sounding calculations for a severe weather event on May 22, 2011 from the Norman, OK sounding.
Imports
import numpy as np
import pandas as pd
import metpy.calc as mpcalc
from metpy.cbook import get_test_data
from metpy.units import units
Effective Shear Algorithm for use in Supercell Composite Calculation
def effective_layer(pressure, temperature, dewpoint, height, height_layer=False):
"""A function that determines the effective inflow layer for a convective sounding.
Uses the default values of Thompason et al. (2004) for CAPE (100 J/kg) and CIN (-250 J/kg).
Input:
- pressure: sounding pressure with units
- temperature: sounding temperature with units
- dewpoint: sounding dewpoint temperature with units
- height: sounding heights with units
Returns:
- pbot/hbot, ptop/htop: pressure/height of the bottom level,
pressure/height of the top level
"""
from metpy.calc import cape_cin, parcel_profile
from metpy.units import units
pbot = None
for i in range(pressure.shape[0]):
prof = parcel_profile(pressure[i:], temperature[i], dewpoint[i])
sbcape, sbcin = cape_cin(pressure[i:], temperature[i:], dewpoint[i:], prof)
if sbcape >= 100 * units('J/kg') and sbcin > -250 * units('J/kg'):
pbot = pressure[i]
hbot = height[i]
bot_idx = i
break
if not pbot:
return None, None
for i in range(bot_idx + 1, pressure.shape[0]):
prof = parcel_profile(pressure[i:], temperature[i], dewpoint[i])
sbcape, sbcin = cape_cin(pressure[i:], temperature[i:], dewpoint[i:], prof)
if sbcape < 100 * units('J/kg') or sbcin < -250 * units('J/kg'):
ptop = pressure[i]
htop = height[i]
break
if height_layer:
return hbot, htop
else:
return pbot, ptop
Obtain Data and Format
Upper air data can be obtained using the siphon package, but for this example we will use some of MetPy’s sample data.
as_file_obj=False), skiprows=5, usecols=[0, 1, 2, 3, 6, 7], names=col_names)
is necessary due to the formatting of the MetPy sample data. This formatting is not needed when using upper air data obtained via Siphon. Obtaining data with Siphon will be covered in a later notebook.
col_names = ['pressure', 'height', 'temperature', 'dewpoint', 'direction', 'speed']
sounding_data = pd.read_fwf(get_test_data('20110522_OUN_12Z.txt', as_file_obj=False),
skiprows=7, usecols=[0, 1, 2, 3, 6, 7], names=col_names)
# Drop any rows with all not a number (NaN) values for temperature, dewpoint, and winds
sounding_data = sounding_data.dropna(subset=('temperature', 'dewpoint', 'direction', 'speed'),
how='all').reset_index(drop=True)
Downloading file '20110522_OUN_12Z.txt' from 'https://github.com/Unidata/MetPy/raw/v1.6.2/staticdata/20110522_OUN_12Z.txt' to '/home/runner/.cache/metpy/v1.6.2'.
Assign Units
We will pull the data out of the example dataset into individual variables and assign units. This is explained in further detain in the Simple Sounding notebook and in the Metpy documentation.
pres = sounding_data['pressure'].values * units.hPa
temp = sounding_data['temperature'].values * units.degC
dewpoint = sounding_data['dewpoint'].values * units.degC
wind_speed = sounding_data['speed'].values * units.knots
wind_dir = sounding_data['direction'].values * units.degrees
u, v = mpcalc.wind_components(wind_speed, wind_dir)
height = sounding_data['height'].values * units.meter
Compute the wind components
u, v = mpcalc.wind_components(wind_speed, wind_dir)
Compute common sounding index parameters
ctotals = mpcalc.cross_totals(pres, temp, dewpoint)
kindex = mpcalc.k_index(pres, temp, dewpoint)
showalter = mpcalc.showalter_index(pres, temp, dewpoint)
total_totals = mpcalc.total_totals_index(pres, temp, dewpoint)
vert_totals = mpcalc.vertical_totals(pres, temp)
Compture the parcel profile for a surface-based parcel
prof = mpcalc.parcel_profile(pres, temp[0], dewpoint[0])
Compute the corresponding Lifted Index (LI), Convective Available Potential Energy (CAPE), Convective Inhibition (CIN) values for a surface parcel
lift_index = mpcalc.lifted_index(pres, temp, prof)
cape, cin = mpcalc.cape_cin(pres, temp, dewpoint, prof)
Determine the Lifted Condensation Level (LCL), Level of Free Convection (LFC), and Equilibrium Level (EL) for our surface parcel
lclp, lclt = mpcalc.lcl(pres[0], temp[0], dewpoint[0])
lfcp, _ = mpcalc.lfc(pres, temp, dewpoint)
el_pressure, _ = mpcalc.el(pres, temp, dewpoint, prof)
Compute the characteristics of a mean layer parcel (50-hPa depth)
ml_t, ml_td = mpcalc.mixed_layer(pres, temp, dewpoint, depth=50 * units.hPa)
ml_p, _, _ = mpcalc.mixed_parcel(pres, temp, dewpoint, depth=50 * units.hPa)
mlcape, mlcin = mpcalc.mixed_layer_cape_cin(pres, temp, prof, depth=50 * units.hPa)
Compute the characteristics of the most unstable parcel (50-hPa depth)
mu_p, mu_t, mu_td, _ = mpcalc.most_unstable_parcel(pres, temp, dewpoint, depth=50 * units.hPa)
mucape, mucin = mpcalc.most_unstable_cape_cin(pres, temp, dewpoint, depth=50 * units.hPa)
Compute the Bunkers Storm Motion vector and use to calculate the critical angle
(u_storm, v_storm), *_ = mpcalc.bunkers_storm_motion(pres, u, v, height)
critical_angle = mpcalc.critical_angle(pres, u, v, height, u_storm, v_storm)
Work on the calculations needed to compute the significant tornado parameter
# Estimate height of LCL in meters from hydrostatic thickness
new_p = np.append(pres[pres > lclp], lclp)
new_t = np.append(temp[pres > lclp], lclt)
lcl_height = mpcalc.thickness_hydrostatic(new_p, new_t)
# Compute Surface-based CAPE
sbcape, _ = mpcalc.surface_based_cape_cin(pres, temp, dewpoint)
# Compute SRH, given a motion vector toward the NE at 9.9 m/s
*_, total_helicity = mpcalc.storm_relative_helicity(height, u, v, depth=1 * units.km,
storm_u=u_storm, storm_v=v_storm)
# Copmute Bulk Shear components and then magnitude
ubshr, vbshr = mpcalc.bulk_shear(pres, u, v, height=height, depth=6 * units.km)
bshear = mpcalc.wind_speed(ubshr, vbshr)
# Use all computed pieces to calculate the Significant Tornado parameter
sig_tor = mpcalc.significant_tornado(sbcape, lcl_height,
total_helicity, bshear).to_base_units()
Compute the supercell composite parameter, if possible
# Determine the top and bottom of the effective layer using our own function
hbot, htop = effective_layer(pres, temp, dewpoint, height, height_layer=True)
# Perform the calculation of supercell composite if an effective layer exists
if hbot:
esrh = mpcalc.storm_relative_helicity(height, u, v, depth=htop - hbot, bottom=hbot)
eubshr, evbshr = mpcalc.bulk_shear(pres, u, v, height=height, depth=htop - hbot, bottom=hbot)
ebshear = mpcalc.wind_speed(eubshr, evbshr)
super_comp = mpcalc.supercell_composite(mucape, esrh[0], ebshear)
else:
super_comp = np.nan
Print Important Sounding Parameters
print('Important Sounding Parameters for KOUN on 22 Mary 2011 12 UTC')
print()
print(f' CAPE: {cape:.2f}')
print(f' CIN: {cin:.2f}')
print(f'LCL Pressure: {lclp:.2f}')
print(f'LFC Pressure: {lfcp:.2f}')
print(f' EL Pressure: {el_pressure:.2f}')
print()
print(f' Lifted Index: {lift_index:.2f}')
print(f' K-Index: {kindex:.2f}')
print(f'Showalter Index: {showalter:.2f}')
print(f' Cross Totals: {ctotals:.2f}')
print(f' Total Totals: {total_totals:.2f}')
print(f'Vertical Totals: {vert_totals:.2f}')
print()
print('Mixed Layer - Lowest 50-hPa')
print(f' ML Temp: {ml_t:.2f}')
print(f' ML Dewp: {ml_td:.2f}')
print(f' ML CAPE: {mlcape:.2f}')
print(f' ML CIN: {mlcin:.2f}')
print()
print('Most Unstable - Lowest 50-hPa')
print(f' MU Temp: {mu_t:.2f}')
print(f' MU Dewp: {mu_td:.2f}')
print(f' MU Pressure: {mu_p:.2f}')
print(f' MU CAPE: {mucape:.2f}')
print(f' MU CIN: {mucin:.2f}')
print()
print('Bunkers Storm Motion Vector')
print(f' u_storm: {u_storm:.2f}')
print(f' v_storm: {v_storm:.2f}')
print(f'Critical Angle: {critical_angle:.2f}')
print()
print(f'Storm Relative Helicity: {total_helicity:.2f}')
print(f'Significant Tornado Parameter: {sig_tor:.2f}')
print(f'Supercell Composite Parameter: {super_comp:.2f}')
Important Sounding Parameters for KOUN on 22 Mary 2011 12 UTC
CAPE: 3223.89 joule / kilogram
CIN: -96.26 joule / kilogram
LCL Pressure: 949.09 hectopascal
LFC Pressure: 735.99 hectopascal
EL Pressure: 194.72 hectopascal
Lifted Index: [-6.96] delta_degree_Celsius
K-Index: 22.10 degree_Celsius
Showalter Index: [-0.08] delta_degree_Celsius
Cross Totals: 17.10 delta_degree_Celsius
Total Totals: 50.20 delta_degree_Celsius
Vertical Totals: 33.10 delta_degree_Celsius
Mixed Layer - Lowest 50-hPa
ML Temp: 20.99 degree_Celsius
ML Dewp: 20.55 degree_Celsius
ML CAPE: 3254.17 joule / kilogram
ML CIN: -138.20 joule / kilogram
Most Unstable - Lowest 50-hPa
MU Temp: 20.40 degree_Celsius
MU Dewp: 20.40 degree_Celsius
MU Pressure: 925.00 hectopascal
MU CAPE: 3693.64 joule / kilogram
MU CIN: -60.50 joule / kilogram
Bunkers Storm Motion Vector
u_storm: 21.85 knot
v_storm: 4.55 knot
Critical Angle: 67.33 degree
Storm Relative Helicity: 279.49 meter ** 2 / second ** 2
Significant Tornado Parameter: [4.60]
Supercell Composite Parameter: [9.07]