# Make the output of plotting commands be displayed inline within the notebook,
%matplotlib inline 
from mpl_toolkits.basemap import Basemap  # import Basemap matplotlib toolkit
import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset  # netcdf4-python module

gfs_fcst = Dataset('http://nomads.ncep.noaa.gov:9090/dods/gens/gens20140123/gep_all_00z')
#gfs_fcst = Dataset('basemap_netcdf4-python.nc') # use local netcdf file instead
print gfs_fcst # get some summary information about the dataset

time = gfs_fcst.variables['time']
print time
from netCDF4 import num2date
valid_dates = num2date(time[:], time.units).tolist()
print [d.strftime('%Y%m%d%H') for d in valid_dates]

from datetime import datetime
nt_superbowl = valid_dates.index(datetime(2014, 2, 3, 0, 0))
print valid_dates[nt_superbowl]

t2mens = gfs_fcst.variables['tmp2m'][:,nt_superbowl,:,:]
t2m = t2mens.mean(axis=0)
print t2m.shape, t2m.min(), t2m.max()
lats = gfs_fcst.variables['lat'][:]; lons = gfs_fcst.variables['lon'][:]
lons, lats = np.meshgrid(lons, lats)  # make lats/lons into 2D arrays (needed for plotting)

fig = plt.figure(figsize=(8,8))
m = Basemap(projection='stere',lon_0=-74,lat_0=41,width=4.e6,height=4.e6)
m.drawcoastlines()
m.drawstates()
m.drawcountries()
x,y = m(lons, lats) # convert lats/lons to map projection coordinates
cs = m.contourf(x,y,t2m,np.linspace(230,300,41),cmap=plt.cm.jet)  # color-filled contours
cb = m.colorbar(cs)  # draw colorbar
parallels = m.drawparallels(np.arange(10,90,10),labels=[1,0,0,0])  # draw parallels, label on left
meridians = m.drawmeridians(np.arange(-10,-100,-10),labels=[0,0,0,1]) # label meridians on bottom
t = plt.title('ens mean t2m fcst from %s for %s' % (valid_dates[0],valid_dates[nt_superbowl]),fontweight='bold')

fig = plt.figure(figsize=(16,14))
for nens in range(1,21):
    ax = plt.subplot(4,5,nens)
    m.drawcoastlines()
    cs = m.contourf(x,y,t2mens[nens],np.linspace(230,300,41),cmap=plt.cm.jet)
    t = plt.title('ens member %s' % nens)

freezeprob = 100.*((t2mens < 273).sum(axis=0))/t2mens.shape[0]
print freezeprob.min(), freezeprob.max(), freezeprob.shape
m.drawcoastlines(color='r')
from mpl_toolkits.basemap import cm as basemapcm
cs = m.contourf(x,y,freezeprob,np.linspace(0,100.,41),cmap=basemapcm.GMT_haxby_r)
cb = m.colorbar()
xpt, ypt = m(-74,41)
nyc = m.plot([xpt],[ypt],'ro') # plot red dot near the Meadowlands.
t = plt.title('ens prob t2m < 273 (go Broncos!)')