Legendgrams

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]:
%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")
%time db = db.to_crs(nad83)
CPU times: user 748 ms, sys: 8.96 ms, total: 757 ms
Wall time: 779 ms
  • Convert into PySAL's geotable (note almost no performance hit, 1.13s for +3,000 polygons!)
In [4]:
%time db = to_df(db)
db.shape
CPU times: user 1.19 s, sys: 22.1 ms, total: 1.21 s
Wall time: 1.48 s
Out[4]:
(3085, 70)
  • Get a feel for what it looks like (also very fast!)
In [5]:
%time geoplot(db)
CPU times: user 1.19 s, sys: 54.3 ms, total: 1.24 s
Wall time: 1.12 s
  • 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)