Search code examples
rdictionarygeolocationgeospatial

Developing Geographic Thematic Maps with R

There are clearly a number of packages in R for all sorts of spatial analysis. That can by seen in the CRAN Task View: Analysis of Spatial Data. These packages are numerous and diverse, but all I want to do is some simple thematic maps. I have data with county and state FIPS codes and I have ESRI shape files of county and state boundaries and the accompanying FIPS codes which allows joining with the data. The shape files could be easily converted to other formats, if needed.

So what's the most straight forward way to create thematic maps with R?

This map looks like it was created with an ESRI Arc product, but this is the type of thing I would like to do with R:

alt text http://www.infousagov.com/images/choro.jpg Map copied from here.

Solution

  • The following code has served me well. Customize it a little and you are done. alt text
    (source: eduardoleoni.com) 

    library(maptools)
substitute your shapefiles here
state.map <- readShapeSpatial("BRASIL.shp")
counties.map <- readShapeSpatial("55mu2500gsd.shp")
## this is the variable we will be plotting
counties.map@data$noise <- rnorm(nrow(counties.map@data))

    heatmap function

    plot.heat <- function(counties.map,state.map,z,title=NULL,breaks=NULL,reverse=FALSE,cex.legend=1,bw=.2,col.vec=NULL,plot.legend=TRUE) {
  ##Break down the value variable
  if (is.null(breaks)) {
    breaks=
      seq(
          floor(min(counties.map@data[,z],na.rm=TRUE)*10)/10
          ,
          ceiling(max(counties.map@data[,z],na.rm=TRUE)*10)/10
          ,.1)
  }
  counties.map@data$zCat <- cut(counties.map@data[,z],breaks,include.lowest=TRUE)
  cutpoints <- levels(counties.map@data$zCat)
  if (is.null(col.vec)) col.vec <- heat.colors(length(levels(counties.map@data$zCat)))
  if (reverse) {
    cutpointsColors <- rev(col.vec)
  } else {
    cutpointsColors <- col.vec
  }
  levels(counties.map@data$zCat) <- cutpointsColors
  plot(counties.map,border=gray(.8), lwd=bw,axes = FALSE, las = 1,col=as.character(counties.map@data$zCat))
  if (!is.null(state.map)) {
    plot(state.map,add=TRUE,lwd=1)
  }
  ##with(counties.map.c,text(x,y,name,cex=0.75))
  if (plot.legend) legend("bottomleft", cutpoints, fill = cutpointsColors,bty="n",title=title,cex=cex.legend)
  ##title("Cartogram")
}

    plot it

    plot.heat(counties.map,state.map,z="noise",breaks=c(-Inf,-2,-1,0,1,2,Inf))