UGRID Remote Subsetting Demo with Time

Much of the power of OPeNDAP is the ability to select and subset just the data you need from a remote dataset. This works great for structured grid data because one can subset on index ranges, but for unstructured grid (e.g. triangular grid) this is impossible. It's therefore necessary to add an additional constraint syntax to the URL and do the subsetting on the server. This is one of the goals of the OPeNDAP-Unidata Linked Servers OPULS project, a collaboration between Unidata, OPeNDAP, Inc. and the University of Washington eScience Institute.

In this demo, we show how the subset expression works, extracting the Galveston Bay region from a model simulation of the entire Gulf of Mexico. The model dataset is UGRID-compliant FVCOM model output from the IOOS Modeling Testbed project. A sample dataset was moved from the testbed THREDDS Data Server to a developmental Hyrax server in the Cloud which has the subsetting enabled. The subsetting is accomplished on the server backend using GridFields.



In [1]:

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



In [10]:

    
# 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'
url='http://ec2-54-242-224-73.compute-1.amazonaws.com:8080/opendap/ebs/Ike/2D_varied_manning_windstress/test_dir-norename.ncml'



In [11]:

    
# 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



In [12]:

    
# construct the subset expression
# here extract every 3rd time step between 3 and 10
expr='_expr_{}{ugr5(0,ua[3:3:10][*],va[3:3:10][*],zeta[3:3:10][*],\"%.6f<lat&lat<%.6f&%.6f<lon&lon<%.6f\")}{}' % (bbox[2],bbox[3],bbox[0],bbox[1])



In [16]:

    
# construct the subset expression
# here extract every 3rd time step between 3 and 10
expr='_expr_{}{ugr5(0,zeta[3:3:10][*],\"%.6f<lat&lat<%.6f&%.6f<lon&lon<%.6f\")}{}' % (bbox[2],bbox[3],bbox[0],bbox[1])



In [17]:

    
# construct the subset url 
url_subset = url + expr



In [18]:

    
print url_subset









    



http://ec2-54-242-224-73.compute-1.amazonaws.com:8080/opendap/ebs/Ike/2D_varied_manning_windstress/test_dir-norename.ncml_expr_{}{ugr5(0,zeta[3:3:10][*],"29.300000<lat&lat<29.800000&-95.000000<lon&lon<-94.400000")}{}



In [19]:

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









    



---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-19-d8ffa62d15c9> in <module>()
      1 time0 = time.time()
----> 2 ncs = netCDF4.Dataset(url_subset)
      3 lons=ncs.variables['ugr_result.lon'][:]
      4 lats=ncs.variables['ugr_result.lat'][:]
      5 lonc=ncs.variables['ugr_result.lonc'][:]

/opt/anaconda/envs/np17py27-1.5/lib/python2.7/site-packages/netCDF4.so in netCDF4.Dataset.__init__ (netCDF4.c:19148)()

RuntimeError: NetCDF: I/O failure



In [9]:

    
print shape(ua)
print shape(zeta)









    



---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-9-dc86fcf08ccd> in <module>()
----> 1 print shape(ua)
      2 print shape(zeta)

NameError: name 'ua' is not defined



In [ ]:

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



In [ ]:

    
# plot subset water level and vectors

figure(figsize=(10,8),frameon=True)
# plot one of the time steps 
time_step=-1
tricontourf(lons,lats,zeta[time_step,:].flatten(),triangles=nvs,levels=arange(-1,1,.1))
gca().set_aspect(1./cos(mean(lats)*pi/180))
colorbar()
Q = quiver(lonc[idv],latc[idv],ua[time_step,idv].flatten(),va[time_step,idv].flatten(),scale=5)
qk = quiverkey(Q,0.82,0.92,0.20,'0.2 m/s',labelpos='W')
axis(bbox)
title('Galveston Bay Water Level (m) and Depth-Averaged Velocity');



In [ ]: