from pylab import *
import matplotlib.tri as Tri
import netCDF4
import datetime as dt

# DAP Data URL
url = 'http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_FVCOM_OCEAN_MASSBAY_FORECAST.nc'
# Open DAP
nc = netCDF4.Dataset(url)
nc.variables.keys()

# take a look at the "metadata" for the variable "u"
print nc.variables['u']

# Desired time for snapshot
# ....right now (or some number of hours from now) ...
start = dt.datetime.utcnow() + dt.timedelta(hours=0)
# ... or specific time (UTC)
start = dt.datetime(2013,3,2,15,0,0)

# Get desired time step  
time_var = nc.variables['time']
itime = netCDF4.date2index(start,time_var,select='nearest')

# Get lon,lat coordinates for nodes (depth)
lat = nc.variables['lat'][:]
lon = nc.variables['lon'][:]
# Get lon,lat coordinates for cell centers (depth)
latc = nc.variables['latc'][:]
lonc = nc.variables['lonc'][:]
# Get Connectivity array
nv = nc.variables['nv'][:].T - 1 
# Get depth
h = nc.variables['h'][:]  # depth 

dtime = netCDF4.num2date(time_var[itime],time_var.units)
daystr = dtime.strftime('%Y-%b-%d %H:%M')
print daystr

tri = Tri.Triangulation(lon,lat, triangles=nv)

# get current at layer [0 = surface, -1 = bottom]
ilayer = 0
u = nc.variables['u'][itime, ilayer, :]
v = nc.variables['v'][itime, ilayer, :]

levels=arange(-32,2,1)   # depth contours to plot
ax= [-70.97, -70.82, 42.25, 42.35] # region to plot

# find velocity points in bounding box
ind = argwhere((lonc >= ax[0]) & (lonc <= ax[1]) & (latc >= ax[2]) & (latc <= ax[3]))

subsample=3
np.random.shuffle(ind)
Nvec = int(len(ind) / subsample)
idv = ind[:Nvec]

# tricontourf plot of water depth with vectors on top
fig1 = figure(figsize=(18,10))
ax1 = fig1.add_subplot(111,aspect=(1.0/cos(mean(lat)*pi/180.0)))
plt.tricontourf(tri, -h,levels=levels,shading='faceted',cmap=plt.cm.gist_earth)
plt.axis(ax)
ax1.patch.set_facecolor('0.5')
cbar=colorbar()
cbar.set_label('Water Depth (m)', rotation=-90)
Q = ax1.quiver(lonc[idv],latc[idv],u[idv],v[idv],scale=20)
qk = quiverkey(Q,0.92,0.08,0.50,'0.5 m/s',labelpos='W')
title('NECOFS Velocity, Layer %d, %s' % (ilayer, daystr))