Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment