Extract data from USGS CMG WCS Service

In [1]:

    

from owslib.wcs import WebCoverageService
import numpy as np
import numpy.ma as ma
endpoint='http://coastalmap.marine.usgs.gov/cmgp/services/EastCoast/Mass_Seafloor/MapServer/WCSServer?request=GetCapabilities&service=WCS'


    

In [2]:

    

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


    

In [3]:

    

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


    

    




Northern Cape Cod Bay Swath backscatter 1m_24
Red Brook Harbor backscatter 1m_25
Buzzards Bay backscatter 1m_26
western Elizabeth Islands Buzzards Bay backscatter 0.5m_27
Boston Harbor backscatter 1m_20
Duxbury-Hull backscatter 1m_21
Duxbury-Hull NOAA H10993 backscatter 1m_22
Northern Cape Cod Bay Klein backscatter 1m_23
western Elizabeth Islands Nashawena Island backscatter 0.5m_28
western Elizabeth Islands Penikese Island backscatter 0.5m_29
Cape Ann - Salisbury bathymetric grid - 5m_1
Boston Harbor Bathymetric grid - 30m_3
Nahant Bathymetric grid - 5m_2
Northern Cape Cod Bay bathymetric grid 5m_5
Duxbury-Hull bathymetric grid 5m_4
Buzzards Bay Bathymetric grid - 5m_7
Red Brook Harbor bathymetric grid 1m_6
Vineyard Sound bathymetric grid 5m_9
western Elizabeth Islands bathymetric grid 2.5m_8
Nahant-Gloucester shaded relief_11
Cape Ann - Salisbury shaded relief_10
Duxbury-Hull shaded relief_13
Boston Harbor shaded relief_12
Buzzards Bay shaded relief 5m_15
Northern Cape Cod Bay shaded relief_14
Vineyard Sound shaded relief_17
western Elizabeth Islands bathymetric shaded relief_16
Nahant-Gloucester backscatter 1m_19
Cape Ann backscatter mosaic_18
Vineyard Sound backscatter 1m_31
western Elizabeth Islands Vineyard Sound backscatter 0.5m_30
Nahant-Gloucester bottom type_32

In [4]:

    

lidar = wcs['3']
print lidar.title
print lidar.boundingBoxWGS84
print lidar.timelimits
print lidar.supportedFormats


    

    




Boston Harbor Bathymetric grid - 30m_3
(-71.05592748611777, 42.256890708126605, -70.81446033774644, 42.43833963977496)
[]
[]

In [5]:

    

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


    

In [6]:

    

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


    

In [7]:

    

from osgeo import gdal
gdal.UseExceptions()


    

In [8]:

    

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


    

In [9]:

    

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

    

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


    

In [11]:

    

print x0,x1,y1,y0


    

    




-71.0559274861 -70.8144603377 42.2568907081 42.4383396398

In [12]:

    

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


    

In [13]:

    

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


    

    




-31.7 2.87094

In [14]:

    

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