Mapping with R

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 cap 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 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)
This entry was posted in Praxes. Bookmark the permalink.