from owslib.wmts import WebMapTileService

url = 'http://map1c.vis.earthdata.nasa.gov/wmts-geo/wmts.cgi'
wmts = WebMapTileService(url)
layers = list(wmts.contents)
print 'Request type: {}'.format(wmts.identification.type)
print 'Request title: {}'.format(wmts.identification.title)
print 'Number of available layers: {}'.format(len(layers))
print 'Example layers: {}'.format(layers[:10])

layer = layers[16]

print 'Layer title: {}'.format(wmts[layer].title)
print 'Bounding box: {}'.format(wmts[layer].boundingBoxWGS84)

print [wmts.tilematrixsets[name].identifier for name in wmts.tilematrixsets]

import cartopy.crs as ccrs
import matplotlib.pyplot as plt

us_extent = [-126.0, -66.0, 24.0, 50.0]

def set_up_axes(proj, coast_colour=None, extent=us_extent):
    fig = plt.figure(figsize=(14, 9))
    ax = plt.axes(projection=proj)
    if coast_colour is not None:
        ax.coastlines(color=coast_colour)
    ax.set_extent(extent, crs=ccrs.PlateCarree())
    return ax

proj = ccrs.PlateCarree()
ax = set_up_axes(proj, coast_colour='w')

ax.add_wmts(url, layer)

plt.show()

from cartopy.io.img_tiles import OSM

osm_tiles = OSM()
proj = osm_tiles.crs
ax = set_up_axes(proj)

ax.add_image(osm_tiles, 5)

plt.show()

import matplotlib.patheffects as mpath

# Set up the axes and add a WMTS layer, as above.
proj = ccrs.PlateCarree()
ax = set_up_axes(proj, coast_colour='w')

ax.add_wmts(url, layer)

# Add the locations of our three cities.
cities = {'Washington DC': [38.928, -76.981],
          'Austin': [30.308, -97.753],
          'San Francisco': [37.758, -122.438],}
city_red = '#b80000'
star_size = 75

for name, latlon in cities.iteritems():
    plt.scatter(latlon[1], latlon[0], c=city_red, s=star_size,
                linewidths=0, marker='*', zorder=5)
    plt.text(latlon[1]-1.0, latlon[0]+1.0, name, color='w', size=12,
             path_effects=[mpath.withStroke(linewidth=2, foreground='k')])

plt.show()

import iris
import iris.quickplot as qplt

year_2015 = iris.Constraint(time=lambda cell: cell.point.year == 2015)
with iris.FUTURE.context(cell_datetime_objects=True):
    cube = iris.load_cube(iris.sample_data_path('A1B_north_america.nc'), year_2015)

# Set up the axes, add a WMTS layer and city locations; all as above.
proj = ccrs.PlateCarree()
ax = set_up_axes(proj, coast_colour='w')

ax.add_wmts(url, layer)

for name, latlon in cities.iteritems():
    plt.scatter(latlon[1], latlon[0], c=city_red, s=star_size,
                linewidths=0, marker='*', zorder=5)
    plt.text(latlon[1]-1.0, latlon[0]+1.0, name, color='w', size=12,
             path_effects=[mpath.withStroke(linewidth=2, foreground='k')])
    
# Add the average air temperature data to the plot as contours.
nlevels = 12
mappable = qplt.contour(cube, nlevels)
plt.clabel(mappable, inline=1, fontsize=12)

plt.show()