choropleth.R

  ##################################################################
  ## Source code for the book: "Displaying time series, spatial and
  ## space-time data with R"
  
  ## Copyright (C) 2013-2012 Oscar Perpiñán Lamigueiro
  
  ## This program is free software you can redistribute it and/or modify
  ## it under the terms of the GNU General Public License as published
  ## by the Free Software Foundation; either version 2 of the License,
  ## or (at your option) any later version.
   
  ## This program is distributed in the hope that it will be useful, but
  ## WITHOUT ANY WARRANTY; without even the implied warranty of
  ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  ## General Public License for more details.
   
  ## You should have received a copy of the GNU General Public License
  ## along with this program; if not, write to the Free Software
  ## Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  ## 02111-1307, USA.
  ####################################################################

  ##################################################################
  ## Initial configuration
  ##################################################################
  ## Clone or download the repository and set the working directory
  ## with setwd to the folder where the repository is located.
  
  library(lattice)
  library(ggplot2)
  library(latticeExtra)
  
  myTheme <- custom.theme.2(pch=19, cex=0.7,
                            region=rev(brewer.pal(9, 'YlOrRd')),
                            symbol = brewer.pal(n=8, name = "Dark2"))
  myTheme$strip.background$col='transparent'
  myTheme$strip.shingle$col='transparent'
  myTheme$strip.border$col='transparent'
  
  xscale.components.custom <- function(...){
      ans <- xscale.components.default(...)
      ans$top=FALSE
      ans}
  yscale.components.custom <- function(...){
      ans <- yscale.components.default(...)
      ans$right=FALSE
      ans}
  myArgs <- list(as.table=TRUE,
                 between=list(x=0.5, y=0.2),
                 xscale.components = xscale.components.custom,
                 yscale.components = yscale.components.custom)
  defaultArgs <- lattice.options()$default.args
  
  lattice.options(default.theme = myTheme,
                  default.args = modifyList(defaultArgs, myArgs))
  

  ##################################################################
  ## Choropleth maps
  ##################################################################

  votes2011 <- read.csv('data/votes2011.csv',
                        colClasses=c('factor', 'factor', 'numeric', 'numeric'))

  ##################################################################
  ## Administrative boundaries
  ##################################################################

  library(sp)
  library(maptools)

  old <- setwd(tempdir())
  download.file('http://goo.gl/TIvr4', 'mapas_completo_municipal.rar')
  system2('unrar', c('e', 'mapas_completo_municipal.rar'))
  espMap <- readShapePoly(fn="esp_muni_0109")
  Encoding(levels(espMap$NOMBRE)) <- "latin1"
  
  provinces <- readShapePoly(fn="spain_provinces_ag_2")
  setwd(old)

  ## dissolve repeated polygons
  espPols <- unionSpatialPolygons(espMap, espMap$PROVMUN) 

  ## Extract Canarias islands from the SpatialPolygons object
  canarias <-  sapply(espPols@polygons, function(x)substr(x@ID, 1, 2) %in% c("35",  "38"))
  peninsulaPols <- espPols[!canarias]
  islandPols <- espPols[canarias]
  
  ## Shift the island extent box to position them at the bottom right corner
  dy <- bbox(peninsulaPols)[2,1] - bbox(islandPols)[2,1]
  dx <- bbox(peninsulaPols)[1,2] - bbox(islandPols)[1,2]
  islandPols2 <- elide(islandPols, shift=c(dx, dy))
  bbIslands <- bbox(islandPols2)
  
  ## Bind Peninsula (without islands) with shifted islands
  espPols <- rbind(peninsulaPols, islandPols2)

  ## Match polygons and data using ID slot and PROVMUN column
  IDs <- sapply(espPols@polygons, function(x)x@ID)
  idx <- match(IDs, votes2011$PROVMUN)
  
  ##Places without information
  idxNA <- which(is.na(idx))
  
  ##Information to be added to the SpatialPolygons object
  dat2add <- votes2011[idx, ]
  
  ## SpatialPolygonsDataFrame uses row names to match polygons with data
  row.names(dat2add) <- IDs
  espMapVotes <- SpatialPolygonsDataFrame(espPols, dat2add)
  
  ## Drop those places without information
  espMapVotes <- espMapVotes[-idxNA, ]

  ##################################################################
  ## Map
  ##################################################################

  library(colorspace)  
  
  classes <- levels(factor(espMapVotes$whichMax))
  nClasses <- length(classes)
  
  qualPal <- rainbow_hcl(nClasses, start=30, end=300)

  ## distance between hues
  step <- 360/nClasses 
  ## hues equally spaced
  hue = (30 + step*(seq_len(nClasses)-1))%%360 
  qualPal <- hcl(hue, c=50, l=70)

pdf(file="figs/whichMax.pdf")
  spplot(espMapVotes["whichMax"], col='transparent', col.regions=qualPal)
dev.off()

pdf(file="figs/pcMax.pdf")
  quantPal <- rev(heat_hcl(16))
  spplot(espMapVotes["pcMax"], col='transparent', col.regions=quantPal)
dev.off()

##################################################################
## Categorical and quantitative variables combined in a multivariate choropleth map
##################################################################

  classes <- levels(factor(espMapVotes$whichMax))
  nClasses <- length(classes)
  step <- 360/nClasses
  multiPal <- lapply(1:nClasses, function(i){
      rev(sequential_hcl(16, h = (30 + step*(i-1))%%360))
      })

  pList <- lapply(1:nClasses, function(i){
      ## Only those polygons corresponding to a level are selected
      mapClass <- espMapVotes[espMapVotes$whichMax==classes[i],]
      pClass <- spplot(mapClass['pcMax'], col.regions=multiPal[[i]],
                       col='transparent',
                       ## labels only needed in the last legend
                       colorkey=(if (i==nClasses) TRUE else list(labels=rep('', 6))),
                       at = seq(0, 100, by=20))
  })
  
  p <- Reduce('+', pList)

  ## Function to add a title to a legend
  addTitle <- function(legend, title){
    titleGrob <- textGrob(title, gp=gpar(fontsize=8), hjust=1, vjust=1)
    ## retrieve the legend from the trellis object
    legendGrob <- eval(as.call(c(as.symbol(legend$fun), legend$args)))
    ## Layout of the legend WITH the title
    ly <- grid.layout(ncol=1, nrow=2,
                      widths=unit(0.9, 'grobwidth', data=legendGrob))
    ## Create a frame to host the original legend and the title
    fg <- frameGrob(ly, name=paste('legendTitle', title, sep='_'))
    ## Add the grobs to the frame
    pg <- packGrob(fg, titleGrob, row=2)
    pg <- packGrob(pg, legendGrob, row=1)
    }
  
  ## Access each trellis object from pList...
  for (i in seq_along(classes)){
    ## extract the legend (automatically created by spplot)...
    lg <- pList[[i]]$legend$right
    ## ... and add the addTitle function to the legend component of each trellis object
    pList[[i]]$legend$right <- list(fun='addTitle',
                                    args=list(legend=lg, title=classes[i]))
  }

  ## List of legends
  legendList <- lapply(pList, function(x){
    lg <- x$legend$right
    clKey <- eval(as.call(c(as.symbol(lg$fun), lg$args)))
    clKey
  })
  
  ## Function to pack the list of legends in a unique legend
  ## Adapted from latticeExtra::: mergedTrellisLegendGrob
  packLegend <- function(legendList){
    N <- length(legendList)
    ly <- grid.layout(nrow = 1,  ncol = N)
    g <- frameGrob(layout = ly, name = "mergedLegend")
    for (i in 1:N) g <- packGrob(g, legendList[[i]], col = i)
    g
  }
  
  ## The legend of p will include all the legends
  p$legend$right <- list(fun = 'packLegend',  args = list(legendList = legendList))
  

png(filename="figs/mapLegends.png")
  canarias <- provinces$PROV %in% c(35, 38)
  peninsulaLines <- provinces[!canarias,]
  
  p +
    layer(sp.polygons(peninsulaLines,  lwd = 0.1)) +
    layer(grid.rect(x=bbIslands[1,1], y=bbIslands[2,1],
                    width=diff(bbIslands[1,]),
                    height=diff(bbIslands[2,]),
                    default.units='native', just=c('left', 'bottom'),
                    gp=gpar(lwd=0.5, fill='transparent')))
dev.off()