Previous

R and GIS

Tags:

Lately, I've been making an effort to gain experience with the R statistical language. I've motivated this by working on a side project involving election data at the precinct level, which I found at http://projects.iq.harvard.edu/eda/data.

So far I've taken the baby step of just mapping out a single state, Pennsylvania, and overlaying the counties on it. Here's my code:


library(maptools)
gpclibPermit()#to get GPCLIB to work
library(maptools)
library(gpclib)
library(maps)
library(classInt)
library(RColorBrewer)
library(sp)
require(graphics)
library(rgdal)
library(lattice)

#load data
pennStateMap <- readShapeSpatial(
"/home/jgoodson/stateData/Pennsylvania/shapeFiles/pa_final.shp") 
#get rid of Lake Erie.
pennStateMap <- subset(pennStateMap, 
                       pennStateMap@data$ALAND10/pennStateMap@data$AWATER10 > 0.01) 

#set up the palette and number of levels
plotvar <- pennStateMap@data$USPP2008
nclr <- 10
plotclr <- brewer.pal(nclr, "RdBu")
class <- classIntervals(plotvar, nclr, style="equal")
poliPallet <- colorRampPalette(c("red1", "blue1"), space="Lab") 

#merge the districts to get counties and set up the overlay
pennCounties <- unionSpatialPolygons(pennStateMap, 
                                     as.numeric(pennStateMap@data$COUNTYFP10)) 
spCounties <- list("sp.polygons", 
                   as(pennCounties, 
                   "SpatialPolygonsDataFrame"), 
                    col.regions="transparent")

#This fiddles with margins
trellis.pars.save <- trellis.par.get()  
trellis.pars <- trellis.par.get() 
layout.height <- trellis.par.get("layout.heights") 
layout.height$top.padding <- 0  
layout.height$bottom.padding <- 0 
layout.width <- trellis.par.get("layout.widths") 
layout.width$left.padding <- 2 
layout.width$right.padding <- 10  
trellis.par.set(layout.heights = layout.height) 
trellis.par.set(layout.widths = layout.width) 
trellis.par.set(axis.line=list(col=NA)) 

png(filename="test.png", width=800,height=600,units="px", 
    pointsize=12,bg="transparent",res=NA)

print(spplot(pennStateMap, "USPP2008", 
             main="2008 Presidential Election Results by Precinct in Pennsylvania", 
             col.regions=poliPallet(10), col="transparent", 
             at=round(class$brks, digits=1), 
             sp.layout = list(spCounties)))

dev.off()

This results in the following image:

Anti-Social Media

Next

Similar Posts