Mapping with R

| No Comments

I have been learning R software for statistical computing. I have used it for all of my homework assignments in EPSY 8261 (Probability and Inference) and 8262 (Regression and the General Linear Model). My prof has an R blog here:

I envision bridging my studies in educational psychology (quantitative methods, evaluation) with my interest in spatial analysis and mapping. I read yesterday that APA's science directorate is trying to encourage greater use of GIS in psychological research. You can read more here:

I decided to explore R's mapping capabilities today. Here are a couple findings:
-The shapefiles package is not really necessary.
-The maptools package can read and display shapefiles, but the author encourages converting shapefiles into sp class or polylist objects.
-Once shapefiles are converted, one can use the maptools plot.polylist polygon helper function.

Pat Bartlein at the University of Oregon recommends the following packages: sp, maptools, rgdal, maps, mapproj, mapdata, classInt, scatterplot3d, and RColorBrewer. He created a great website for folks like me who are starting out using R for geography: Thanks, Pat!

Applied researchers in the Twin Cities unfamiliar with GIS should know that the Metropolitan Council and its MetroGIS initiative have done a great job of collaborating with local entities to make a ton of data available at no charge:

I prefer ArcGIS for making maps and performing spatial analyses. I use MapInfo at work because our IT/geographer person prefers it. Mapping with R is similar to statistical computing in R--it's a blank portrait waiting for color and texture, whereas SPSS and commercial GIS software are like the Etch A Sketch.

Click on the thumbnail for a large version of the Twin Cities map I created with R:

Twin Cities Metro Area 041208.png

Here is the script:


counties=readShapePoly(file.choose(),proj4string=CRS("+proj=utm +zone=15 +ellps=GRS80 +datum=NAD83")) #choose MetroGIS' county.shp
counties=spTransform(counties, CRS=CRS("+proj=longlat"))


#capitalize first letter of county names
capwords <- function(s) {
cap <- function(s) paste( toupper(substring(s,1,1)),
tolower(substring(s,2)), sep = "", collapse = " " )
sapply(strsplit(s, split = " "), cap, USE.NAMES = !is.null(names(s)))

display.brewer.all() #check out the color options
cols <- c(brewer.pal(7,"Pastel2"))

water=readShapePoly(file.choose(),proj4string=CRS("+proj=utm +zone=15 +ellps=GRS80 +datum=NAD83")) #choose MetroGIS' hydro_luse2000.shp
water=spTransform(water, CRS=CRS("+proj=longlat"))

coordinates(counties) #centroids

plot(counties,col=cols,axes=T,cex.axis=.75) #plot function in maptools package
text(coordinates(counties), labels=names, cex=.75)
title(main="Twin Cities Metro Area",sub="Created with R",font.sub=2)

Leave a comment