contextily
¶%matplotlib inline
import contextily as ctx
import geopandas as gpd
import rasterio as rio
import numpy as np
import matplotlib.pyplot as plt
# Data
from pysal.examples import get_path
In this notebook, we will show the basic functionality available in contextily
, a package to work with web-tiles for background maps. To do that, we will use additional data to illustrate contextily
can be integrated with other libraries such as geopandas
and rasterio
. However, it is important to note that some of those are NOT required by contextily
.
Let us first load up from PySAL
the polygons of the 48 contiguous states using geopandas
, and pick out Texas, for instance:
tx = gpd.read_file(get_path('us48.shp')).set_index('STATE_ABBR').loc[['TX'], 'geometry']
tx.crs = {'init': 'epsg:4326'}
tx['TX']
We can use this polygon to get its bounding box:
w, s, e, n = tx['TX'].bounds
w, s, e, n
(-106.6495132446289, 25.845197677612305, -93.50721740722656, 36.49387741088867)
At this point, we can use those bounds to download the tiles for that part of the world. Before that, however, it is convenient to make sure the download will not be too heavy (i.e. not too many tiles). This is easily checked by using the howmany
utility:
_ = ctx.howmany(w, s, e, n, 6, ll=True)
Using zoom level 6, this will download 9 tiles
At the zoom level 6, we need to download 9 tiles, which is not too onerous. Let us then go ahead and download them:
%time img, ext = ctx.bounds2img(w, s, e, n, 6, ll=True)
CPU times: user 32 ms, sys: 8 ms, total: 40 ms Wall time: 2.53 s
That's it! Just under three seconds and we have a map under our fingertips! Note how bounds2img
also returns the extent that the tile covers. This will come in handy later on when we want to align it with more data.
IMPORTANT: the tile extent is always returned in its original CRS, the Spherical Mercator (EPSG:3857).
Let us quickly visualize it with matplotlib
:
plt.imshow(img, extent=ext);
Sometimes, we know we will be working with an area for a while and it is more convenient to download the tiles only once and then load them up on-demand. To do that, we can use the function bounds2raster
, which writes the image into a GeoTIFF raster file:
%time _ = ctx.bounds2raster(w, s, e, n, 6, 'tx.tif', ll=True)
CPU times: user 52 ms, sys: 4 ms, total: 56 ms Wall time: 2.69 s
At this point, using the tile map is exactly the same as using any other raster file. For this, we will use the fantastic library rasterio
.
Let us see how we can load it back up and plot it (all of this is standard rasterio
operations, check out its documentation for more detail if interested).
rtr = rio.open('tx.tif')
# NOTE the transpose of the image data
img = np.array([ band for band in rtr.read() ]).transpose(1, 2, 0)
# Plot
plt.imshow(img, extent=rtr.bounds);
Sometimes, raster files can be (very) large and we are not particularly interested in loading the entire file. For those cases, rasterio
includes functionality to load up a specific part -window- of the raster. The function bb2wdw
in contextily
makes it easy to create that window from Mercator coordinates:
# Mercator coordinates for Houston area
hou = (-10676650.69219051, 3441477.046670125, -10576977.7804825, 3523606.146650609)
# Window
wdw = ctx.bb2wdw(hou, rtr)
# Raster subset
sub = np.array([ rtr.read(band, window=wdw) \
for band in range(1, rtr.count+1)]).transpose(1, 2, 0)
# Plot
plt.imshow(sub, extent=(hou[0], hou[2], hou[1], hou[3]));
One of the most interesting applications of using these tiles is to employ them as basemaps to overlay additional data on top of them. This can be easily done with matplotlib
, provided all the data are in the same CRS at the time of plotting.
Let us see how one would go about plotting the polygon for the state of Texas on top of the raster file:
# Shortify the bound box named tuple
bb = rtr.bounds
# Set up the figure
f, ax = plt.subplots(1, figsize=(9, 9))
# Load the tile raster (note the re-arrangement of the bounds)
ax.imshow(img, extent=(bb.left, bb.right, bb.bottom, bb.top))
# Overlay the polygon on top (note we reproject it to the raster's CRS)
tx.to_crs(rtr.crs).plot(edgecolor='none', ax=ax)
# Remove axis for aesthetics
ax.set_axis_off()
# Show
plt.show()
sources = [i for i in dir(ctx.tile_providers) if i[0]!='_']
sources
['ST_TERRAIN', 'ST_TERRAIN_BACKGROUND', 'ST_TERRAIN_LABELS', 'ST_TERRAIN_LINES', 'ST_TONER', 'ST_TONER_BACKGROUND', 'ST_TONER_HYBRID', 'ST_TONER_LINES', 'ST_TONER_LITE', 'ST_WATERCOLOR']
You can set them on bounds2img
and bounds2raster
using the argument url
. Checkout the documentation for more details.
Just because we can, let us get a feel for what they look like:
f, axs = plt.subplots(2, 5, figsize=(25, 10))
axs = axs.flatten()
for src, ax in zip(sources, axs):
img, ext = ctx.bounds2img(w, s, e, n, 6, url=getattr(ctx.sources, src), ll=True)
ax.imshow(img, extent=ext)
ax.set_title(src)
ax.set_axis_off()
plt.show()
NOTE Please always remember to give proper attribution to the map provider. See here for the proper way to do it, but essentially it is:
Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL.
Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under CC BY SA.