WMTS (Web Map Tile Service) is an OGC (Open Geospatial Consortium) standard for requesting geo-located map tiles over HTTP. The OGC provide overview documents of the WMTS standard.
A WMTS request to a WMTS server contains layer (a map layer containing a given geographic phenomenon) and area of interest information. The information in the request defines the response from the WMTS server to which the request was made. The response is composed of geo-located map tiles that can be displayed by the application making the request.
The map tiles in a WMTS response are raster images, typically provided in PNG
or JPEG
format. Increased detail is available at higher zoom levels through a tree structure of tiles. Tiles are laid out in an (x, y)
grid to cover the spatial extent of the data with a z
axis for increased zoom levels.
An important part of OGC standards services is the get capabilities request. This returns an XML document that specifies the capabilities (metadata) about the service. Information included in the get capabilities response includes:
This notebook demonstrates using Cartopy to:
Let's demonstrate making an WMTS request with Cartopy and plotting the map image contained in the response. Let's use data from the NASA blue marble dataset exposed as a WMTS service.
To make a WMTS request in Cartopy, as with every WMTS request, we provide a URL to a WMTS server. This server provides a web map tile service, which may provide one or more layers of geographic phenomena. To complete the request we need to also provide the name of the layer we're interested in. To find the names of the available layers we need to refer to the get capabilities response from the WMTS server.
We can use the Python library OWSLib
to retrieve this information from the get capabilities response from the WMTS server:
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])
Request type: OGC WMTS Request title: NASA Global Imagery Browse Services for EOSDIS Number of available layers: 99 Example layers: ['MODIS_Aqua_Cloud_Top_Temp_Night', 'MODIS_Aqua_SurfaceReflectance_Bands721', 'MODIS_Combined_Value_Added_AOD', 'AMSRE_Brightness_Temp_89H_Day', 'MLS_HNO3_46hPa_Night', 'GMI_Brightness_Temp_Asc', 'GMI_Rain_Rate_Dsc', 'MODIS_Terra_Brightness_Temp_Band31_Day', 'MODIS_Terra_Cloud_Top_Pressure_Day', 'MODIS_Terra_Land_Surface_Temp_Day']
Let's take a more in-depth look at a layer of interest:
layer = layers[16]
print 'Layer title: {}'.format(wmts[layer].title)
print 'Bounding box: {}'.format(wmts[layer].boundingBoxWGS84)
Layer title: VIIRS_EarthAtNight_2012 Bounding box: (-180.0, -90.0, 180.0, 90.0)
With a bit of effort we can also find the CRS options for the layers provided by this WMTS server:
print [wmts.tilematrixsets[name].identifier for name in wmts.tilematrixsets]
['EPSG4326_250m', 'EPSG4326_2km', 'EPSG4326_16km', 'EPSG4326_31m', 'EPSG4326_500m', 'EPSG4326_1km']
In this case, the WMTS Server only provides tiles in EPSG-4326 (WGS-84), or Plate Carrée. Currently Cartopy can only plot WMTS tiles in a CRS that they are supplied in by the WMTS Server.
Making a WMTS request is completely handled by Cartopy behind the scenes. The above examples are just to demonstrate finding layer names and other useful information contained in a WMTS response.
With this information from the get capabilities to hand we can construct a WMTS request. Let's import some more useful libraries and define a convenience function that returns a common set of axes onto which to plot our WMTS layer and other data.
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
Now we can make a plot using map tiles from the NASA blue marble dataset:
proj = ccrs.PlateCarree()
ax = set_up_axes(proj, coast_colour='w')
ax.add_wmts(url, layer)
plt.show()
As an aside, web mapping services such as Google maps and Open Street Map (OSM) also use WMTS to serve the map tiles that are used to produce the map in an application. Cartopy has built-in map tile retrieval for Google tiles and OSM tiles. For example, to retrieve the OSM tiles for the USA:
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()
We can add value to our blue marble plot above by adding data to it. For example, we can add the location of major US cities to the plot and add air temperature data taken from an Iris cube.
Let's add the locations of three US cities to the plot, along with a text label of the city's name:
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()
The Iris sample data contains data on the predicted yearly average air temperature over the US according to the A1B scenario. Let's add the predicted air temperature for 2015 to the plot:
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()