#!/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
#
#
This work is licensed under a Creative Commons Attribution 4.0 International License.