adjacency.R

library(rgeos)
library(sp) ## the best package for handling shapefiles
library(spdep) ## for the function poly2nb
library(maptools) ## for the function readShapePoly
library(RColorBrewer) ## for pretty colour palettes
### Read in the shapefile
x <- readShapePoly("data/westminster_const_region.shp")

## Add the Press Association number
name2pa <- read.csv("data/name2pa.csv",header=T)
x$refno <- name2pa$refno[match(as.character(x$NAME),name2pa$ShapeFileName)]

# Create spatial polygons object
gpclibPermit()
sp.const <- unionSpatialPolygons(x, x$refno) ## may take a while
# Lose some detail
sp.const <- gSimplify(sp.const,100)
# Find out how many neighbours each constituency has
neighbouring.districts <- poly2nb(sp.const)

## Let's turn this into a data frame we can use
nb.df <- data.frame(refno = names(sp.const),
        numNeighbours = card(neighbouring.districts))
## Add some other info
nb.df <- merge(nb.df,name2pa,all.x=T,all.y=F)
## Check islands as a sanity check
nb.df[which(nb.df$numNeighbours==0),]

# Create spatial polygon data frame
# Use function SpatialPolygonsDataFrame from package "sp" (Pebesma 2005). 
# This matches external data in with geodata from sp.const:
# Rownames must match
rownames(nb.df) <- nb.df$refno
area.spdf <- SpatialPolygonsDataFrame(sp.const, nb.df)

# Now plot this
manual.colours <- cut(area.spdf$numNeighbours,
        breaks = quantile(area.spdf$numNeighbours,probs=seq(0,1,length.out=7)))
manual.colours.numeric <- as.numeric(manual.colours)
my.palette <- brewer.pal(7,name="PuBuGn")
manual.colours.hex <- my.palette[manual.colours.numeric]
plot(area.spdf, col = manual.colours.hex)