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 [21]:
# 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/hyrax/ebs/Ike/2D_varied_manning_windstress/test_dir-norename.ncml'
In [3]:
# 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 [4]:
# 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['depth'][:]
etime = time.time()- time0
print 'Elapsed time to read full grid: %.2f seconds' % etime
In [5]:
# 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)')
Out[5]:
In [6]:
# construct the subset expression
expr='_expr_{}{ugr5(0,ua,\"%.6f<lat&lat<%.6f&%.6f<lon&lon<%.6f\")}{}' % (bbox[2],bbox[3],bbox[0],bbox[1])
In [7]:
# construct the subset expression
expr='?ugr5(0,ua,\"%.6f<lat&lat<%.6f&%.6f<lon&lon<%.6f\")' % (bbox[2],bbox[3],bbox[0],bbox[1])
In [23]:
# construct the subset url
url_subset = url + '.dods' + expr
In [24]:
print url_subset
In [10]:
url='http://ec2-54-242-224-73.compute-1.amazonaws.com:8080/opendap/ebs/Ike/2D_varied_manning_windstress/test_dir-norename.ncml?ugr5(0,zeta[3:3:30][*],\"29.3<lat&lat<29.8&-95.0<lon&lon<-94.4\")'
print url
In [25]:
nc = netCDF4.Dataset(url_subset)
In [ ]:
url='http://ec2-54-242-224-73.compute-1.amazonaws.com:8080/opendap/ebs/Ike/2D_varied_manning_windstress/test_dir-norename.ncml'
nc = netCDF4.Dataset(url)
In [ ]:
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.depth'][:]
etime = time.time()- time0
print 'Elapsed time to read subset grid: %.2f seconds' % etime
In [ ]:
# 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()
In [ ]:
# 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()
In [ ]:
# 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
In [ ]:
Image(url=wms_get_map)
In [ ]:
ncs.variables.keys()
In [ ]:
print nc.variables['ua']
In [ ]:
print ncs.variables['ugr_result.ua']
In [ ]:
ua=ncs.variables['ugr_result.ua'][:]
In [ ]:
va=ncs.variables['ugr_result.va'][:]
In [ ]:
lonc=mean(lons[nvs],axis=0)
latc=mean(lats[nvs],axis=0)
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 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)')
In [ ]: