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

# # Legendgrams
# 
# * [Dani Arribas-Bel](http://darribas.org) [`[@darribas]`](http://twitter.com/darribas)
# 
# This document present Python code to generate choropleth legends in the form of histograms. In passing, it also shows how easy and quick plotting geographic information is with `PySAL`'s new `geoplot` and how `PySAL`'s data structures can be combined and interfaced with those in `geopandas`.
# 

# In[1]:


get_ipython().run_line_magic('matplotlib', 'inline')

import matplotlib.pyplot as plt
import palettable as pltt
import numpy as np
import pysal as ps
import geopandas as gpd
from pysal.contrib import pdio
from pysal.contrib.geotable.utils import to_df, to_gdf
from pysal.contrib.viz.mapping import geoplot


# * Read `NAT.shp` dataset from `PySAL` examples into a `GeoDataFrame`

# In[2]:


db = gpd.read_file(ps.examples.get_path('NAT.shp'))


# * Project into `NAD83`

# In[3]:


db.crs = {'init':'epsg:4326'}
nad83 = ("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 "\
         "+lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs")
get_ipython().run_line_magic('time', 'db = db.to_crs(nad83)')


# * Convert into `PySAL`'s `geotable` (note almost no performance hit, 1.13s for +3,000 polygons!)

# In[4]:


get_ipython().run_line_magic('time', 'db = to_df(db)')
db.shape


# * Get a feel for what it looks like (also very fast!)

# In[5]:


get_ipython().run_line_magic('time', 'geoplot(db)')


# * Select Viridis as the color palette

# In[6]:


pal = pltt.matplotlib.Viridis_12


# * Quantile map

# In[7]:


geoplot(db, col='HR60', classi='Quantiles', 
        k=12, palette=pal)


# * Equal interval map

# In[8]:


geoplot(db, col='HR60', classi='equal_interval', 
        k=12, palette=pal)


# * Add fancy legend/histogram (legendgram???)

# In[9]:


def legendgram(f, ax, y, breaks, pal, rescale=False):
    '''
    Add a histogram in a choropleth with colors aligned with map
    ...
    
    Arguments
    ---------
    f           : Figure
    ax          : AxesSubplot
    y           : ndarray/Series
                  Values to map
    breaks      : list
                  Sequence with breaks for each class (i.e. boundary values
                  for colors)
    rescale     : False/tuple
                  [Optional. Default=False] If a tuple, limits to set the X
                  axis of the histogram
    '''
    k = len(breaks)
    p = ax.get_position()
    histpos = [p.x0 + p.width*0.01, \
            p.y0 + p.height*0.1,  \
            p.width * 0.27, \
            p.height * 0.2]
    histax = f.add_axes(histpos)
    N, bins, patches = histax.hist(y, bins=50, color='0.1')
    #---
    pl = pal.get_mpl_colormap()
    bucket_breaks = [0]+[np.searchsorted(bins, i) for i in breaks]
    for c in range(k):
        for b in range(bucket_breaks[c], bucket_breaks[c+1]):
            patches[b].set_facecolor(pl(c/k))
    #---
    if rescale:
        histax.set_xlim(rescale)
    histax.set_frame_on(False)
    histax.get_yaxis().set_visible(False)
    histax.tick_params(labelsize=12)
    return None


# In[10]:


y = db['HR60']
f, ax = plt.subplots(1, figsize=(12, 9))
geoplot(db, col='HR60', classi='Quantiles', 
        k=12, palette=pal, ax=ax)
breaks = ps.Quantiles(y, 12).bins
legendgram(f, ax, y, breaks, pal)


# In[11]:


f, ax = plt.subplots(1, figsize=(12, 12))
geoplot(db, col='HR60', classi='equal_interval', 
        k=12, palette=pal, ax=ax)
breaks = ps.Equal_Interval(y, 12).bins
legendgram(f, ax, y, breaks, pal)


# 
# # License
# 
# <a rel="license" href="http://creativecommons.org/licenses/by/4.0/"><img alt="Creative Commons License" style="border-width:0" src="https://i.creativecommons.org/l/by/4.0/88x31.png" /></a><br />This work is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by/4.0/">Creative Commons Attribution 4.0 International License</a>.