In [1]:

    

%matplotlib inline


    

In [2]:

    

import matplotlib.pyplot as plt

import cartopy.crs as ccrs
from cartopy.io import shapereader
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

def make_map(projection=ccrs.PlateCarree()):
    fig, ax = plt.subplots(figsize=(9, 13),
                           subplot_kw=dict(projection=projection))
    gl = ax.gridlines(draw_labels=True)
    gl.xlabels_top = gl.ylabels_right = False
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    return fig, ax


    

In [3]:

    

from cartopy.feature import NaturalEarthFeature

extent = [-39, -38.25, -13.25, -12.5]

coast = NaturalEarthFeature(category='physical', scale='10m',
                            facecolor='none', name='coastline')

fig, ax = make_map(projection=ccrs.PlateCarree())

ax.set_extent(extent)

feature = ax.add_feature(coast, edgecolor='gray')


    

    




/home/usgs/miniconda/envs/ioos/lib/python2.7/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if self._edgecolors == str('face'):






    

In [4]:

    

from cartopy.io import img_tiles

tiles = img_tiles.OSM()

ax = plt.subplot(111, projection=ccrs.Mercator())
ax.set_extent(extent)
ax.add_image(tiles, 9)
ax.coastlines("10m")


    

    
Out[4]:





<cartopy.mpl.feature_artist.FeatureArtist at 0x7f2a34355fd0>






    

In [5]:

    

import netCDF4


    

In [ ]:

    




    

In [6]:

    

from owslib.wcs import WebCoverageService
import numpy as np
import numpy.ma as ma
endpoint='http://geoport.whoi.edu/thredds/wcs/bathy/gom03_v1_0?service=WCS&version=1.0.0&request=GetCapabilities'


    

In [7]:

    

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


    

In [8]:

    

wcs.contents


    

    
Out[8]:





{'topo': <owslib.coverage.wcs100.ContentMetadata at 0x7f2a342795d0>}

In [9]:

    

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


    

    




Topography

In [10]:

    

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


    

In [11]:

    

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


    

In [12]:

    

from osgeo import gdal
gdal.UseExceptions()
ds = gdal.Open('test.tif')


    

In [13]:

    

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

    

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


    

In [15]:

    

from cartopy.feature import NaturalEarthFeature
extent = [-71.5, -70.5, 42.0, 42.7]


coast = NaturalEarthFeature(category='physical', scale='10m',
                            facecolor='none', name='coastline')

fig, ax = make_map(projection=ccrs.PlateCarree())



ax.set_extent(extent)

feature = ax.add_feature(coast, edgecolor='gray')


    

In [15]:

    




    

In [ ]:

    

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