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
.
%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
NAT.shp
dataset from PySAL
examples into a GeoDataFrame
db = gpd.read_file(ps.examples.get_path('NAT.shp'))
NAD83
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
PySAL
's geotable
(note almost no performance hit, 1.13s for +3,000 polygons!)%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
(3085, 70)
%time geoplot(db)
CPU times: user 1.19 s, sys: 54.3 ms, total: 1.24 s Wall time: 1.12 s
pal = pltt.matplotlib.Viridis_12
geoplot(db, col='HR60', classi='Quantiles',
k=12, palette=pal)
geoplot(db, col='HR60', classi='equal_interval',
k=12, palette=pal)
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
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)
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)
This work is licensed under a Creative Commons Attribution 4.0 International License.