experimenting with maps in R

A number of resources exist for making elegant maps in R, especially from Hadley Wickham (http://had.co.nz/ggplot2/). The maps package and ggplot2 work well together to make national, state, or county borders. I need school district borders for a specific exercise. I’m no expert in mapping software (e.g. tiger files, shape files, etc.), however I know these are standard file formats for most mapping datasets. The Census website (http://www.census.gov/cgi-bin/geo/shapefiles2010/main) houses publicly available shapefiles for school districts. So, I used ggplot2 to play around with with these files. Here is the code and output:

library(maptools)
library(gpclib)
library(ggplot2)
d1 <- readShapeSpatial("tl_2010_39_unsd10.shp")
plot(d1)
meta <- as.data.frame(d1)
gpclibPermit()
d1ddf <- fortify(d1, region="NAME10")
id <- d1ddf$id
id <- id[!duplicated(id)]
y <- rnorm(length(id))
perf <- data.frame(id, y)
d2ddf <- merge(d1ddf, perf, by="id")
q1 <- quantile(d2ddf$y, c(.2, .4, .6, .8), na.rm=T)[1]
q2 <- quantile(d2ddf$y, c(.2, .4, .6, .8), na.rm=T)[2]
q3 <- quantile(d2ddf$y, c(.2, .4, .6, .8), na.rm=T)[3]
q4 <- quantile(d2ddf$y, c(.2, .4, .6, .8), na.rm=T)[4]
d2ddf$y_cat[d2ddf$y < q1] <- "20th Percentile or Below"
d2ddf$y_cat[d2ddf$y >= q1 & d2ddf$y<=q2] <- "21st-40th Percentile"
d2ddf$y_cat[d2ddf$y > q2 & d2ddf$y<=q3] <- "41st-60th Percentile"
d2ddf$y_cat[d2ddf$y > q3 & d2ddf$y<=q4] <- "61st-80th Percentile"
d2ddf$y_cat[d2ddf$y > q4] <- "Above 80th Percentile"
d2ddf$y_cat <- ordered(factor(d2ddf$y_cat, levels=c("20th Percentile or Below",
                                                    "21st-40th Percentile",
                                                    "41st-60th Percentile",
                                                    "61st-80th Percentile",
                                                    "Above 80th Percentile")))
p <- ggplot(d2ddf)
p1 <- p + geom_polygon(aes(long, lat, group=group))
p2 <- p1 + geom_polygon(aes(fill=y_cat,long,lat,group=group)) +
           geom_polygon(data = d1, colour = alpha("white", 1/2),
           size = 0.2, fill = NA) +
           scale_fill_brewer(pal="PuRd", name="Outcome") +
           geom_path(aes(long,lat,group=group),colour="white") +
           opts(title="OH Outcomes") +
           xlab("Longitude") +
           ylab("Latitude")
p2

This is probably not the most efficient way to produce maps like this (especially since this outcome is arbitrary) but I’m pleased with the output:

This entry was posted in Uncategorized. Bookmark the permalink.

Leave a Reply

Your email address will not be published. Required fields are marked *

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>