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: http://rguide.blogspot.com/.

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: http://www.apa.org/science/psa/feb08gis.html.

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: http://geography.uoregon.edu/GeogR/index.html. 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: http://www.metrogis.org/.

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:

library(maptools)
library(rgdal)

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"))

summary(counties)

#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)))
}
names=capwords(as.vector(counties$CO_NAME))
names

library(RColorBrewer)
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
plot(water,col='#6BAED6',border=NA,add=T)
text(coordinates(counties), labels=names, cex=.75)
title(main="Twin Cities Metro Area",sub="Created with R",font.sub=2)
title(xlab="Longitude",ylab='Latitude',cex.lab=.75,line=2.25)

Leave a comment