The EDEX modelsounding plugin creates 64-level vertical profiles from GFS and ETA (NAM) BUFR products distirubted over NOAAport. Paramters which are requestable are pressure, temperature, specHum, uComp, vComp, omega, cldCvr.
from awips.dataaccess import DataAccessLayer
import matplotlib.tri as mtri
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
from math import exp, log
import numpy as np
from metpy.calc import get_wind_components, lcl, dry_lapse, parcel_profile, dewpoint
from metpy.calc import get_wind_speed,get_wind_dir, thermo, vapor_pressure
from metpy.plots import SkewT, Hodograph
from metpy.units import units, concatenate
DataAccessLayer.changeEDEXHost("edex-cloud.unidata.ucar.edu")
request = DataAccessLayer.newDataRequest()
request.setDatatype("modelsounding")
forecastModel = "GFS"
request.addIdentifier("reportType", forecastModel)
request.setParameters("pressure","temperature","specHum","uComp","vComp","omega","cldCvr")
locations = DataAccessLayer.getAvailableLocationNames(request)
locations.sort()
list(locations)
['KFRM', 'KFSD', 'KMHE', 'KYKN', 'OAX']
request.setLocationNames("KFRM")
cycles = DataAccessLayer.getAvailableTimes(request, True)
times = DataAccessLayer.getAvailableTimes(request)
try:
fcstRun = DataAccessLayer.getForecastRun(cycles[-1], times)
list(fcstRun)
response = DataAccessLayer.getGeometryData(request,[fcstRun[0]])
except:
print('No times available')
exit
Construct arrays for each parameter to plot (temperature, pressure, moisutre (spec. humidity), wind components, and cloud cover)
tmp,prs,sh = np.array([]),np.array([]),np.array([])
uc,vc,om,cld = np.array([]),np.array([]),np.array([]),np.array([])
for ob in response:
tmp = np.append(tmp,ob.getNumber("temperature"))
prs = np.append(prs,ob.getNumber("pressure"))
sh = np.append(sh,ob.getNumber("specHum"))
uc = np.append(uc,ob.getNumber("uComp"))
vc = np.append(vc,ob.getNumber("vComp"))
om = np.append(om,ob.getNumber("omega"))
cld = np.append(cld,ob.getNumber("cldCvr"))
print("parms = " + str(ob.getParameters()))
print("site = " + str(ob.getLocationName()))
print("geom = " + str(ob.getGeometry()))
print("datetime = " + str(ob.getDataTime()))
print("reftime = " + str(ob.getDataTime().getRefTime()))
print("fcstHour = " + str(ob.getDataTime().getFcstTime()))
print("period = " + str(ob.getDataTime().getValidPeriod()))
sounding_title = forecastModel + " " + str(ob.getLocationName()) + "("+ str(ob.getGeometry())+")" + str(ob.getDataTime())
parms = ['uComp', 'cldCvr', 'temperature', 'vComp', 'pressure', 'omega', 'specHum'] site = KFRM geom = POINT (-94.41999816894531 43.65000152587891) datetime = 1970-01-18 13:45:50.400000 (0) reftime = Jan 18 70 13:45:50 GMT fcstHour = 15 period = (Jan 18 70 13:46:05 , Jan 18 70 13:46:05 )
Because the modelsounding plugin does not return dewpoint values, we must calculate the profile ourselves. Here are three examples of dewpoint calculated from specific humidity, including a manual calculation following NCEP AWIPS/NSHARP.
1) MetPy calculated mixing ratio and vapor pressure
t = (tmp-273.15) * units.degC
p = prs/100 * units.mbar
u,v = uc*1.94384,vc*1.94384 # m/s to knots
spd = get_wind_speed(u, v) * units.knots
dir = get_wind_dir(u, v) * units.deg
rmix = (sh/(1-sh)) *1000 * units('g/kg')
e = vapor_pressure(p, rmix)
td = dewpoint(e)
2) metpy calculated assuming spec. humidity = mixing ratio
td2 = dewpoint(vapor_pressure(p, sh))
3) NCEP AWIPS soundingrequest plugin
based on GEMPAK/NSHARP, from https://github.com/Unidata/awips2-ncep/blob/unidata_16.2.2/edex/gov.noaa.nws.ncep.edex.plugin.soundingrequest/src/gov/noaa/nws/ncep/edex/plugin/soundingrequest/handler/MergeSounding.java#L1783
# new arrays
ntmp = tmp
# where p=pressure(pa), T=temp(C), T0=reference temp(273.16)
rh = 0.263*prs*sh / (np.exp(17.67*ntmp/(ntmp+273.15-29.65)))
vaps = 6.112 * np.exp((17.67 * ntmp) / (ntmp + 243.5))
vapr = rh * vaps / 100
dwpc = np.array(243.5 * (np.log(6.112) - np.log(vapr)) / (np.log(vapr) - np.log(6.112) - 17.67)) * units.degC
/Users/mjames/miniconda2/envs/python-awips/lib/python2.7/site-packages/ipykernel_launcher.py:8: RuntimeWarning: divide by zero encountered in log /Users/mjames/miniconda2/envs/python-awips/lib/python2.7/site-packages/ipykernel_launcher.py:8: RuntimeWarning: invalid value encountered in divide
%matplotlib inline
plt.rcParams['figure.figsize'] = (12, 14)
# Create a skewT plot
skew = SkewT()
# Plot the data
skew.plot(p, t, 'r', linewidth=2)
skew.plot(p, td, 'b', linewidth=2)
skew.plot(p, td2, 'y')
skew.plot(p, dwpc, 'g', linewidth=2)
skew.plot_barbs(p, u, v)
skew.ax.set_ylim(1000, 100)
skew.ax.set_xlim(-40, 60)
plt.title(sounding_title)
# An example of a slanted line at constant T -- in this case the 0 isotherm
l = skew.ax.axvline(0, color='c', linestyle='--', linewidth=2)
# Draw hodograph
ax_hod = inset_axes(skew.ax, '40%', '40%', loc=2)
h = Hodograph(ax_hod, component_range=get_wind_speed(u, v).max())
h.add_grid(increment=20)
h.plot_colormapped(u, v, spd)
# Show the plot
plt.show()