In [1]:

    

import matplotlib.tri as tri
import datetime as dt
import matplotlib.pyplot as plt
import numpy as np


    

In [2]:

    

import cartopy.crs as ccrs
%matplotlib inline


    

    




---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-2-272d320a3b86> in <module>()
----> 1 import cartopy.crs as ccrs
      2 get_ipython().magic(u'matplotlib inline')

ImportError: No module named cartopy.crs

In [3]:

    

import iris
import pyugrid


    

In [4]:

    

iris.FUTURE.netcdf_promote = True


    

In [5]:

    

#ADCIRC
#url =  'http://comt.sura.org/thredds/dodsC/data/comt_1_archive/inundation_tropical/UND_ADCIRC/Hurricane_Ike_3D_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'


    

In [6]:

    

ug = pyugrid.UGrid.from_ncfile(url)

print "There are %i nodes"%ug.nodes.shape[0]
print "There are %i faces"%ug.faces.shape[0]


    

    




There are 592761 nodes
There are 1166747 faces

In [7]:

    

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


    

    




/home/usgs/miniconda/envs/ioos/lib/python2.7/site-packages/iris/fileformats/_pyke_rules/compiled_krb/fc_rules_cf_fc.py:1300: UserWarning: Ignoring netCDF variable u'ele' invalid units u'non-dimensional'
  warnings.warn(msg.encode('ascii', errors='backslashreplace'))
/home/usgs/miniconda/envs/ioos/lib/python2.7/site-packages/iris/fileformats/_pyke_rules/compiled_krb/fc_rules_cf_fc.py:1300: UserWarning: Ignoring netCDF variable u'Cs' invalid units u'non-dimensional'
  warnings.warn(msg.encode('ascii', errors='backslashreplace'))
/home/usgs/miniconda/envs/ioos/lib/python2.7/site-packages/iris/fileformats/_pyke_rules/compiled_krb/fc_rules_cf_fc.py:1300: UserWarning: Ignoring netCDF variable u'sigma' invalid units u'non-dimensional'
  warnings.warn(msg.encode('ascii', errors='backslashreplace'))

In [8]:

    

print cube


    

    




sea_surface_height_above_geoid / (m) (time: 120; -- : 592761)
     Dimension coordinates:
          time                            x         -
     Auxiliary coordinates:
          latitude                        -         x
          longitude                       -         x
     Attributes:
          Conventions: CF-1.0
          NCO: 4.0.8
          cdm_data_type: ugrid
          comment: David Forrest drf@vims.edu 2012-12-13
          history: ncml aggregation of .../whole_domain/ datafiles
          id: inundation_tropical.VIMS_SELFE.Hurricane_Ike_2D_final_run_with_waves
          institution: Virginia Institute of Marine Science -- http://web.vims.edu/physical/
          location: node
          mesh: selfe_mesh
          ncmlFile: /data/comt_1_archive/inundation_tropical/VIMS_SELFE/Hurricane_Ike_2D_f...
          nco_openmp_thread_number: 1
          references: testbed.sura.org:/data/ftp/upload/Inundation/vims/selfe_tropical/runs/Ike/2D_varied_manning_windstress/new_10min/Inundation_Model_Metadata_Template_v3_for_2D_windstress_sv_Ike_10min.docx...
          source: SELFE run  (version v.vvv  compiled yyyy-mm-dd) for Hurricane Ike on the...
          summary: A 2D hindcast of Hurricane Ike (2008) using WWMII+SELFE on the ULLR mesh...
          title: Inundation Tropical : VIMS : SELFE : Hurricane Ike 2D final run with w...

In [9]:

    

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


    

In [10]:

    

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


    

In [11]:

    

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


    

In [12]:

    

ind = -1 # last time index
zcube = cube[ind]


    

In [13]:

    

plt.figure(figsize=(12,12))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([-90, -60, 5, 50])
ax.coastlines()
levs = np.arange(-1,5,.2)
plt.tricontourf(triang, zcube.data, levels=levs)
plt.colorbar()
plt.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
plt.title('%s: Elevation (m): %s' % (zcube.attributes['title'],tstr));


    

In [13]: