heatmap
as contextual map tiles with contextily
¶This document quickly demonstrates how to source form the Strava heatmap
to obtain map tiles that can be easily integrated in any modern Python geo-data workflow.
# Display on the notebook
%matplotlib inline
# Import contextily
import contextily as ctx
First we will set the link for the Strava tiles (taken from here):
src = 'https://heatmap-external-a.strava.com/tiles/all/hot/tileZ/tileX/tileY.png'
Then it's all a matter of using the Place
API (or any other way that contextily
allows) to get the map delivered to your notebook. Let's play with Liverpool (UK), for example:
lvl = ctx.Place('Liverpool', url=src)
ctx.plot_map(lvl.im);
import geopandas as gpd
import numpy as np
import matplotlib.pyplot as plt
import rasterio as rio
from rasterio import mask
from rasterio.plot import show as rioshow
from shapely.geometry import mapping as shply2json
src = 'https://heatmap-external-a.strava.com/tiles/all/bluered/tileZ/tileX/tileY.png'
Because contextily
allows you to combine webtile maps with any other data you may have, you can easily build more sophisticated maps. In this case, we will recreate the London boroughs example with Strava data. Here is what we will attempt to replicate:
You can download the original borough data from here and a reprojected GeoJSON
from here, which is what we will use:
brs = gpd.read_file('boroughs.geojson')
brs.plot();
In order to render the images faster and not have to rely on the remote server to pull the tiles, we will first download them all at once and store them as a tiff
raster file (keep in mind this might take a little bit to run):
%%time
raster_link = 'london.tiff'
minX, minY, maxX, maxY = brs.to_crs(epsg=4326).total_bounds
_ = ctx.bounds2raster(minX, minY, maxX, maxY, 12, raster_link, ll=True, url=src)
CPU times: user 2.69 s, sys: 292 ms, total: 2.98 s Wall time: 41.5 s
Just to get a sense, this is what the entire area of London looks like through the lens of Strava data:
r_src = rio.open(raster_link)
f, ax = plt.subplots(1, figsize=(12, 12))
rioshow(r_src.read(), ax=ax)
ax.set_axis_off()
plt.show()
Now all we need to do is to be able to load the data in raster file for the geometry of each borough. The following function does exactly that (plus it colors empty values in white for aesthetic purposes).
def crop_borough(brg, raster='london.tiff'):
'''
Crop borough expressed as Shapely polygon (`brg`) from `raster`
and return as image array
...
Arguments
---------
brg : Polygon
Shapely polygon to be cropped
raster : str
[Optional. Default='london.tiff'] Path to raster file
containing image data to be cropped
'''
mi, mt = mask.mask(rio.open(raster), \
[shply2json(brg)], crop=True)
mi.data[np.where(mi.data==0)] = 255
img = mi.data.swapaxes(2, 0).swapaxes(1, 0)
return img
As an example, if we want to see what the City of London looks like, we can just run:
plt.imshow(
crop_borough(
brs.set_index('NAME')\
.loc['City of London', 'geometry']
));
Now we are in a position where we can replicate the original plot with Strava tiles:
f, axs = plt.subplots(6, 6, figsize=(24, 24))
axs = axs.flatten()
for i in range(brs.shape[0]):
ax = axs[i]
ax.imshow(
crop_borough(
brs.loc[i, 'geometry']
))
for i in range(36):
axs[i].set_axis_off()
#plt.savefig('/home/dani/Desktop/boroughs.png')
plt.show()
Can you tell which is which? In case you need further cues to identify each borough, here's a version that includes a small map of London highlighting the borough displayed! To do this, we'll also write a small function that adds the map first:
uno = gpd.GeoSeries(brs.unary_union)
def add_london(ax, poly_id, brs):
p = ax.get_position()
mapos = [p.x0 + p.width*0.05, \
p.y0 + p.height*0.05, \
p.width * 0.27, \
p.height * 0.2]
mapax = f.add_axes(mapos)
uno.plot(ax=mapax, color='k', alpha=0.25)
brs.loc[[poly_id], :].plot(ax=mapax, color='#1BFF1B')
mapax.set_axis_off()
mapax.axis('equal')
return mapax
See what it would look like for the City of London:
f, ax = plt.subplots(1)
ax.imshow(
crop_borough(
brs.set_index('NAME')\
.loc['City of London', 'geometry']
))
_ = add_london(ax, 32, brs)
plt.show()
And roll it out to the entire set:
f, axs = plt.subplots(6, 6, figsize=(24, 24))
axs = axs.flatten()
for i in range(brs.shape[0]):
ax = axs[i]
ax.imshow(
crop_borough(
brs.loc[i, 'geometry']
))
_ = add_london(ax, i, brs)
for i in range(36):
axs[i].set_axis_off()
plt.show()