In [1]:

    

import matplotlib.tri as tri
import datetime as dt


    

In [2]:

    

import cartopy.crs as ccrs
import iris
iris.FUTURE.netcdf_promote = True
import pyugrid


    

In [3]:

    

#ADCIRC
url =  'http://comt.sura.org/thredds/dodsC/data/comt_1_archive/inundation_tropical/UND_ADCIRC/Hurricane_Ike_2D_final_run_with_waves'

#FVCOM
#url = 'http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_GOM3_FORECAST.nc'

#SELFE
#url = 'http://comt.sura.org/thredds/dodsC/data/comt_1_archive/inundation_tropical/VIMS_SELFE/Hurricane_Ike_2D_final_run_with_waves'

#WW3
url='http://comt.sura.org/thredds/dodsC/data/comt_2/pr_inundation_tropical/EMC_ADCIRC-WW3/Dec2013Storm_2D_preliminary_run_1_waves_only'


    

In [4]:

    

cube = iris.load_cube(url,'sea_surface_height_above_geoid')


    

    




/home/testbed/anaconda/lib/python2.7/site-packages/iris/fileformats/_pyke_rules/compiled_krb/fc_rules_cf_fc.py:1195: UserWarning: Ignoring netCDF variable 'Fp3D' invalid units 'log10(m4s+0.01)'
  warnings.warn(msg.format(msg_name, msg_units))

In [5]:

    

print cube


    

    




sea_surface_height_above_geoid / (m) (time: 8; -- : 254565)
     Dimension coordinates:
          time                            x       -
     Auxiliary coordinates:
          latitude                        -       x
          longitude                       -       x
     Attributes:
          Conventions: CF-1.6, UGRID-0.9
          WAVEWATCH_III_switches: F90 NOGRB UG NC3 LRB4 MPI DIST NOPA PR3 ST4 STAB0 FLX0 NL1 BT4 DB1 MLIM...
          WAVEWATCH_III_version_number: 4.08
          _DODS_Unlimited_Dimension: time
          altitude_resolution: n/a
          area: EC2001
          cdm_data_type: 
          easternmost_longitude: -60.0400009
          globwave_name: sea_surface_height_above_sea_level
          id: pr_inundation_tropical.EMC_ADCIRC-WW3.Dec2013Storm_2D_preliminary_run_...
          latitude_resolution: n/a
          location: node
          longitude_resolution: n/a
          maximum_altitude: 9000 m
          mesh: ww3_mesh
          minimum_altitude: -12000 m
          northernmost_latitude: 45.8388062
          product_name: ww3.20131006_cur.nc
          southernmost_latitude: 7.9111028
          start_date: 2013-10-06T00:00:00Z
          stop_date: 2013-10-06T21:00:00Z
          summary: A 2D test using ADCIRC+WW3 
          title: Pr Inundation Tropical : EMC : ADCIRC-WW3 : Dec2013Storm 2D preliminary...
          westernmost_longitude: -97.8568878

In [6]:

    

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


    

In [7]:

    

ug = pyugrid.UGrid.from_ncfile(url)

# What's in there?
#print "There are %i nodes"%ug.nodes.shape[0]
#print "There are %i edges"%ug.edges.shape[0]
#print "There are %i faces"%ug.faces.shape[0]


    

In [8]:

    

cube.mesh = ug
cube.mesh_dimension = 1  # (0:time,1:node)


    

In [9]:

    

lon = cube.mesh.nodes[:,0]
lat = cube.mesh.nodes[:,1]
nv = cube.mesh.faces


    

In [10]:

    

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


    

In [11]:

    

# skip trying to find the closest time index to requested time, because it's messy
ind = -1 # just take the last time index for now
zcube = cube[ind]


    

In [12]:

    

figure(figsize=(12,12))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([lon.min(), lon.max(), lat.min(), lat.max()])
ax.coastlines()
levs=arange(-1,5,.2)
tricontourf(triang, zcube.data, levels=levs)
colorbar()
tricontour(triang, zcube.data, colors='k',levels=levs)
tvar = cube.coord('time')
tstr = tvar.units.num2date(tvar.points[ind])
gl = ax.gridlines(draw_labels=True)
gl.xlabels_top = False
gl.ylabels_right = False
title('%s: Elevation (m): %s' % (zcube.attributes['title'],tstr));


    

    




/home/testbed/anaconda/lib/python2.7/site-packages/cartopy/io/__init__.py:261: DownloadWarning: Downloading: http://www.nacis.org/naturalearth/110m/physical/ne_110m_coastline.zip
  warnings.warn('Downloading: {}'.format(url), DownloadWarning)






    

In [12]: