LISA cluster maps with PySAL¶
- Dani Arribas-Bel (@darribas)
This notebook is available as a gist, or viewable online as static html.
Local Indicators of Spatial Association (LISAs, Anselin, 1995) are a common tool to explore spatial heterogeneity, identify spatial concentration of (dis)similar values and test for the probability such agglomerations originate from a random process.
PySAL has had functionality to run LISAs since its very beginning. However, the output of such computation is not very useful just by itself as, being a local statistic, LISAs imply running a test for every single observation. The usual approach has been to visualize significant clusters on a map using a simple color coding for each class (High-High, HH; Low-Low, LL; High-Low, HL; Low-High, LH). Such visualiations have been available for a long time in packages such as GeoDa, but this would imply breaking the workflow to switch between Python and GeoDa, plus with the annoyance that the workflow cannot thus be automated, a convenient feature for reproducible science.
In this notebook, we show how to leverage the visualization layer that is being built in PySAL to create generic and custom LISA cluster maps.
Set up¶
For this example, we will use the NAT dataset, included in PySAL by default, exploring the variable HR90, homicide rates in 1990 at the county level.
Let's start by importing the required code. This is all code that will be made available in the upcoming release of PySAL in July 2015 (if you're reading this afterwards, it's all in there by default!).
%matplotlib inline
import pysal as ps
import numpy as np
from pysal.contrib.viz import mapping as maps
shp_link = ps.examples.get_path('NAT.shp')
print 'Reading from ', shp_link
hr90 = np.array(ps.open(shp_link.replace('.shp', '.dbf')).by_col('HR90'))
Reading from /Users/dani/code/pysal_darribas/pysal/examples/nat/NAT.shp
Before any further analysis, it is always good practice to visualize the distribution of values. We can now conveniently do that in PySAL:
_ = maps.plot_choropleth(shp_link, hr90, 'quantiles', cmap='Greens', figsize=(9, 6))
Before we can plot a map, we need to run a LISA. In the spirit of flexibility and modularity envisioned in PySAL, the two operations (computation, visualization) are decoupled. We will use a simple contiguity matrix for this.
w = ps.queen_from_shapefile(shp_link)
lisa = ps.Moran_Local(hr90, w, permutations=9999)
LISA maps in one line of code¶
The simplest way to obtain a quick visualization of LISA results is with the "one-liner" plot_lisa_cluster:
_ = maps.plot_lisa_cluster(shp_link, lisa, figsize=(9, 6))