Extract data from USACE WCS Service

In [1]:

    

import matplotlib.pyplot as plt
from owslib.wcs import WebCoverageService
%matplotlib inline


    

In [2]:

    

endpoint='http://cloud.insideidaho.org/ArcGIS/services/climatologyMeteorologyAtmosphere/climatologyMeteorologyAtmosphere/MapServer/WCSServer?request=GetCapabilities&service=WCS'


    

In [3]:

    

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


    

In [4]:

    

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


    

    




Agroclimate Zones_1
Koppen Climate Classification_2

In [5]:

    

wcs.contents


    

    
Out[5]:





{'1': <owslib.coverage.wcs100.ContentMetadata at 0x7fc3a408e5d0>,
 '2': <owslib.coverage.wcs100.ContentMetadata at 0x7fc3a408e610>}

In [6]:

    

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


    

    




Agroclimate Zones_1
(-117.54391664332157, 41.93095065270037, -110.64024328374633, 49.054486221179275)
[]
[]

In [25]:

    




    

In [7]:

    

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


    

In [8]:

    

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


    

In [9]:

    

from osgeo import gdal
gdal.UseExceptions()


    

In [10]:

    

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


    

In [11]:

    

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

    

import cartopy.crs as ccrs


    

In [13]:

    

print x0,x1,y1,y0


    

    




-117.543916643 -110.640243284 41.9309506527 49.0544862212

In [14]:

    

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