R Choropleth with shape file

rowlingson.R

###
### Barry Rowlingson, Lancaster University
###

## needed for shapefiles:
require(rgdal)

## needed for colour mapping - not on CRAN:
## http://r-forge.r-project.org/projects/colourscheme/
## try:
## install.packages("colourschemes", repos="http://R-Forge.R-project.org")
require(colourschemes)

## pretty colours:
require(RColorBrewer)

## read the unemployment data:
unem=read.table("http://datasets.flowingdata.com/unemployment09.csv",sep=",",
  as.is=TRUE,
  colClass=c(rep("character",8),"numeric"),
  col.names=c("code","f1","f2","name","yr","x1","x2","x3","percentage")
  )

## compute the FIPS code:
unem$fips=paste(unem$f1,unem$f2,sep="")

## read the county data:
## source: http://www.census.gov/geo/www/cob/co2000.html
## this is in lat-long, with Alaska, Hawaii, and Puerto Rico
##   in their proper geographic place. So I'm just going to
##   plot the lower 48

download.file("http://www2.census.gov/geo/tiger/GENZ2014/shp/cb_2014_us_county_20m.zip", "cb_2014_us_county_20m.zip",
              method="curl", extra="-L")
unzip("cb_2014_us_county_20m.zip")
county = readOGR("cb_2014_us_county_20m","cb_2014_us_county_20m")
county$fips = paste(county$STATEFP,county$COUNTYFP,sep="")

## work out how the fips codes match up
m = match(county$fips,unem$fips)

## add the unemployment data to the county map:
county$unem = unem$percentage[m]

## fix a couple of missing data:
county$unem[is.na(county$unem)]=0

## use the purple-red palette:
colours = brewer.pal(6,"PuRd")

## make a colour scheme:
sd = data.frame(col=colours,values=c(1,3,5,7,9,11))
sc = nearestScheme(sd)

## set the plot region to the lower 48:
## the aspect ratio may not be right. Meh. This is a
## quick exercise! If I was doing this for real I'd find
## a better shapefile!

plot(c(-129,-61),c(21,53),type="n",axes=FALSE,xlab="",ylab="")

## add the counties
plot(county,col=sc(county$unem),add=TRUE,border="white",lwd=0.2)