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

# # `contextily`

# In[2]:


get_ipython().run_line_magic('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:

# In[4]:


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:

# In[6]:


w, s, e, n = tx['TX'].bounds
w, s, e, n


# ## Download tiles
# 
# 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:

# In[8]:


_ = ctx.howmany(w, s, e, n, 6, ll=True)


# At the zoom level 6, we need to download 9 tiles, which is not too onerous. Let us then go ahead and download them:

# In[8]:


get_ipython().run_line_magic('time', 'img, ext = ctx.bounds2img(w, s, e, n, 6, ll=True)')


# 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`:

# In[9]:


plt.imshow(img, extent=ext);


# ## Save tiles into raster files
# 
# 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:

# In[10]:


get_ipython().run_line_magic('time', "_ = ctx.bounds2raster(w, s, e, n, 6, 'tx.tif', ll=True)")


# ## Read raster tiles and combine with vector data
# 
# 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](https://mapbox.github.io/rasterio/) for more detail if interested).

# In[24]:


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:

# In[12]:


# 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]));


# ## Combine layers
# 
# 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:

# In[29]:


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


# ## Different backgrounds
# 
# `contextily` gives access to several tile maps, all from the awesome people at [Stamen](http://stamen.com). The full list is:

# In[32]:


sources = [i for i in dir(ctx.tile_providers) if i[0]!='_']
sources


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

# In[35]:


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](http://maps.stamen.com/#watercolor/12/37.7706/-122.3782) for the proper way to do it, but essentially it is:
# 
# * Toner and Terrain: 
# 
# > Map tiles by <a href="http://stamen.com">Stamen Design</a>, under <a href="http://creativecommons.org/licenses/by/3.0">CC BY 3.0</a>. Data by <a href="http://openstreetmap.org">OpenStreetMap</a>, under <a href="http://www.openstreetmap.org/copyright">ODbL</a>.
# 
# * Watercolor:
# 
# > Map tiles by <a href="http://stamen.com">Stamen Design</a>, under <a href="http://creativecommons.org/licenses/by/3.0">CC BY 3.0</a>. Data by <a href="http://openstreetmap.org">OpenStreetMap</a>, under <a href="http://creativecommons.org/licenses/by-sa/3.0">CC BY SA</a>.