In [1]:
from owslib.wcs import WebCoverageService
import numpy as np
import numpy.ma as ma
#endpoint='http://gis.sam.usace.army.mil/server/services/JALBTCX/NCMP_BareEarth_1m/ImageServer/WCSServer?request=GetCapabilities&service=WCS'
endpoint='http://gis.sam.usace.army.mil/server/services/JALBTCX/NCMP_2005_MA_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
In [4]:
wcs.contents
Out[4]:
In [5]:
lidar = wcs['1']
print lidar.title
print lidar.boundingBoxWGS84
print lidar.timelimits
print lidar.supportedFormats
In [6]:
# try Plum Island Sound Region
bbox = (-70.825,42.671,-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
from cartopy.io.img_tiles import MapQuestOpenAerial, MapQuestOSM, OSM
In [12]:
print x0,x1,y1,y0
In [13]:
elevation=ma.masked_equal(elevation,0)
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);