Extract data from USGS ScienceBase



In [1]:

    
%matplotlib inline



In [2]:

    
from owslib.wcs import WebCoverageService
import numpy as np
import numpy.ma as ma
endpoint='https://www.sciencebase.gov/catalogMaps/mapping/ows/5638cf1fe4b0d6133fe73040?service=wcs&request=getcapabilities&version=1.1.1'



In [3]:

    
wcs = WebCoverageService(endpoint,version='1.1.1',timeout=60)









    



/home/usgs/miniconda/envs/ioos/lib/python2.7/site-packages/owslib/coverage/wcs110.py:84: FutureWarning: The behavior of this method will change in future versions. Use specific 'len(elem)' or 'elem is not None' test instead.
  elem=self._capabilities.find(self.ns.OWS('ServiceProvider')) or self._capabilities.find(self.ns.OWS('ServiceProvider'))



In [4]:

    
wcs.contents









    Out[4]:





{'bh_30mbathy': <owslib.coverage.wcs110.ContentMetadata at 0x7f8a8d615310>}



In [5]:

    
for k,v in wcs.contents.iteritems():
    print v.title









    



bh_30mbathy



In [6]:

    
g = wcs['bh_30mbathy']
print g.title
print g.boundingBoxWGS84
print g.supportedFormats









    



bh_30mbathy
(-71.05592748611865, 42.25694835365398, -70.81459673354432, 42.43833963976761)
[]



In [7]:

    
# try Boston Harbor
bbox = (-71.05592748611777, 42.256890708126605, -70.81446033774644, 42.43833963977496)
output = wcs.getCoverage(identifier="bh_30mbathy",bbox=bbox,crs='EPSG:4326',format='GeoTIFF',
                         resx=0.0003, resy=0.0003)









    



---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-7-a1e534669cc5> in <module>()
      2 bbox = (-71.05592748611777, 42.256890708126605, -70.81446033774644, 42.43833963977496)
      3 output = wcs.getCoverage(identifier="bh_30mbathy",bbox=bbox,crs='EPSG:4326',format='GeoTIFF',
----> 4                          resx=0.0003, resy=0.0003)

/home/usgs/miniconda/envs/ioos/lib/python2.7/site-packages/owslib/coverage/wcs110.pyc in getCoverage(self, identifier, bbox, time, format, store, rangesubset, gridbaseCRS, gridtype, gridCS, gridorigin, gridoffsets, method, **kwargs)
    158             method=self.ns.WCS_OWS('Get')
    159         try:
--> 160             base_url = next((m.get('url') for m in self.getOperationByName('GetCoverage').methods if m.get('type').lower() == method.lower()))
    161         except StopIteration:
    162             base_url = self.url

/home/usgs/miniconda/envs/ioos/lib/python2.7/site-packages/owslib/coverage/wcs110.pyc in getOperationByName(self, name)
    208             if item.name == name:
    209                 return item
--> 210         raise KeyError("No operation named %s" % name)
    211 
    212 class Operation(object):

KeyError: 'No operation named GetCoverage'



In [ ]:

    
wcs.



In [ ]:

    
f=open('test.tif','wb')
f.write(output.read())
f.close()



In [ ]:

    
from osgeo import gdal
gdal.UseExceptions()



In [ ]:

    
ds = gdal.Open('test.tif')



In [ ]:

    
band = ds.GetRasterBand(1)
elevation = band.ReadAsArray()
nrows, ncols = elevation.shape

# I'm making the assumption that the image isn't rotated/skewed/etc. 
# This is not the correct method in general, but let's ignore that for now
# If dxdy or dydx aren't 0, then this will be incorrect
x0, dx, dxdy, y0, dydx, dy = ds.GetGeoTransform()

if dxdy == 0.0:
    x1 = x0 + dx * ncols
    y1 = y0 + dy * nrows



In [ ]:

    
import cartopy.crs as ccrs
from cartopy.io.img_tiles import MapQuestOpenAerial, MapQuestOSM, OSM



In [ ]:

    
print x0,x1,y1,y0



In [ ]:

    
elevation=ma.masked_less_equal(elevation,-1.e5)



In [ ]:

    
print elevation.min(), elevation.max()



In [ ]:

    
plt.figure(figsize=(8,10))
ax = plt.axes(projection=ccrs.PlateCarree())
tiler = MapQuestOpenAerial()
ax.add_image(tiler, 14)
plt.imshow(elevation, cmap='jet', extent=[x0, x1, y1, y0],
           transform=ccrs.PlateCarree(),alpha=0.6,zorder=2);
ax.gridlines(draw_labels=True,zorder=3);



In [ ]: