In [4]:

    

from pylab import *
import collections

import matplotlib.tri as Tri
from shapely import geometry
import wkb2shp
import netCDF4
import datetime
import time


    

In [5]:

    

titl='ADCIRC';
url='http://testbedapps-dev.sura.org/thredds/dodsC/auto/in.und.adcirc.ike.ultralite.lr.vardrag.wave.2d'
vname='zeta'
start=datetime.datetime(2008, 9, 13, 6, 0, 0)


    

In [6]:

    

time0=time.time()
nc=netCDF4.Dataset(url).variables
nc.keys()
lon = nc['x'][:]
lat = nc['y'][:]
nv = nc['element'][:,:] -1
time_var = nc['time']
dtime = netCDF4.num2date(time_var[:],time_var.units)
istart = netCDF4.date2index(start,time_var,select='nearest')
var = nc[vname][istart,:]    
print 'elapsed time= %d seconds' % (time.time()-time0)


    

    




elapsed time= 17 seconds

In [7]:

    

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


    

In [11]:

    

figure(figsize=(10,10))
subplot(111,aspect=(1.0/cos(mean(lat)*pi/180.0)))
levels = linspace(0.,5.,11)
contour = tricontourf(tri, var,levels=levels,shading='faceted')
colorbar(orientation='horizontal')


    

    
Out[11]:





<matplotlib.colorbar.Colorbar instance at 0x562f998>






    

In [9]:

    

geoms = []
vals = [] # tuples of vmin,vmax

for colli,coll in enumerate(contour.collections):
    vmin,vmax = contour.levels[colli:colli+2]
    
    for p in coll.get_paths():
        p.simplify_threshold = 0.0
        polys = p.to_polygons()
        
        geoms.append( geometry.Polygon(polys[0],polys[1:] ) )
        vals.append( (vmin,vmax) )
         
wkb2shp.wkb2shp("contour-output.shp", geoms, overwrite=True,
                fields =array( vals, dtype=[('min',float64),
                                            ('max',float64)] ))


    

    




Removing the old to make way for the new