Calculate derived atmospheric quantities using metpy.calc — thermodynamics, soundings, kinematics, boundary layer, smoothing, and more.
Help the user calculate derived atmospheric quantities with MetPy. Reference: https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.html
conda install -c conda-forge metpy
MetPy requires quantities to carry units via pint. Always attach units before calling metpy.calc functions.
import metpy.calc as mpcalc
from metpy.units import units
import numpy as np
T = np.array([20., 15., 10.]) * units.degC
p = np.array([1000., 850., 700.]) * units.hPa
rh = np.array([80., 70., 60.]) * units.percent
With xarray, use to attach units from CF attributes:
.metpy.quantify()import xarray as xr
ds = xr.open_dataset("file.nc").metpy.quantify()
T = ds["temperature"]
# Potential temperature
theta = mpcalc.potential_temperature(p, T)
# Temperature from potential temperature
T = mpcalc.temperature_from_potential_temperature(p, theta)
# Dry lapse rate (parcel)
T_parcel = mpcalc.dry_lapse(p, T[0])
# Air density
rho = mpcalc.density(p, T, mixing_ratio)
# Dry/moist static energy
dse = mpcalc.dry_static_energy(height, T)
mse = mpcalc.moist_static_energy(height, T, specific_humidity)
# Layer thickness (hypsometric)
dz = mpcalc.thickness_hydrostatic(p, T)
# Geopotential <-> height
z = mpcalc.geopotential_to_height(geopotential)
phi = mpcalc.height_to_geopotential(z)
# Pressure weighted mean
T_mean = mpcalc.mean_pressure_weighted(p, T, height=z, bottom=z[0], depth=100*units.hPa)
# Saturation vapor pressure
es = mpcalc.saturation_vapor_pressure(T)
# Dewpoint
Td = mpcalc.dewpoint(vapor_pressure)
Td = mpcalc.dewpoint_from_relative_humidity(T, rh)
Td = mpcalc.dewpoint_from_specific_humidity(specific_humidity, p)
# Relative humidity
rh = mpcalc.relative_humidity_from_dewpoint(T, Td)
rh = mpcalc.relative_humidity_from_mixing_ratio(mixing_ratio, T, p)
rh = mpcalc.relative_humidity_from_specific_humidity(specific_humidity, T, p)
# Mixing ratio
w = mpcalc.mixing_ratio(vapor_pressure, p)
w = mpcalc.mixing_ratio_from_relative_humidity(T, p, rh)
w = mpcalc.saturation_mixing_ratio(p, T)
# Specific humidity
q = mpcalc.specific_humidity_from_mixing_ratio(w)
q = mpcalc.specific_humidity_from_dewpoint(p, Td)
# Equivalent / virtual potential temperature
theta_e = mpcalc.equivalent_potential_temperature(p, T, Td)
theta_es = mpcalc.saturation_equivalent_potential_temperature(p, T, Td)
theta_v = mpcalc.virtual_potential_temperature(p, T, w)
# Virtual temperature
Tv = mpcalc.virtual_temperature(T, w)
Tv = mpcalc.virtual_temperature_from_dewpoint(p, T, Td)
# Wet bulb temperature
Tw = mpcalc.wet_bulb_temperature(p, T, Td)
# Precipitable water
pw = mpcalc.precipitable_water(p, Td)
# Vertical velocity conversion
w = mpcalc.vertical_velocity(omega, p, T) # omega -> w (m/s)
om = mpcalc.vertical_velocity_pressure(w, p, T) # w -> omega (Pa/s)
# Parcel profile
prof = mpcalc.parcel_profile(p, T[0], Td[0])
prof, p_lcl, T_lcl = mpcalc.parcel_profile_with_lcl(p, T, Td)
# Key levels
lcl_p, lcl_T = mpcalc.lcl(p[0], T[0], Td[0])
lfc_p, lfc_T = mpcalc.lfc(p, T, Td)
el_p, el_T = mpcalc.el(p, T, Td)
# CAPE / CIN
cape, cin = mpcalc.cape_cin(p, T, Td, prof)
cape, cin = mpcalc.surface_based_cape_cin(p, T, Td)
cape, cin = mpcalc.mixed_layer_cape_cin(p, T, Td)
cape, cin = mpcalc.most_unstable_cape_cin(p, T, Td)
dcape = mpcalc.downdraft_cape(p, T, Td)
# Mixed layer parcel
p_ml, T_ml, Td_ml = mpcalc.mixed_parcel(p, T, Td)
T_ml = mpcalc.mixed_layer(p, T)
# Stability indices
k = mpcalc.k_index(p, T, Td)
tt = mpcalc.total_totals_index(p, T, Td)
ct = mpcalc.cross_totals(p, T, Td)
vt = mpcalc.vertical_totals(p, T)
sw = mpcalc.sweat_index(p, T, Td, u, v)
ssi = mpcalc.showalter_index(p, T, Td)
li = mpcalc.lifted_index(p, T, prof)
# Wind / shear
u, v = mpcalc.wind_components(speed, direction)
shear = mpcalc.bulk_shear(p, u, v, height=z, depth=6*units.km)
srh, _, _ = mpcalc.storm_relative_helicity(z, u, v, depth=3*units.km)
rm, lm, _ = mpcalc.bunkers_storm_motion(p, u, v, z)
mcs_up, mcs_dn = mpcalc.corfidi_storm_motion(p, u, v)
# Composite parameters
stp = mpcalc.significant_tornado(sbcape, lcl_height, srh, shear)
scp = mpcalc.supercell_composite(mucape, srh, shear)
Requires grid spacing (use mpcalc.lat_lon_grid_deltas for lat/lon grids):
dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)
# Wind
spd = mpcalc.wind_speed(u, v)
dirn = mpcalc.wind_direction(u, v)
# Vorticity / divergence
zeta = mpcalc.vorticity(u, v, dx=dx, dy=dy)
abvort = mpcalc.absolute_vorticity(u, v, dx=dx, dy=dy, latitude=lat)
div = mpcalc.divergence(u, v, dx=dx, dy=dy)
# Deformation
shear_def = mpcalc.shearing_deformation(u, v, dx=dx, dy=dy)
stretch_def = mpcalc.stretching_deformation(u, v, dx=dx, dy=dy)
total_def = mpcalc.total_deformation(u, v, dx=dx, dy=dy)
# Geostrophic / ageostrophic wind
ug, vg = mpcalc.geostrophic_wind(height, dx=dx, dy=dy, latitude=lat)
uag, vag = mpcalc.ageostrophic_wind(height, u, v, dx=dx, dy=dy)
# Advection
T_adv = mpcalc.advection(T, u=u, v=v, dx=dx, dy=dy)
# Frontogenesis
fronto = mpcalc.frontogenesis(theta, u, v, dx=dx, dy=dy)
# Potential vorticity
pv = mpcalc.potential_vorticity_baroclinic(theta, u, v, p)
pv = mpcalc.potential_vorticity_barotropic(height, u, v, dx=dx, dy=dy)
# Q-vectors
qx, qy = mpcalc.q_vector(u, v, T, p, dx=dx, dy=dy)
# Coriolis
f = mpcalc.coriolis_parameter(lat)
# Brunt-Vaisala frequency
N = mpcalc.brunt_vaisala_frequency(z, T, p)
N2 = mpcalc.brunt_vaisala_frequency_squared(z, T, p)
# Richardson number
Ri = mpcalc.gradient_richardson_number(z, T, u, v)
# TKE / friction velocity (from time series)
tke = mpcalc.tke(u, v, w, perturbation=True)
ustar = mpcalc.friction_velocity(u, w)
hi = mpcalc.heat_index(T, rh)
wc = mpcalc.windchill(T, speed)
at = mpcalc.apparent_temperature(T, rh, speed)
p_std = mpcalc.height_to_pressure_std(z)
z_std = mpcalc.pressure_to_height_std(p)
p_sl = mpcalc.altimeter_to_sea_level_pressure(altimeter, z)
T_smooth = mpcalc.smooth_gaussian(T_2d, n=3)
T_smooth = mpcalc.smooth_n_point(T_2d, n=9, passes=1)
T_smooth = mpcalc.smooth_rectangular(T_2d, size=(3,3), passes=2)
T_smooth = mpcalc.smooth_circular(T_2d, radius=2, passes=1)
theta_levels = np.array([300, 310, 320]) * units.K
ds_isentropic = mpcalc.isentropic_interpolation_as_dataset(
theta_levels, T, p, u, v
)
# Grid spacing from lat/lon
dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)
# Layer extraction
p_layer, T_layer = mpcalc.get_layer(p, T, height=z, bottom=z[0], depth=3*units.km)
# Find intersections (e.g. parcel/environment profile)
x_int, y_int = mpcalc.find_intersections(p, T, prof)
# Perturbation from mean
T_prime = mpcalc.get_perturbation(T_ts)