Extract data from USACE WCS Service

In [1]:

    

from owslib.wcs import WebCoverageService
endpoint='http://gis.sam.usace.army.mil/server/services/JALBTCX/NCMP_BareEarth_1m/ImageServer/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


    

    




C042_NCMP_BareEarth_1m_1

In [4]:

    

wcs.contents


    

    
Out[4]:





{'1': <owslib.coverage.wcs100.ContentMetadata at 0x7f1241c18050>}

In [5]:

    

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


    

    




C942_NCMP_BareEarth_1m_1
(-124.65235145264401, 33.707727648767, -70.73480695664401, 47.941841904766996)
[]
[]

In [6]:

    

# try Plum Island Sound Region
bbox = (-70.825,42.701,-70.7526,42.762)
output = wcs.getCoverage(identifier="1",bbox=bbox,crs='EPSG:4326',format='GeoTIFF',
                         resx=0.0001, resy=0.0001)


    

In [7]:

    

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


    

In [8]:

    

from osgeo import gdal
gdal.UseExceptions()


    

In [9]:

    

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


    

In [10]:

    

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

    

import cartopy.crs as ccrs


    

In [12]:

    

print x0,x1,y1,y0


    

    




-70.825 -70.7526 42.701 42.762

In [13]:

    

plt.figure(figsize=(8,8))
ax = plt.axes(projection=ccrs.PlateCarree())
plt.imshow(elevation, cmap='jet', extent=[x0, x1, y1, y0],transform=ccrs.PlateCarree());
ax.gridlines(draw_labels=True);


    

In [13]: