from pylab import *
import time
# The netCDF4-Python library accesses OPeNDAP using the Unidata NetCDF-C library:
import netCDF4

# the full dataset OPeNDAP URL:
#url='http://ec2-54-245-151-123.us-west-2.compute.amazonaws.com:8080/opendap/ugrids/fvcom_1step.nc'
url='http://ec2-54-242-224-73.compute-1.amazonaws.com:8080/opendap/ebs/fvcom_1step.nc'

# the region we want to subset, in Matlab style bbox:
bbox=[-95.0, -94.4, 29.3, 29.8]  # [lonmin lonmax latmin latmax] Galveston Bay

# open and read a field from the full dataset:
time0 = time.time()
nc = netCDF4.Dataset(url)
nc.variables
lon = nc.variables['lon'][:]
lat = nc.variables['lat'][:]
# read connectivity array, and convert to pythonic 0-based indexing
nv = nc.variables['nv'][:] - nc.variables['nv'].start_index
h = nc.variables['h'][:]
etime = time.time()- time0
print 'Elapsed time to read full grid: %.2f seconds' % etime

# plot full bathy:
figure(figsize=(10,6),frameon=True)
tricontourf(lon,lat,-h,triangles=nv,
   levels=[-5000,-4000,-3000,-2000,-1000,-500,-100,-50,-20,-10,0])
colorbar()
gca().set_aspect(1./cos(30.0*pi/180))
title('FVCOM Bathy (m)')

# construct the subset expression
expr='_expr_{}{ugr3(0,ua,va,h,\"%.6f<lat&lat<%.6f&%.6f<lon&lon<%.6f\")}{}' % (bbox[2],bbox[3],bbox[0],bbox[1])

# construct the subset url 
url_subset = url + expr

print url_subset

time0 = time.time()
ncs = netCDF4.Dataset(url_subset)
lons=ncs.variables['ugr_result.lon'][:]
lats=ncs.variables['ugr_result.lat'][:]
# read connectivity array, and convert to pythonic 0-based indexing
nvs = ncs.variables['ugr_result.nv'][:] - ncs.variables['ugr_result.nv'].start_index
hs=ncs.variables['ugr_result.h'][:]
etime = time.time()- time0
print 'Elapsed time to read subset grid: %.2f seconds' % etime

# plot subset bathy:
figure(figsize=(8,6),frameon=True)
tricontourf(lons,lats,-hs,triangles=nvs,levels=range(-20,0))
gca().set_aspect(1./cos(30.0*pi/180))
axis(bbox)
title('Galveston Bay subset bathy (m)')
colorbar()

# plot subset bathy with mesh drawn
figure(figsize=(8,6),frameon=True)
tricontourf(lons,lats,-hs,triangles=nvs,levels=range(-20,0))
triplot(lons,lats,triangles=nvs)
gca().set_aspect(1./cos(30.0*pi/180))
axis(bbox)
title('Galveston Bay subset bathy with mesh indicated(m)')
colorbar()

# Plot WMS bathy image from NOAA Coastal Relief Model
from IPython.core.display import Image
wms = 'http://geoport-dev.whoi.edu/thredds/wms/bathy/crm_vol5.nc?'
wms_query = 'LAYERS=topo&ELEVATION=0&TIME=2012-12-07T08%3A42%3A13Z&TRANSPARENT=true&STYLES=boxfill%2Frainbow&CRS=EPSG%3A4326&COLORSCALERANGE=-20%2C0&NUMCOLORBANDS=32&LOGSCALE=false&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&EXCEPTIONS=application%2Fvnd.ogc.se_inimage&FORMAT=image%2Fpng&SRS=EPSG%3A4326&WIDTH=360&HEIGHT=360'
wms_bbox = '&BBOX=%.6f,%.6f,%.6f,%.6f' % (bbox[0],bbox[2],bbox[1],bbox[3])
wms_get_map = wms + wms_query + wms_bbox
print wms_get_map

Image(url=wms_get_map)

ncs.variables.keys()

print nc.variables['ua']

print ncs.variables['ugr_result.ua']

ua=ncs.variables['ugr_result.ua'][:]

va=ncs.variables['ugr_result.va'][:]

lonc=mean(lons[nvs],axis=0)
latc=mean(lats[nvs],axis=0)

# subsample velocity vectors 
ind=range(len(lonc))
subsample=3
np.random.shuffle(ind)
Nvec = int(len(ind) / subsample)
idv = ind[:Nvec]

# plot subset bathy:
figure(figsize=(10,8),frameon=True)
tricontourf(lons,lats,-hs,triangles=nvs,levels=range(-20,1))
gca().set_aspect(1./cos(mean(lats)*pi/180))
colorbar()
Q = quiver(lonc[idv],latc[idv],ua[idv],va[idv],scale=5)
qk = quiverkey(Q,0.82,0.92,0.20,'0.2 m/s',labelpos='W')
axis(bbox)
title('Galveston Bay subset bathy (m)')