Extract data from The National Map NLCD WCS Service

In [ ]:

    

from owslib.wcs import WebCoverageService
import numpy as np
import numpy.ma as ma
endpoint='http://raster.nationalmap.gov/ArcGIS/services/LandCover/USGS_EROS_LandCover_NLCD/MapServer/WCSServer?request=GetCapabilities&service=WCS'
import matplotlib.pyplot as plt
%matplotlib inline


    

In [2]:

    

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


    

In [3]:

    

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


    

    




11 Impervious_Surface_2001_HI_11
10 Impervious_Surface_2001_AK_10
13 Impervious_Surface_2001_PR_13
12 Impervious_Surface_2001_CONUS_12
15 Canopy_2001_HI_15
14 Canopy_2001_AK_14
17 Canopy_2001_PR_17
16 Canopy_2001_CONUS_16
18 Land_Cover_1992_18
1 Land_Cover_2011_CONUS_1
3 Canopy_2011_CONUS_3
2 Impervious_Surface_2011_CONUS_2
5 Impervious_Surface_2006_CONUS_5
4 Land_Cover_2006_CONUS_4
7 Land_Cover_2001_HI_7
6 Land_Cover_2001_AK_6
9 Land_Cover_2001_CONUS_9
8 Land_Cover_2001_PR_8

In [4]:

    

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


    

    




Land_Cover_2011_CONUS_1
(-130.23284428732586, 21.74226357397646, -63.671993647433084, 52.87727272194017)
[]
[]

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()


    

    




0 129

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