Note: requires OSWLib 0.7_dev or higher
In [1]:
from pylab import *
from owslib.csw import CatalogueServiceWeb
from owslib import fes
import netCDF4
import pandas as pd
In [2]:
from IPython.core.display import HTML
HTML('<iframe src=http://geoport.whoi.edu/geoportal/ width=950 height=400></iframe>')
Out[2]:
In [3]:
# connect to CSW, explore it's properties
#endpoint = 'http://www.ngdc.noaa.gov/geoportal/csw' # NGDC Geoportal
#endpoint = 'http://www.nodc.noaa.gov/geoportal/csw' # NODC Geoportal: granule level
#endpoint = 'http://data.nodc.noaa.gov/geoportal/csw' # NODC Geoportal: collection level
#endpoint = 'http://geodiscover.cgdi.ca/wes/serviceManagerCSW/csw' # NRCAN CUSTOM
#endpoint = 'http://geoport.whoi.edu/gi-cat/services/cswiso' # USGS Woods Hole GI_CAT
#endpoint = 'http://cida.usgs.gov/gdp/geonetwork/srv/en/csw' # USGS CIDA Geonetwork
endpoint = 'http://geoport.whoi.edu/geoportal/csw'
csw = CatalogueServiceWeb(endpoint,timeout=30)
csw.version
Out[3]:
In [4]:
[op.name for op in csw.operations]
Out[4]:
In [5]:
# hopefully something like this will be implemented in fes soon
def dateRange(start_date='1900-01-01',stop_date='2100-01-01',constraint='overlaps'):
if constraint == 'overlaps':
start = fes.PropertyIsLessThanOrEqualTo(propertyname='startDate', literal=stop_date)
stop = fes.PropertyIsGreaterThanOrEqualTo(propertyname='endDate', literal=start_date)
elif constraint == 'within':
start = fes.PropertyIsGreaterThanOrEqualTo(propertyname='startDate', literal=start_date)
stop = fes.PropertyIsLessThanOrEqualTo(propertyname='endDate', literal=stop_date)
return start,stop
In [6]:
# Perform the CSW query, using Kyle's cool new filters on ISO queryables
# find all datasets in a bounding box and temporal extent that have
# specific keywords and also can be accessed via OPeNDAP
box = [-70.5, 43.0, -68.0, 44.0]
bbox = fes.BBox(box)
start,stop = dateRange('1993-01-01','2001-09-01')
long_name = 'ATTENUATION'
#std_name = 'sea_water_temperature'
#keywords = fes.PropertyIsLike(propertyname='anyText', literal=std_name)
serviceType = fes.PropertyIsLike(propertyname='apiso:ServiceType', literal='*opendap*')
# apply all the filters using the "and" syntax: [[filter1,filter2]]
csw.getrecords2(constraints=[[start,stop,serviceType,bbox]],maxrecords=150,esn='full')
#csw.get_filtered_records([[keywords,start,stop,serviceType,bbox]],maxrecords=5)
csw.records.keys()
Out[6]:
In [7]:
for rec,item in csw.records.iteritems():
print item.title
In [8]:
# get specific ServiceType URL from records
def service_urls(records,service_string='urn:x-esri:specification:ServiceType:odp:url'):
urls=[]
for key,rec in records.iteritems():
#create a generator object, and iterate through it until the match is found
#if not found, gets the default value (here "none")
url = next((d['url'] for d in rec.references if d['scheme'] == service_string), None)
if url is not None:
urls.append(url)
return urls
In [9]:
dap_urls = service_urls(csw.records,service_string='urn:x-esri:specification:ServiceType:odp:url')
print ".html\n".join(dap_urls)
In [10]:
def standard_names(nc):
'''
get dictionary of variables with standard_names
'''
d={}
for k,v in nc.iteritems():
try:
standard_name=v.getncattr('standard_name')
try:
d[standard_name]=[d[standard_name],[k]]
except:
d[standard_name]=[k]
except:
pass
return d
In [11]:
def long_names(nc):
'''
get dictionary of variables with long_names
'''
d={}
for k,v in nc.iteritems():
try:
long_name=v.getncattr('long_name').strip()
try:
d[long_name]=[d[long_name],[k]]
except:
d[long_name]=[k]
except:
pass
return d
In [12]:
def names(nc):
'''
get dictionary of variables with names
'''
d={}
for k,v in nc.iteritems():
try:
name=v.getncattr('name').strip()
try:
d[name]=[d[name],[k]]
except:
d[name]=[k]
except:
pass
return d
In [13]:
dap_urls = [url for url in dap_urls if '-a1h_1d' in url]
print ".html\n".join(dap_urls)
In [14]:
long_name='ATTENUATION'
In [15]:
for url in dap_urls:
ncd = netCDF4.Dataset(url)
nc = ncd.variables
# make a dictionary containing all data from variables that matched the standard_name
# find list of variables for each standard_name
d = long_names(nc)
# find all the variables matching standard_name=std_name
if long_name in d:
# read all the data into a dictionary
lat = nc['lat'][:]
lon = nc['lon'][:]
depth = nc['depth'][:]
try:
water_depth = str(ncd.WATER_DEPTH)
except:
try:
water_depth = str(ncd.water_depth)
except:
water_depth = 'unknown'
time_var = nc['time']
dtime = netCDF4.num2date(time_var[:],time_var.units)
data_dict={}
for v in d[long_name]:
data_dict[v]=nc[v][:].flatten()
print url
# Create Pandas data frame, with time index
ts = pd.DataFrame.from_dict(data_dict)
ts.index=dtime
ts.plot(figsize=(12,4));
title('%s, %s, depth=%4.2f, water_depth=%s' % (url[-30:],long_name, depth, water_depth))
In [16]:
ncd = netCDF4.Dataset(url)
In [17]:
ncd.WATER_DEPTH
Out[17]:
In [18]:
name='SDP'
for url in dap_urls:
ncd = netCDF4.Dataset(url)
nc = ncd.variables
# make a dictionary containing all data from variables that matched the standard_name
# find list of variables for each standard_name
d = names(nc)
# find all the variables matching standard_name=std_name
if name in d:
# read all the data into a dictionary
lat = nc['lat'][:]
lon = nc['lon'][:]
depth = nc['depth'][:]
try:
water_depth = str(ncd.WATER_DEPTH)
except:
try:
water_depth = str(ncd.water_depth)
except:
water_depth = 'unknown'
time_var = nc['time']
dtime = netCDF4.num2date(time_var[:],time_var.units)
data_dict={}
for v in d[name]:
data_dict[v]=nc[v][:].flatten()
print url
# Create Pandas data frame, with time index
ts = pd.DataFrame.from_dict(data_dict)
ts.index=dtime
ts.plot(figsize=(12,4));
title('%s, %s, depth=%4.2f, water_depth=%s' % (url[-30:],name, depth, water_depth))
In [19]:
import iris
import iris.plot as iplt
import iris.quickplot as qplt
bathy_url='http://geoport.whoi.edu/thredds/dodsC/bathy/gom03_v1_0'
c = iris.load(bathy_url)
In [20]:
slice=c[0].extract(iris.Constraint(
longitude=lambda cell: box[0] < cell < box[2],
latitude=lambda cell: box[1] < cell < box[3]))
In [21]:
url='http://geoport-dev.whoi.edu/thredds/dodsC/ECOHAB_II/6226att-a1h_1d.nc'
ncd=netCDF4.Dataset(url)
nc = ncd.variables
In [22]:
lon=nc['lon'][:]
lat=nc['lat'][:]
In [23]:
figure(figsize=(12,8))
latb = slice.coord(axis='Y').points
lonb = slice.coord(axis='X').points
subplot(111,aspect=(1.0/cos(mean(latb)*pi/180.0)))
pcolormesh(lonb,latb,ma.masked_invalid(slice.data),vmin=-200,vmax=0)
colorbar()
contour(lonb,latb,ma.masked_invalid(slice.data),[-100, -50, 0],colors=['k','k','k'])
plot(lon,lat,'mo',markersize=12)
Out[23]:
In [24]:
long_name='STAND. DEV. (PRESS)'
In [25]:
lat
Out[25]:
In [26]:
url='http://geoport.whoi.edu/thredds/dodsC/ECOHAB_II/5211psd-a1h.nc'
ncd = netCDF4.Dataset(url)
nc = ncd.variables
In [ ]:
nc.keys()
In [ ]:
print nc['SDP_850']
In [ ]:
dap_urls
In [ ]:
ncd
In [ ]:
('WATER_DEPTH' in ncd.ncattrs()) or ('water_depth' in ncd.ncattrs())
In [ ]: