Sounding Calculation Examples

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

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

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)

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

u, v = mpcalc.wind_components(wind_speed, wind_dir)

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)

prof = mpcalc.parcel_profile(pres, temp[0], dewpoint[0])

lift_index = mpcalc.lifted_index(pres, temp, prof)
cape, cin = mpcalc.cape_cin(pres, temp, dewpoint, prof)

lclp, lclt = mpcalc.lcl(pres[0], temp[0], dewpoint[0])
lfcp, _ = mpcalc.lfc(pres, temp, dewpoint)
el_pressure, _ = mpcalc.el(pres, temp, dewpoint, prof)

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)

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)

(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)

# 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()

# 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 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}')