Extract data from USGS ScienceBase

In [1]:

    

%matplotlib inline


    

In [11]:

    

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 [12]:

    

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


    

In [13]:

    

wcs.contents


    

    




---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-13-fb40f44705ea> in <module>()
----> 1 wcs.contents

AttributeError: 'NoneType' object has no attribute 'contents'

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 [ ]: