In [2]:

    

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


    

In [3]:

    

# DAP Data URL
url = 'http://testbedapps-dev.sura.org/thredds/dodsC/in/usf/fvcom/ike/ultralite/vardrag/nowave/3d'


    

In [4]:

    

# open with NetCDF4
nc = netCDF4.Dataset(url).variables
nc.keys()


    

    
Out[4]:





[u'time',
 u'Times',
 u'lon',
 u'lat',
 u'lonc',
 u'latc',
 u'siglay',
 u'h',
 u'nv',
 u'zeta',
 u'u',
 u'v',
 u'ua',
 u'va',
 u'maxele',
 u'fvcom_mesh']

In [5]:

    

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


    

    




<type 'netCDF4.Variable'>
float32 u(u'time', u'siglay', u'nele')
    long_name: Eastward water velocity
    units: meters s-1
    type: data
    missing_value: -999.0
    field: u, scalar
    standard_name: eastward_sea_water_velocity
    coordinates: time siglay latc lonc
    mesh: fvcom_mesh
    location: face
    coverage_content_type: modelResult
unlimited dimensions = (u'time',)
current size = (361, 10, 826866)

In [6]:

    

# Desired time for snapshot
# ....right now (or some number of hours from now) ...
start = dt.datetime.utcnow() + dt.timedelta(hours=33)
# ... or specific time (UTC)
start = dt.datetime(2008,9,13,6,0,0)


    

In [7]:

    

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


# Get depth
h = nc['h'][:]  # depth


    

In [8]:

    

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


    

    




2008-Sep-13 06:00

In [9]:

    

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


    

In [10]:

    

meshvar = nc[nc['u'].mesh]


    

In [11]:

    

print meshvar


    

    




<type 'netCDF4.Variable'>
int32 fvcom_mesh()
    cf_role: mesh_topology
    topology_dimension: 2
    node_coordinates: lon lat
    face_coordinates: lonc latc
    face_node_connectivity: nv
unlimited dimensions = ()
current size = ()

In [12]:

    

# Get lon,lat coordinates for nodes (depth)
lat = nc['lat'][:]
lon = nc['lon'][:]
# Get lon,lat coordinates for cell centers (depth)
latc = nc['latc'][:]
lonc = nc['lonc'][:]


    

In [13]:

    

trivar=nc[meshvar.face_node_connectivity]


    

In [14]:

    

print trivar


    

    




<type 'netCDF4.Variable'>
float32 nv(u'three', u'nele')
    long_name: nodes surrounding element
    cf_role: face_node_connectivity
    start_index: 1
unlimited dimensions = ()
current size = (3, 826866)

In [15]:

    

nv = trivar[:] - trivar.start_index
[m,n]=shape(nv)
if (m < n):
    nv = nv.T


    

In [16]:

    

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


    

In [17]:

    

#woods hole
levels=arange(-30,2,1)
ax = [-70.7, -70.6, 41.48, 41.55]
maxvel = 1.0
subsample = 2


    

In [18]:

    

#boston harbor
levels=arange(-34,2,1)   # depth contours to plot
ax= [-70.97, -70.82, 42.25, 42.35] # 
maxvel = 0.5
subsample = 3


    

In [19]:

    

#galveston bay
levels=arange(-34,2,1)   # depth contours to plot
ax=[-95.2519,  -94.3856,   29.0,   30.0]
maxvel = 6.0
subsample = 6


    

In [20]:

    

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


    

In [21]:

    

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


    

In [22]:

    

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


    

In [22]: