#!/usr/bin/env python
# coding: utf-8

# # The Strava `heatmap` as contextual map tiles with `contextily`
# 
# This document quickly demonstrates how to source form the [Strava `heatmap`](https://labs.strava.com/heatmap/#7.00/-120.90000/38.36000/hot/all) to obtain map tiles that can be easily integrated in any modern Python geo-data workflow.

# In[13]:


# Display on the notebook
get_ipython().run_line_magic('matplotlib', 'inline')

# Import contextily
import contextily as ctx


# ## Basic example

# First we will set the link for the Strava tiles (taken from [here](https://qms.nextgis.com/geoservices/1447/)):

# In[10]:


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](https://github.com/darribas/contextily/blob/master/contextily_guide.ipynb)) to get the map delivered to your notebook. Let's play with Liverpool (UK), for example:

# In[11]:


lvl = ctx.Place('Liverpool', url=src)

ctx.plot_map(lvl.im);


# ## A more sophisticated illustration

# In[28]:


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](http://nbviewer.jupyter.org/gist/darribas/9a0d3b6177b7ca6be007/london_boroughs.ipynb) with Strava data. Here is what we will attempt to replicate:
# 
# ![](https://github.com/darribas/bits/raw/gh-pages/cards/boroughs.png)
# 
# You can download the original borough data from [here](https://data.london.gov.uk/dataset/statistical-gis-boundary-files-london) and a reprojected `GeoJSON` from here, which is what we will use:

# In[15]:


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):

# In[29]:


get_ipython().run_cell_magic('time', '', "raster_link = 'london.tiff'\nminX, minY, maxX, maxY = brs.to_crs(epsg=4326).total_bounds\n_ = ctx.bounds2raster(minX, minY, maxX, maxY, 12, raster_link, ll=True, url=src)\n")


# Just to get a sense, this is what the entire area of London looks like through the lens of Strava data:

# In[30]:


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).

# In[31]:


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:

# In[32]:


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:

# In[33]:


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:

# In[34]:


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:

# In[35]:


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:

# In[36]:


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()