In [1]:
from pylab import *
from matplotlib.collections import PolyCollection
from matplotlib.collections import TriMesh
import matplotlib.tri as Tri
from mpl_toolkits.basemap import Basemap
import datetime as dt
import netCDF4
In [2]:
import matplotlib
matplotlib.__version__
Out[2]:
In [3]:
url = 'http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_FVCOM_OCEAN_MASSBAY_FORECAST.nc'
#url = 'http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_GOM2_FORECAST.nc'
nc = netCDF4.Dataset(url)
# read node locations
lat = nc.variables['lat'][:]
lon = nc.variables['lon'][:]
# read element centroid locations
latc = nc.variables['latc'][:]
lonc = nc.variables['lonc'][:]
# read connectivity array
nv = nc.variables['nv'][:].T - 1
time_var = nc.variables['time']
In [4]:
# create a triangulation object, specifying the triangle connectivity array
tri = Tri.Triangulation(lon,lat, triangles=nv)
In [5]:
# plot depth using tricontourf
h = nc.variables['h'][:]
fig=figure(figsize=(12,12))
ax=fig.add_subplot(111,aspect=1.0/cos(latc.mean() * pi / 180.0))
tricontourf(tri,-h,levels=range(-300,10,10))
colorbar()
Out[5]:
In [6]:
# get velocity nearest to current time
start = dt.datetime.utcnow()+ dt.timedelta(hours=0)
istart = netCDF4.date2index(start,time_var,select='nearest')
layer = 0 # surface layer
u = nc.variables['u'][istart, layer, :]
v = nc.variables['v'][istart, layer, :]
mag = numpy.sqrt((u*u)+(v*v))
In [7]:
# Now try plotting speed and vectors with Basemap using a PolyCollection
m = Basemap(projection='merc', llcrnrlat=lat.min(), urcrnrlat=lat.max(),
llcrnrlon=lon.min(), urcrnrlon=lon.max(), lat_ts=lat.mean(), resolution=None)
# project from lon,lat to mercator
xnode, ynode = m(lon, lat)
xc, yc = m(lonc, latc)
# create a TRI object with projected coordinates
tri = Tri.Triangulation(xnode, ynode, triangles=nv)
# make a PolyCollection using triangles
verts = concatenate((tri.x[tri.triangles][..., None],
tri.y[tri.triangles][..., None]), axis=2)
collection = PolyCollection(verts)
collection.set_edgecolor('none')
In [8]:
timestamp=start.strftime('%Y-%m-%d %H:%M:%S')
In [9]:
# set the magnitude of the polycollection to the speed
collection.set_array(mag)
collection.norm.vmin=0
collection.norm.vmax=0.5
In [10]:
fig=figure(figsize=(12,12))
ax=fig.add_subplot(111)
m.drawmapboundary(fill_color='0.3')
#m.drawcoastlines()
#m.fillcontinents()
# add the speed as colored triangles
ax.add_collection(collection) # add polygons to axes on basemap instance
# add the vectors
Q = m.quiver(xc,yc,u,v,scale=30)
# add a key for the vectors
qk = plt.quiverkey(Q,0.1,0.1,0.20,'0.2 m/s',labelpos='W')
title('FVCOM Surface Current speed at %s UTC' % timestamp)
Out[10]:
In [11]:
# try using the TriMesh collection: can't figure this out
collection2 = TriMesh(tri)
fig=figure(figsize=(12,12))
ax = fig.add_subplot(111)
ax.add_collection(collection2)
Out[11]: