This is a write up of a set of functions written to create LISA maps - maps of Local Moran’s I values (or other spatial autocorrelation indicators if you so wish) in a non-standard way.

The main contribution of these maps is that both the Moran’s I values and their significances are plotted on a single map using poing shape/size in addition to colour palettes.

Additionally the maps incorporate a Moran’s scatterplot, which can of course be used separately.

A final bonus is a mapping function for a simple cloropleth map that has a density plot subset.

The code here is completely self contained and has been tested with the example data for England and Wales and some UK census 2001 variables for demonstration. (There is also a quick example at the end with a map of Germany, just to test it works).

The original .Rmd file for this, as well as a clean R file and the required data to run this example are all available at https://github.com/majazaloznik/VisualVignettes.

For comments, suggestions, issues and bugs: maja.zaloznik@gmail.com or report an Issue on the github repository.

First written @ Max Planck Institute for Demographic Research, February 2012. Presented @ AQMEN, Glasgow, November 2014 as part of the Visual Vignettes talk (repository in progress).

Outline

Data Setup

Palette Functions

Legend Functions

Moran Scatterplot

Maps

Bonus test example

First load all the required libraries:

library(spdep)
library(RColorBrewer)
library(classInt)
library(rgdal)
library(maptools)
gpclibPermit()
library(dplyr)
oldpar <- par(no.readonly=TRUE)
par(family = "serif", font = 5)

Data Setup

This example is demonstrated using England and Wales maps with two variables from the 2001 UK census Small area microdata collection. The shape file, a list of weights based on first order queen neighborhoods and the two variables are loaded here:

load("../Data/EnglandWales.R")

If you want to use your own data, simply swap it here as appropriate:

shape <- shape.LA
list.of.weights <- lw.foq
mydata <- data.frame(variable=distwork)

The next bit sets up the data frame that will be used throughout, by calculating:

lisa.data <-  mutate(mydata, 
                     lagged = lag.listw(list.of.weights,variable),
                     z = variable-mean(variable),
                     z = z/sqrt((sum(z^2)/n())),
                     laggedz = lag.listw(list.of.weights, z),
                     lisa = localmoran(variable, list.of.weights, alternative="two.sided")[,1],
                     sig = localmoran(variable, list.of.weights, alternative="two.sided")[,5],
                     quadrant = ifelse(z >= 0, ifelse(laggedz > 0,1,2), ifelse(laggedz >= 0,4,3)))

Palette Functions

We will be using several palettes in the Lisa Maps, the functions for which are defined as follows:

Regular Colour Palette - For cloropleth map

This palette takes the lisa.data dataframe prepared above. Other things that can be tweaked are self explanatory:

  • bin.n is the number of bins;
  • col.scheme is selected from the RColorBrewer palettes (use display.brewer.all() to see the options);
  • style is for classIntervals, choose from “fixed”, “sd”, “equal”, “pretty”, “quantile”, “kmeans”, “hclust”, “bclust”, “fisher”, or “jenks”;
  • z-score is logical, whether or not to use the raw variable or z-scores (this is only important for the legend, obviously doesn’t affect the colour palette)
FunPalette <- function (df=lisa.data, bin.n=6, col.scheme="YlOrRd", style="pretty", z.score=FALSE){
  if (z.score == TRUE) {df <- lisa.data$z} else {df <- lisa.data$variable} 
  bins <- classIntervals(df, n=bin.n, style=style)
  pal <- brewer.pal(length(bins$brks)-1, col.scheme)
  col <- findColours(bins, pal)
  pal.names <- paste(round(bins$brks[-length(bins$brks)],2), 
                     round(bins$brks[-1], 2), sep= " - ")
  pal.n.from <- round(bins$brks[-length(bins$brks)],2)
  pal.n.to <- round(bins$brks[-1],2)
  return(list(col, pal, pal.names, pal.n.from, pal.n.to, bins$brks, z.score))
}

Palette of point shapes and sizes - For significances in Moran map

This palette takes the lisa.data dataframe prepared above and returns point sizes and shapes for the map, as well as legend labels.

Simple adjustments can be made just by changing the scale factor, but of course if you want different significance bins, sizes and shapes these can be fixed inside the function as well.

FunSigSize <- function(df=lisa.data, scale=1) {
  sigz <- df[,6]
  bins <- c(0, 0.001, 0.005, 0.01, 0.05,  1 )
  cex.legend <- c(2.60, 2.34, 1.95, 1.30, 0.91)*scale
  pch.legend <- c(19, 19, 19, 19, 4)       # change this for different symbols
  cexez <- cex.legend[findInterval(sigz ,bins, rightmost.closed=T)]
  cex.names <- paste(round(bins[-length(bins)],4),"-",round(bins[-1],4))
  pchz <- pch.legend[findInterval(sigz ,bins, rightmost.closed=T)]
  return(list(cexez, cex.legend, cex.names, pchz, pch.legend,bins))
}

Palette for Local Moran’s I

This palette determines the colours for each area based on it’s Local Moran’s I value AND the quadrant it is in. It takes again the lisa.data dataframe prepared above and returns a colour vector as well as other legend stuff.

  • steps determines the number of colour steps within each subpalette;
  • four.way is logical which determines wheter each quadrant gets its own palette, or (default) whether positive values (HH & LL) get one and negative ones (HL & LH) get another;
  • pal.XX are RColorBrewer palettes (use display.brewer.all() to see the options);
FunPaletteMoran <- function (df=lisa.data,
                            steps = 3,
                            four.way = FALSE,
                            pal.hh = "PuBu", 
                            pal.hl = "Oranges",
                            pal.ll = "RdPu",
                            pal.lh = "Greens"){
  if (steps > 7) stop("Sorry, can't have more than seven steps !")
  if (four.way == FALSE) {pal.ll <- pal.hh 
                          pal.lh <- pal.hl}
  bins.pos <- seq(0, max(df$lisa), length=steps+1)
  bins.neg <- seq(min(df$lisa),0, length=steps+1)
  pal.hh <- brewer.pal(steps+2, pal.hh)[-c(1:2)]
  pal.hl <- rev(brewer.pal(steps+2, pal.hl)[-c(1:2)])
  pal.ll <- brewer.pal(steps+2, pal.ll)[-c(1:2)]
  pal.lh <- rev(brewer.pal(steps+2, pal.lh)[-c(1:2)])
  df <- df %>%
    mutate (int = ifelse(lisa >= 0, findInterval(lisa, bins.pos, rightmost.closed=TRUE), 
                         findInterval(lisa, bins.neg, rightmost.closed=TRUE)),
            colz = ifelse(quadrant==1, pal.hh[int], 
                          ifelse(quadrant == 2, pal.hl[int],
                                 ifelse(quadrant == 3, pal.ll[int], pal.lh[int]))))
  pal.n.from <- c(round(bins.neg[-length(bins.neg)],2),round(bins.pos[-length(bins.pos)],2) )
  pal.n.to <- c(round( bins.neg[-1],2), round( bins.pos[-1],2))
  pal.names <- paste( pal.n.from, pal.n.to,sep= " - ")                        
  return(list(df$colz, pal.hh, pal.hl, pal.ll, pal.lh, pal.names))
}

Legend Functions

Although the legends can be generally done using the legend() function, the moran one, where you have four colour palettes is a bit trickier. Also, different positioning of legends on different devices can be tricky, so instead I’ve done these helper functions for legends, which can be placed in subplots using par(fig).

FunLegend <- function(palette=FunPalette){
  plot(c(0,1), c(0,1),  col = NA, ann = FALSE, axes = FALSE)
  points(rep(0.3,length(palette[[2]])), 1/length(palette[[2]])*
           seq(1, length(palette[[2]])),pch=15, col=palette[[2]], cex=2)   
  points(rep(0.3,length(palette[[2]])),1/length(palette[[2]])*
           seq(1, length(palette[[2]])),pch=0, lwd=1, cex=2)
  text(rep(0.6,length(palette[[2]])), 1/length(palette[[2]])*
         seq(1, length(palette[[2]])), palette[[3]])
}

FunLegendSig <- function(sizes=FunSigSize()){
  par(mar=c(0,0,2,0)) 
  plot(c(0,1), c(0,1),  col = NA, ann = FALSE, axes = FALSE)
  title(main="Significance Levels", font.main=1, cex.main=1)
  points(rep(0.2,length(sizes[[2]])), 1/length(sizes[[2]])*
           seq(1, length(sizes[[2]])),pch=sizes[[5]], cex=sizes[[2]], col="gray", lwd=2)   
  points(rep(0.2,length(sizes[[2]])), 1/length(sizes[[2]])*
           seq(1, length(sizes[[2]])),  pch=c(rep(1,length(sizes[[5]])-1), NA),
         cex=sizes[[2]][sizes[[5]]!=4], lwd=0.5)
  text(rep(0.7,length(sizes[[2]])), 1/length(sizes[[2]])*
         seq(1, length(sizes[[2]])), sizes[[3]])}

FunLegendMoran <- function(palette=FunPaletteMoran()){
  par(mar=c(0,0,2,0)) 
  plot(c(0,1), c(0,1),  col = NA, ann = FALSE, axes = FALSE)
  title(main="Local Moran's I", font.main=1, cex.main=1)
  points(rep(0.3,2*length(palette[[2]])), 1/(2*length(palette[[2]]))*
           seq(1, 2*length(palette[[2]])),pch=19,  col=c(palette[[3]], palette[[2]]), cex=1.5)
  points(rep(0.3,2*length(palette[[2]])), 1/(2*length(palette[[2]]))*
           seq(1, 2*length(palette[[2]])),  cex=1.5)
  if (!identical(palette[[3]], palette[[5]])){
  points(rep(0.2,2*length(palette[[2]])), 1/(2*length(palette[[2]]))*
           seq(1, 2*length(palette[[2]])),pch=19,  col=c(palette[[5]], palette[[4]]), cex=1.5)
  points(rep(0.2,2*length(palette[[2]])), 1/(2*length(palette[[2]]))*
           seq(1, 2*length(palette[[2]])),  cex=1.5)}
  text(rep(0.7,2*length(palette[[2]])), 1/(2*length(palette[[2]]))*
         seq(1, 2*length(palette[[2]])), palette[[6]])
}

Moran Scatterplot

Can be used independently, or called from FunLisaMap (see bellow). It takes again the lisa.data dataframe prepared above, as well as:

ToDo: Right now you can manually change inside the function if you want either:

Vaule: For plotting legend convenience, the function returns the palettes used.

FunMoranScatterplot <- function(df=lisa.data, z.score=FALSE, 
                                palette= FunPaletteMoran(),
                                sizes = FunSigSize(),  ...) {  
  if (z.score == TRUE) {x <- df$z
                        lx <- df$laggedz
                        xlab <- "z-score"
                        ylab <- "lagged z-score"} else {
                          x <- df$variable
                          lx <- df$lagged
                          xlab <- "variable"
                          ylab <- "lagged variable"} 
  plot(x, lx, col=palette[[1]], cex=sizes[[1]], pch=sizes[[4]], xlab=xlab, ylab=ylab, ...)
  points(x[sizes[[4]]!=4], lx[sizes[[4]]!=4], cex=sizes[[1]][sizes[[4]]!=4])
  abline(v= mean(x), lty=2, lwd=0.5)
  ## option1? x=lx and mean(x) on both
  abline(a=0, b=1, lwd=0.5) 
  abline(h= mean(x), lty=2, lwd=0.5)  
  ## option2? linear model and mean of lagged
  #abline(lm(lx~x), lwd=0.5) 
  #abline(h= mean(lx), lty=2, lwd=0.5)
  return(list(palette, sizes))
}

So here’s an example of the Moran scatterplot, using the four.way option as well as z-scores and with the associated legends along side. The variable is the percent of people with distance to work < 5km.

oldpar <- par(no.readonly=TRUE)
par(mar=c(5.1, 4.1, 4.1, 9.1))
lgnd <- FunMoranScatterplot(z.score=TRUE, palette=FunPaletteMoran(four.way=TRUE, steps=4) )
leftX <- 0.80 # bottom legend
rightX <- 0.99
bottomY <- 0.1
topy <- 0.5
par(mar=c(5.1, 4.1, 4.1, 2.1))
par(fig=c(leftX,rightX,bottomY, topy), new =TRUE)
FunLegendMoran(lgnd[[1]])
bottomY <- 0.5 # top legend
topy <- 0.9
par(fig=c(leftX,rightX,bottomY, topy), new =TRUE)
FunLegendSig(lgnd[[2]])

par(oldpar)

Maps

Regular cloropleth map of the variable (or z-scores), with inset, density plot and legend

Takes again the lisa.data dataframe prepared above, as well as a spatialPolygon object and a regular palette from FunPalette.

The default inset is for London. The centering and size of the inset can be changed in the code:

  • CoL is the center of the inset box
  • h and w are its height and width - aspect ratio is the same as the overall plot. So change the value 30 to another one for larger or smaller coverage.
  • leftX, rightX, bottomY are the positions of the zoomed inset. Again this can be moved around, and it set to have the same aspect ratio as the overall plot.

The same principles aply to changing the size and position of the density distribution plot and the legend.

FunCloroplethMap <- function(df=lisa.data, shp=shape, palette=FunPalette()){
  oldpar <- par(no.readonly=TRUE)
  # X11(width=10, height=7)  
  par(mar=c(0,0,0,0)) 
  plot(shp, col=palette[[1]])
  par("usr")
  # London inset, change CoL and w1, and h1 to suit
  CoL <- coordinates(shp)[32,] # center on London City (00BK)
  w1 = abs(par("usr")[1] - par("usr")[2])/30
  h1 = abs(par("usr")[3] - par("usr")[4])/30
  rect( CoL[1]-w1, CoL[2]-h1,CoL[1]+w1, CoL[2]+h1, lwd=2)
  leftX <- 0.78 # position London Inset
  rightX <- 0.99
  bottomY <- 0.2
  par(fig=c(leftX,rightX,bottomY, bottomY+rightX-leftX), new =TRUE)
  plot(shp, col=palette[[1]], xlim=c(CoL[1]-w1, CoL[1]+w1),
       ylim=c(CoL[2]-h1, CoL[2]+h1))
  box(lwd=2)  
  leftX <- 0 # position variable density distribution
  rightX <- 0.25
  bottomY <- 0.35
  par(fig=c(leftX,rightX,bottomY, bottomY+rightX-leftX), new =TRUE,mar=c(2,1,0,0))
  if (palette[[7]] == TRUE) {x <- df$z} else {x <- df$variable} 
  plot(density(x, na.rm=TRUE), axes=FALSE, ann=FALSE, lwd=2)
  axis(1)
  abline(v=palette[[6]], col="red", lty=2) 
  bottomY <- bottomY+rightX-leftX # position legend above dens.
  par(fig=c(leftX,rightX,bottomY, bottomY+rightX-leftX),mar=c(0,0,0,0), new =TRUE)
  FunLegend(palette)
  par(oldpar)
}

An example, where the only default changed is the number of bins:

FunCloroplethMap(palette=FunPalette(bin.n=6))

Moran map, with inset, moran scatterplot and legends

Takes again the lisa.data dataframe prepared above, as well as a spatialPolygon object. Then the special pallete and size/shape palette from the funcitons FunPaletteMoran and FunSigSize.

The default inset is for London. The centering and size of the inset can be changed in the code:

  • `CoL is the center of the inset box
  • h and w are its height and width - aspect ratio is the same as the overall plot. So change the value 30 to another one for larger or smaller coverage.
  • leftX, rightX, bottomY are the positions of the zoomed inset. Again this can be moved around, and it set to have the same aspect ratio as the overall plot.

The same principles aply to changing the size and position of the two legends and the scatterplot.

FunMoranMap <- function(df=lisa.data, shp=shape, palette=FunPaletteMoran(),
                          sizes=FunSigSize(), z.score=FALSE){
  oldpar <- par(no.readonly=TRUE)
  #X11(width=10, height=7)  
  par(mar=c(0,0,0,0)) 
  plot(shp,  lwd=0.5)
  points(coordinates(shp), col=palette[[1]], 
         cex=sizes[[1]],pch=sizes[[4]], lwd=2 )
  points(coordinates(shp)[sizes[[4]]!=4,], cex=sizes[[1]][sizes[[4]]!=4] )
  par("usr")
  # London inset, change CoL and w1, and h1 to suit
  CoL <- coordinates(shp)[32,] # center on London City (00BK)
  w1 = abs(par("usr")[1] - par("usr")[2])/30
  h1 = abs(par("usr")[3] - par("usr")[4])/30
  rect( CoL[1]-w1, CoL[2]-h1,CoL[1]+w1, CoL[2]+h1, lwd=2)
  leftX <- 0.78 # position London Inset
  rightX <- 0.99
  bottomY <- 0.2
  par(fig=c(leftX,rightX,bottomY, bottomY+rightX-leftX), new =TRUE)
  plot(shape,  lwd=0.5,xlim=c(CoL[1]-w1, CoL[1]+w1),
       ylim=c(CoL[2]-h1, CoL[2]+h1))
  points(coordinates(shp), col=palette[[1]], 
         cex=sizes[[1]],pch=sizes[[4]], lwd=2 )
  points(coordinates(shp)[sizes[[4]]!=4,], cex=sizes[[1]][sizes[[4]]!=4] )
  box(lwd=2)  
  leftX <- 0 # position legend significances
  rightX <- 0.34
  bottomY <- 0.3
  par(fig=c(leftX,rightX,bottomY, bottomY+rightX-leftX), new =TRUE,mar=c(2,1,0,0))
  FunLegendSig(sizes)
  bottomY <- bottomY+rightX-leftX# position legend Morans I
  par(fig=c(leftX,rightX,bottomY, bottomY+rightX-leftX), new =TRUE,mar=c(2,1,0,0))
  FunLegendMoran(palette)
  par(oldpar)
  leftX <- 0.70 # position moran scatterplot
  rightX <- 0.99
  bottomY <- 0.65
  par(fig=c(leftX,rightX,bottomY, bottomY+rightX-leftX), new =TRUE,mar=c(0,0,0,0))
  FunMoranScatterplot(df, palette=palette, sizes=sizes, z.score=z.score)
  par(oldpar) 
}

Let’s test out the default version:

FunMoranMap()

OK, lets try out some other options:

  • another variable - proportion of people two or more cars
  • four way palette,
  • z-scores in the scatteplot,
  • scale down the point sizes a bit
mydata <- data.frame(variable=cars.data)
lisa.data <-  mutate(mydata, 
                     lagged = lag.listw(list.of.weights,variable),
                     z = variable-mean(variable),
                     z = z/sqrt((sum(z^2)/n())),
                     laggedz = lag.listw(list.of.weights, z),
                     lisa = localmoran(variable, list.of.weights, alternative="two.sided")[,1],
                     sig = localmoran(variable, list.of.weights, alternative="two.sided")[,5],
                     quadrant = ifelse(z >= 0, ifelse(laggedz > 0,1,2), ifelse(laggedz >= 0,4,3)))

FunMoranMap(palette=FunPaletteMoran(four.way=TRUE ), sizes=FunSigSize(scale=0.8), z.score=TRUE)

Bonus Test Example

Just a quick test with some other data. Downloadint the shape file for Germany, quick fix of neighbourhood list, weights, and just a random normal variable to plot:

temp <- tempfile()
download.file("http://biogeo.ucdavis.edu/data/diva/adm/DEU_adm.zip",temp)
unzip(temp)
shape.DEU <- readOGR("DEU_adm3.shp", layer="DEU_adm3")
## OGR data source with driver: ESRI Shapefile 
## Source: "DEU_adm3.shp", layer: "DEU_adm3"
## with 434 features and 19 fields
## Feature type: wkbPolygon with 2 dimensions
unlink(temp)

##create neighbourhood list and fix Ruegen
shape.nb <- poly2nb(shape.DEU)
which(card(shape.nb)==0)
## [1] 250
shape.nb[[250]] <- as.integer(245)
shape.nb[245] <- list(sort(c(shape.nb[[245]], as.integer(250))))

## create list of weights
lw.DEU <- nb2listw(shape.nb)

## Setup data
shape <- shape.DEU
list.of.weights <- lw.DEU
set.seed(42)
mydata <- data.frame(variable= rnorm(length(shape)))

lisa.data <-  mutate(mydata, 
                     lagged = lag.listw(list.of.weights,variable),
                     z = variable-mean(variable),
                     z = z/sqrt((sum(z^2)/n())),
                     laggedz = lag.listw(list.of.weights, z),
                     lisa = localmoran(variable, list.of.weights, alternative="two.sided")[,1],
                     sig = localmoran(variable, list.of.weights, alternative="two.sided")[,5],
                     quadrant = ifelse(z >= 0, ifelse(laggedz > 0,1,2), ifelse(laggedz >= 0,4,3)))

FunCloroplethMap(palette=FunPalette(col.scheme="YlGnBu"))

OK, that seems to work fine, obviously I should comment out the part of the function for the inset, or change it to something more useful (ToDo).

And now for the Moran Plot, this time with some different palettes:

FunMoranMap(palette=FunPaletteMoran(four.way=TRUE, steps=4, pal.hh = "RdBu", 
                            pal.hl = "Greys",
                            pal.ll = "BuPu",
                            pal.lh = "YlOrBr"), sizes=FunSigSize(scale=0.8), z.score=TRUE)