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.
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
Elapsed time to read full grid: 17.56 seconds
# 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)')
<matplotlib.text.Text at 0x2889b10>
# 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
http://ec2-54-242-224-73.compute-1.amazonaws.com:8080/opendap/ebs/fvcom_1step.nc_expr_{}{ugr3(0,ua,va,h,"29.300000<lat&lat<29.800000&-95.000000<lon&lon<-94.400000")}{}
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
Elapsed time to read subset grid: 8.02 seconds
# 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()
<matplotlib.colorbar.Colorbar instance at 0x349fa28>
# 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()
<matplotlib.colorbar.Colorbar instance at 0x3c38b48>
# 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
http://geoport-dev.whoi.edu/thredds/wms/bathy/crm_vol5.nc?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&BBOX=-95.000000,29.300000,-94.400000,29.800000
Image(url=wms_get_map)
ncs.variables.keys()
[u'ugr_result.lon', u'ugr_result.lat', u'ugr_result.lonc', u'ugr_result.latc', u'ugr_result.nv', u'ugr_result.fvcom_mesh', u'ugr_result.ua', u'ugr_result.va', u'ugr_result.h']
print nc.variables['ua']
<type 'netCDF4.Variable'> float32 ua(u'time', u'nele') long_name: Vertically Averaged x-velocity units: meters s-1 type: data missing_value: -999.0 field: ua, scalar coverage_content_type: modelResult standard_name: barotropic_eastward_sea_water_velocity coordinates: time latc lonc mesh: fvcom_mesh location: face _Netcdf4Dimid: 0 unlimited dimensions = (u'time',) current size = (1, 826866)
print ncs.variables['ugr_result.ua']
<type 'netCDF4.Variable'> float64 ugr_result.ua(u'time', u'nele') long_name: Vertically Averaged x-velocity units: meters s-1 type: data missing_value: -999.0 field: ua, scalar coverage_content_type: modelResult standard_name: barotropic_eastward_sea_water_velocity coordinates: time latc lonc mesh: fvcom_mesh location: face _Netcdf4Dimid: 0 unlimited dimensions = (u'time',) current size = (1, 11352)
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)')
--------------------------------------------------------------------------- IndexError Traceback (most recent call last) <ipython-input-21-23bd41272003> in <module>() 4 gca().set_aspect(1./cos(mean(lats)*pi/180)) 5 colorbar() ----> 6 Q = quiver(lonc[idv],latc[idv],ua[idv],va[idv],scale=5) 7 qk = quiverkey(Q,0.82,0.92,0.20,'0.2 m/s',labelpos='W') 8 axis(bbox) /home/local/python27_epd/lib/python2.7/site-packages/numpy/ma/core.pyc in __getitem__(self, indx) 2939 # raise IndexError(msg) 2940 _data = ndarray.view(self, ndarray) -> 2941 dout = ndarray.__getitem__(_data, indx) 2942 # We could directly use ndarray.__getitem__ on self... 2943 # But then we would have to modify __array_finalize__ to prevent the IndexError: index 10297 is out of bounds for axis 0 with size 1