Design philosophy:
pint
)pint
, which provides a unit registry of all the units it knows about.Quantity
object that wraps numpy
arrays to provide unit-aware math operations.metpy.constants
module .# Simple example. The units name imported here is the registry.
from metpy.units import units
import metpy.constants as consts
dist = 10 * units.miles
time = 15 * units.minutes
avg_speed = dist / time
avg_speed
# Convert to something more meaningful to me
avg_speed.to(units.mph)
Can also give pint
a unit expression that it parses:
avg_speed.to('m/s')
Could also do something a bit more useful for physics:
m = 102 * units.kg
a = 9.81 * units('m/s^2')
F = m * a
F.to('N')
Variety of ways to add units:
3e8 * units.meter / units.second
3e8 * units('m/s')
units.Quantity(3e8, units.m / units.sec)
units.Quantity(3e8, 'm/s')
One of the factors that led to the choice of pint
was the fact that it supports temperature units, both changes in temperature and absolute temperature values (which are not strictly multiplicative).
T = 30 * units.degF
P = 500 * units.mbar
P0 = 1000 * units.mbar
T * (P0 / P) ** consts.kappa
T * 2
(T * (1000 / 500) ** consts.kappa).to('K')
What's the density of air at 925 mb and 82 F?
$$P = \rho R T$$$R$ for dry air is in metpy.constants
as Rd
.
Bigger problem: How much energy does it take to raise the temperature of liquid water:
$$Q = c_p m \Delta T$$Assume: 60 kg of water, from 32F to 85F. $C_p$ for liquid water is in metpy.constants
as Cp_l
.
It's 86F outside here (pressure 860 mb), with a relative humidity of 40%. What's the dewpoint?
import metpy.calc as mcalc
temp = 86 * units.degF
press = 860. * units.mbar
humidity = 40 / 100.
dewpt = mcalc.dewpoint_rh(temp, humidity).to('degF')
dewpt
What's the water vapor mixing ratio for those conditions?
e = mcalc.saturation_vapor_pressure(dewpt)
mcalc.mixing_ratio(e, press).to('g/kg')
MetPy also has more sounding-focused calculations. For instance, what's the lifted condensation level (LCL) for those conditions?
mcalc.lcl(press, temp, dewpt)
Or we can even calculate the parcel profile (temperature at a series of pressure levels) for those conditions:
import numpy as np
pressure_levels = np.array([860., 850., 700., 500., 300.]) * units.mbar
mcalc.parcel_profile(pressure_levels, temp, dewpt)
It's 95F outside with a dewpoint of 70F. What's the heat index?
Note: $$RH = e / e_s$$
axvline
directlyfrom siphon.catalog import TDSCatalog
best_gfs = TDSCatalog('http://thredds.ucar.edu/thredds/catalog/grib/NCEP/GFS/Global_0p5deg/catalog.xml?dataset=grib/NCEP/GFS/Global_0p5deg/Best')
best_ds = list(best_gfs.datasets.values())[0]
# Set up a class to access
from siphon.ncss import NCSS
ncss = NCSS(best_ds.access_urls['NetcdfSubset'])
# Get today's date
from datetime import datetime, timedelta
now = datetime.utcnow()
now = datetime(2015, 7, 18, 12, 0, 0)
# Re-use the NCSS access from earlier, but make a new query
query = ncss.query()
query.lonlat_point(-105, 40).time(now + timedelta(hours=12)).accept('csv')
query.variables('Temperature_isobaric', 'Relative_humidity_isobaric',
'u-component_of_wind_isobaric', 'v-component_of_wind_isobaric')
data = ncss.get_data(query)
# Pull out data with some units. Reversing arrays to put higher pressure first.
p = (data['vertCoord'][::-1] * units('Pa')).to('mbar')
T = data['Temperature_isobaric'][::-1] * units('kelvin')
Td = mcalc.dewpoint_rh(T, data['Relative_humidity_isobaric'][::-1] / 100.)
u = data['ucomponent_of_wind_isobaric'][::-1] * units('m/s')
v = data['vcomponent_of_wind_isobaric'][::-1] * units('m/s')
# Create a skewT using matplotlib's default figure size
import matplotlib.pyplot as plt
%matplotlib inline
fig = plt.figure(figsize=(8, 8))
from metpy.plots import SkewT
skew = SkewT(fig)
# Plot the data using normal plotting functions, in this case using
# log scaling in Y, as dictated by the typical meteorological plot
skew.plot(p, T.to('degC'), 'r')
skew.plot(p, Td.to('degC'), 'g')
skew.plot_barbs(p, u.to('knots'), v.to('knots'))
skew.ax.set_ylim(1000, 100)
skew.ax.set_title(data['date'][0]);
Nice, but we could use some of the typical lines:
# Add the relevant special lines
skew.plot_dry_adiabats()
skew.plot_moist_adiabats()
skew.plot_mixing_lines()
fig
Can also use "vertical" lines to highlight isotherms.
skew.ax.axvline(0)
skew.ax.axvline(-20)
fig
Or plot the parcel profile on the Skew-T
prof = mcalc.parcel_profile(p, T[0], Td[0]).to('degC')
skew.plot(p, prof, 'k', linewidth=2)
fig
And can use matplotlib to fill the positive and negative areas.
skew.ax.fill_betweenx(p, T.to('degC'), prof, where=T>=prof, facecolor='blue', alpha=0.4)
skew.ax.fill_betweenx(p, T.to('degC'), prof, where=T<prof, facecolor='red', alpha=0.4)
fig
Plot a different sounding from the HRRR:
from siphon.catalog import TDSCatalog
best_gfs = TDSCatalog('http://thredds-jumbo.unidata.ucar.edu/thredds/catalog/grib/HRRR/CONUS_3km/wrfprs/catalog.xml?dataset=grib/HRRR/CONUS_3km/wrfprs/Best')
best_ds = list(best_gfs.datasets.values())[0]
from siphon.ncss import NCSS
ncss = NCSS(best_ds.access_urls['NetcdfSubset'])
query = ncss.query()
query.lonlat_point(-105, 40).time(now + timedelta(hours=12)).accept('csv')
query.variables('Temperature_isobaric', 'Relative_humidity_isobaric',
'u-component_of_wind_isobaric', 'v-component_of_wind_isobaric')
data = ncss.get_data(query)
# Pull out data with some units. Reversing arrays to put higher pressure first.
p = (data['vertCoord'][::-1] * units('Pa')).to('mbar')
T = data['Temperature_isobaric'][::-1] * units('kelvin')
Td = mcalc.dewpoint_rh(T, data['Relative_humidity_isobaric'][::-1] / 100.)
u = data['ucomponent_of_wind_isobaric'][::-1] * units('m/s')
v = data['vcomponent_of_wind_isobaric'][::-1] * units('m/s')
# Plot the sounding here
# Add some of the standard lines
# Highlight the -10C isotherm with a blue dashed line
# If you're really feeling ambitious, plot the parcel profile, or the LCL.
from IPython.display import Image
Image('images/station-plot.png')
Image('images/hodograph.png')
Image('images/radar-plots.png')
MetPy is open-source and all development done in the open:
We want to encourage everyone to join us in making MetPy as useful as possible: