In [2]:

    

# Extract subsetted snapshots of wind from 4km NAM model

# This approach would work for any CF-Compliant NetCDF or OPeNDAP dataset


    

In [1]:

    

from pylab import *
import netCDF4
import datetime as dt
from mpl_toolkits.basemap import Basemap


    

In [43]:

    

# wind now
dnow = dt.datetime.utcnow() - dt.timedelta(hours=1500)
# or wind at specific time (UTC)
dnow = dt.datetime(2013,6,13,21)
daystr = dnow.strftime('%Y%m%d')
print daystr


    

    




20130613

In [44]:

    

cycle=max(0,(dnow.hour/6-1)*6)  # rounding by 6 hours to last cycle
print cycle


    

    




12

In [45]:

    

#NAM 4km
url='http://nomads.ncep.noaa.gov:9090/dods/nam/nam%s/nam_conusnest_%2.2dz' % (daystr,cycle)
url='http://nomads.ncep.noaa.gov:9090/dods/nam/nam%s/nam_conusnest_%2.2dz' % ('20130613',6)

lon_range=[-74.5, -67.];lat_range=[39., 44.];  # 4km lon is [-180 180]

# NAM 12km
#url=sprintf('http://nomads.ncep.noaa.gov:9090/dods/nam/nam%s/nam_%2.2dz',str,d,h);
#lon_range=[-71.5 -63];lat_range=[41 46];  % 12km lon is [-180 180]

# NAM 32km
#url=sprintf('http://nomads.ncep.noaa.gov:9090/dods/nam/nam%s/nam_na_%2.2dz',str,d,h);
#lon_range=[-71.5 -63]+360;lat_range=[41 46]; % 32km lon is [0 360]

print url
bbox=[lon_range[0], lat_range[0], lon_range[1], lat_range[1]]
print bbox


    

    




http://nomads.ncep.noaa.gov:9090/dods/nam/nam20130613/nam_conusnest_06z
[-74.5, 39.0, -67.0, 44.0]

In [46]:

    

# NetCDF4-Python can open OPeNDAP dataset just like a local NetCDF file
nc = netCDF4.Dataset(url)


    

In [47]:

    

#nc.variables.keys()


    

In [48]:

    

uvar = nc.variables['ugrd10m']
vvar = nc.variables['vgrd10m']
pvar = nc.variables['pressfc']


    

In [49]:

    

def nam_subset(nc,bbox):
    lat = nc.variables['lat'][:]
    lon = nc.variables['lon'][:]
#    lon = lon - 360.0 
    ibox =(lon>=bbox[0])&(lon<=bbox[1])
    jbox =(lat>=bbox[2])&(lat<=bbox[3])
    return lon[ibox],lat[jbox],ibox,jbox


    

In [50]:

    

lon,lat,ibox,jbox = nam_subset(nc,bbox)

time_var = nc.variables['time']


    

In [51]:

    

# Extract desired time.  
start = dnow + dt.timedelta(hours=0)
istart = netCDF4.date2index(start,time_var,select='nearest')
dtime = netCDF4.num2date(time_var[istart],time_var.units)
print start.strftime('%Y/%m/%d %H:%M UTC')
print dtime.strftime('%Y/%m/%d %H:%M UTC')


    

    




2013/06/13 21:00 UTC
2013/06/13 21:00 UTC

In [52]:

    

# Get snapshot of variable [vname] at time = istart
u = uvar[istart,jbox,ibox]
v = vvar[istart,jbox,ibox]
p = pvar[istart,jbox,ibox]/1000.   # kPa
spd = sqrt(u*u + v*v)
#hs1 = ma.masked_where(isnan(hs1),hs1)
# tim = dtime[istart:]


    

In [53]:

    

nsub=4 # subsample vectors
scale=0.002 # length of vectors


    

In [54]:

    

fig=figure(figsize=(12,12)) 
m = Basemap(resolution='i', llcrnrlon=bbox[0],llcrnrlat=bbox[1],
            urcrnrlon=bbox[2],urcrnrlat=bbox[3],lat_0=bbox[1],lon_0=bbox[0],lat_ts=bbox[1])
m.drawcoastlines()
#ax=fig.add_subplot(111,aspect=1.0/cos(min(lat) * pi / 180.0))
pc3=pcolormesh(lon,lat,p,vmin=99.,vmax=102.)
quiver(lon[::nsub],lat[::nsub],u[::nsub,::nsub],v[::nsub,::nsub],scale=1.0/scale, zorder=1e35, width=0.002)
title('NAM Wind Speed in the Gulf of Maine: ' + dtime.strftime('%Y/%m/%d %H:%M UTC'))
#parallels = np.arange(-75.,-63,5.)
#    meridians = np.arange(-100.,-50.,5.)
#    m.drawparallels(parallels,labels=[1,0,0,0],fontsize=10);
#    m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10);
cbax3 = fig.add_axes([0.95, 0.3, 0.03, 0.4])
cb3 = plt.colorbar(pc3, cax=cbax3,  orientation='vertical')
cb3.set_label('kPa')
plt.savefig(dtime.strftime('%Y%m%dT%H%M00Z')+'.png')


    

In [17]: