On Tue, Jul 12, 2016 at 11:21 AM, Bryan Guarente <guarente@xxxxxxxx> wrote: > I have worked through my own version of your sounding and advanced > sounding examples using python-awips and metpy. However, I never saw any > ability to make a model sounding. Is there an example of making a model > sounding with MetPy and python-awips? Initially, I think this would be > simple, but then I realized we would be using a completely different > filetype and thus would likely have to use a different library to decode > it. > > If you have any examples, can you post them to github? > Bryan, Michael James (cc'ed, maintainer of python-awips) came up with the attached and had this to add: "I believe this is correct, assuming specific humidity is equal to vapor pressure e (modelsounding doesn't return dewpoint). Also note that the available locations are limited to areas around OAX due to localization. I'll have to update this to allow global access in the next release." I'm currently working on a gallery of notebooks that work across all of the meteorology tools, and I'll add these when that's ready (I don't want MetPy's own examples pulling from external, online data sources.) Ryan -- Ryan May, Ph.D. Software Engineer UCAR/Unidata Boulder, CO
Attachment:
Model_Soundings.ipynb
Description: Binary data
# coding: utf-8
# In[156]:
from awips.dataaccess import DataAccessLayer
DataAccessLayer.changeEDEXHost("edex-cloud.unidata.ucar.edu")
request = DataAccessLayer.newDataRequest()
request.setDatatype("modelsounding")
# In[157]:
availableLocs = DataAccessLayer.getAvailableLocationNames(request)
availableLocs.sort()
for loc in availableLocs: print loc
# In[158]:
request.addIdentifier("reportType", "ETA")
request.setLocationNames("OAX")
# In[159]:
availableParms = DataAccessLayer.getAvailableParameters(request)
availableParms.sort()
for parm in availableParms: print parm
# In[160]:
request.setParameters("pressure","temperature","specHum","uComp","vComp")
# In[161]:
cycles = DataAccessLayer.getAvailableTimes(request, True)
print "using ", str(cycles[-1]) # 0 for FIRST time, -1 for LAST
# In[162]:
allTimes = DataAccessLayer.getAvailableTimes(request)
# Build one complete model run
fcstRun = []
for time in allTimes:
if str(time)[:19] == str(cycles[-1]):
fcstRun.append(time)
#for time in fcstRun: print time
# In[163]:
response = DataAccessLayer.getGeometryData(request,times=[fcstRun[0]])
# In[164]:
print "parms = " + str(response[0].getParameters())
print "site = " + response[0].getLocationName()
print "datetime = " + str(response[0].getDataTime())
print "geom = " + str(response[0].getGeometry())
print ""
print fcstRun[0]
tmp,prs,sh,uc,vc = [],[],[],[],[]
for ob in response:
#print dir(ob)
print ob.getDataTime()
tmp.append(float(ob.getString("temperature")))
prs.append(float(ob.getString("pressure")))
print float(ob.getString("pressure"))
sh.append(float(ob.getString("specHum")))
uc.append(float(ob.getString("uComp")))
vc.append(float(ob.getString("vComp")))
# In[165]:
import matplotlib.tri as mtri
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
from math import exp
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
# we can use units.* here...
t = (np.array(tmp)-273.15) * units.degC
p = (np.array(prs)/100) * units.mbar
s = np.array(sh)
u, v = np.array(uc),np.array(vc)
spd = get_wind_speed(u, v) * units.knot
dir = get_wind_dir(u, v) * units.deg
# assuming spec. humidity equals vapor pressure (e)
td = dewpoint(vapor_pressure(p, s))
print min(p), max(p)
print min(t), max(t)
print min(s), max(s)
print min(spd), max(spd)
print min(dir), max(dir)
print min(td), max(td)
# In[167]:
get_ipython().magic(u'matplotlib inline')
plt.rcParams['figure.figsize'] = (12, 14)
# Create a skewT plot
skew = SkewT()
# 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, 'r')
skew.plot(p, td, 'g')
skew.plot_barbs(p, u, v)
skew.ax.set_ylim(1000, 100)
skew.ax.set_xlim(-40, 60)
# Calculate LCL height and plot as black dot
l = lcl(p[0], t[0], td[0])
lcl_temp = dry_lapse(concatenate((p[0], l)), t[0])[-1].to('degC')
skew.plot(l, lcl_temp, 'ko', markerfacecolor='black')
# Calculate full parcel profile and add to plot as black line
#prof = parcel_profile(p, t[0], td[0]).to('degC')
#skew.plot(p, prof, 'k', linewidth=2)
# Example of coloring area between profiles
#skew.ax.fill_betweenx(p, t, prof, where=t>=prof, facecolor='blue', alpha=0.4)
#skew.ax.fill_betweenx(p, t, prof, where=t<prof, facecolor='red', alpha=0.4)
# 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=3)
h = Hodograph(ax_hod, component_range=80.)
h.add_grid(increment=20)
h.plot_colormapped(u, v, spd)
# Show the plot
plt.show()
# In[ ]:
python-users archives: