Monthly Archives: September 2011

Spplot with US counties

Wow! It’s been a while! I guess I’ve been super busy, but that probably just seems like an excuse for not posting in over a year. Anyway, enough chitchatting, I wanted to create a post demonstrating how to use spplot() to plot the contiguous United States, and how to overlay US county borders. Finally, I will show how you can add data to a spatialPolygon object (counties in my case) and plot the data as a color gradient (otherwise known as a colorramp). Though my track record isn’t great, I do have a bunch of posts in mind and hope to use to this blog as resource for others (and myself) to help keep track of all my R plotting scripts.

Here goes! As always, let me know if you have any questions or comments!

First, you are going to want to download a .shp (ArcGIS shapefile) of all the US counties. There are probably a lot of different websites to get this data, but I prefer the US Census Bureau data. Make sure to save this in your working directory.

Now for some R fun:

#set working directory
setwd("C:/location_of_your_data") #Change location as this will be unique to your file directory!!!

#include the following libraries
gpclibPermit() # Required to use unionSpatialPolygons

# project for county data
aea.proj <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-102
             +x_0=180 +y_0=50 +ellps=GRS80 +datum=NAD83 +units=m"

# read in US county shapefiles
usa_counties <- readOGR("co99_d00.shp", "co99_d00")

# change projection of counties for nice plotting
proj4string(usa_counties) <- CRS("+proj=longlat +ellps=WGS84")
usa_counties <- spTransform(usa_counties, CRS(aea.proj))

# create county borders as separate feature
county_borders <- as(usa_counties, "SpatialLinesDataFrame")
county_borders <- spTransform(county_borders, CRS(aea.proj))

# create state borders by combining counties
state_borders <- unionSpatialPolygons(usa_counties,

#Plot US state borders!
#Set bounding box of plot area (I found these to work for plotting the contiguous US)
xbox = c(-2000000, 2809644)
ybox = c(-1412456,  1733192)

#NOTE: You need to use as() to convert state_borders to SpatialPolygonsDataFrame to use spplot!
spplot(as(state_borders, "SpatialPolygonsDataFrame"),
xlim=xbox, ylim=ybox, colorkey = F) #suppress plotting colorkey

Now I will show how to add data to the usa_counties SpatialPolygonsDataFrame so that you can plot county level data using a color gradient.

#First I need to create the fips code to later match the data with the correct county.
#Make FIPS code by combining usa_counties$State and usa_counties$County attributes
usa_counties$fips <- paste(usa_counties$STATE, usa_counties$COUNTY, sep="")

#Now lets select some counties.  For this example, I take the first 1000 and add our "data" (just random numbers for this example).
county_data <- data.frame("fips" = usa_counties$fips[1:1000],
"data" = rnorm(1000, 0.8, 0.3))
#match index of "usa_counties" fips with the fips from our new data.frame "county_data"
m = match(as.numeric(usa_counties$fips), county_data$fips)

# add county_data$data to usa_counties SpatialPolygonDataFrame
usa_counties$data <- county_data$data[m]
county_borders$data <- county_data$data[m]

#Create a SpatialPolygonDataFrame that is a subset of usa_counties that includes the "data"
county_subset <-  usa_counties[$data) == F,]

#Plot it up!  Note that we have new plotting arguments here. We are going to plot county_subset, but first we will use the "panel" argument to tell R to first plot the state borders (sp.polygons(...)), then we will plot the county_subset data (panel.polygonsplot(...)) on top of the state borders.  The order of the functions inside of the "panel" argument controls the order of plotting, so this is VERY IMPORTANT!
spplot(county_subset, zcol = "data", col.regions = heat.colors(100),
xlim=xbox, ylim=ybox, colorkey = T,
panel = function(x, y, ...){
sp.polygons(state_borders, zcol = "STATE", fill = "grey80")
panel.polygonsplot(x,y, ...)

UPDATE: A reader posted a comment asking how to get the borders of the county polygons to match the fill color of the polygons.  I played around with the spplot() function a bit, and found one “solution”, though probably far from the most elegant.  In fact, as it stands right now, it can match colors close, but not exactly, the same as the fill colors.  I’ve actually posted a question on R-sig-Geo seeing if anyone else has come across this problem or has a different solution.  Well here is my attempt:

Building off the example, I created a SpatialLinesDataFrame from the subset of counties that we plotted above. I then ranked them from lowest to highest, and scaled the ranking to between 1 and 100 (because this is how many heat.colors(100) we have; ie. 100 items in the vector).  Don’t forget to convert any 0s in your ranking to 1s, so that it can be used as an index!:

# Create a SpatialLinesDataFrame by converting county_subset using as()
county_borders_subset <- as(county_subset, "SpatialLinesDataFrame")

# Create scaled rank to use as index (don't forget to convert 0s to 1s!)
rank_county_data = round(rank(county_borders_subset@data$data)/
rank_county_data[rank_county_data == 0] = 1

Next, I re-plotted the map using spplot(), but added a for() loop within the panel function.  This for() loop takes each county border in “county_borders_subset” (by matching the FIPS code in the SpatialLinesDataFrame with the FIPS code vector “sub_fips”), and plots the border with the color determined by the scaled ranking of county data we added to the data frame above:

# FIPS code vector to help with subsetting county_borders_subset
sub_fips = county_borders_subset@data$fips

# Plotting
spplot(county_subset, zcol = "data", col.regions = heat.colors(100),
       xlim=xbox, ylim=ybox, colorkey = T,
       panel = function(x, y, ...){
      sp.polygons(state_borders, zcol = "STATE", fill = "grey80")
      panel.polygonsplot(x,y, ...)
      for(i in 1:length(sub_fips)){
         sub_county = subset(county_borders_subset, county_borders_subset@data$fips %in% sub_fips[i])
         sp.lines(sub_county, col = heat.colors(100)[rank_county_data[i]])

This is the final result.  It’s close, but not perfect.  One thought that did occur to me is that you can add a sp.polygons() function within the loop in addition to the sp.lines() function, and have them both use the color determined by the scaled ranking I created.  The one downside to this method, is that it’s now not completely consistent with the color ramp on the figure, but it will be very close (how close, I’m not totally sure…).  Have any other suggestions???

US counties with border colors ALMOST the same as the fill colors.

UPDATE 2:  Mathias on R-sig-Geo pointed out an obvious solution the question a commenter had:  set the “col” parameter of the spplot() function to “transparent”!  I did think about this (and came across a few other examples of people using col = “transparent”), but I was hesistant to use it for fear it would create visible gaps between the county polygons.  In fact, it does quit the opposite, and appear to only create a very, very thin gap between the polygons (only noticeable of you zoom waaaay in).  So, thanks again to Mathias.  Here is his suggestion:

# Plotting
 spplot(county_subset, zcol = "data", col.regions = heat.colors(100),
 xlim=xbox, ylim=ybox, colorkey = T, col = "transparent",
 panel = function(x, y, ...){
 sp.polygons(state_borders, zcol = "STATE", fill = "grey80")
 panel.polygonsplot(x,y, ...)

Simple, huh?  Here is the output:

US Counties using col = “transparent”