This is the first post of a short series to show some code I have learnt to produce maps with R.
NOTE: Although the procedure described in this post is valid, there is a newer code version in one of the chapters of the book «Displaying time series, spatial and space-time data with R«
Some time ago I found this infographic from The New York Times (via this page) and I wondered how a multivariate choropleth map could be produced with R. Here is the code I have arranged to show the results of the last Spanish general elections in a similar fashion.
Some packages are needed:
library(maps)
library(maptools)
## EDITED: if you have rgeos installed you won't need
## gpclibPermit() below.
library(sp)
library(lattice)
library(latticeExtra)
library(colorspace)
Let’s start with the data, which is available here (thanks to Emilio Torres, who «massaged» the original dataset, available from here).
Each region of the map will represent the percentage of votes obtained by the predominant political option. Besides, only four groups will be considered: the two main parties («PP» and «PSOE»), the abstention results («ABS»), and the rest of parties («OTH»).
load(url('http://dl.dropbox.com/u/40293713/spainVotes/votos2011.rda'))
votesData <- votos2011[, 12:1023] ##I don't need all the columns
votesData$ABS <- with(votos2011, Total.censo.electoral - Votos.validos) ##abstention
Max <- apply(votesData, 1, max)
whichMax <- apply(votesData, 1, function(x)names(votesData)[which.max(x)])
## OTH for everything but PP, PSOE and ABS
whichMax[!(whichMax %in% c('PP', 'PSOE', 'ABS'))] <- 'OTH'
## Finally, I calculate the percentage of votes with the electoral census
pcMax <- Max/votos2011$Total.censo.electoral * 100
The Spanish administrative boundaries are available as shapefiles at the INE webpage (~70Mb):
espMap <- readShapePoly(fn="mapas_completo_municipal/esp_muni_0109")
Encoding(levels(espMap$NOMBRE)) <- "latin1"
##There are some repeated polygons which can be dissolved with
## unionSpatialPolygons.
## EDITED: gpclibPermit() is needed for unionSpatialPolygons to work
## but can be ommited if you have rgeos installed
## (recommended, see comment of Roger Bivand below).
espPols <- unionSpatialPolygons(espMap, espMap$PROVMUN)
(EDITED, following the question of Sandra). Spanish maps are commonly displayed with the Canarian islands next to the peninsula. First we have to extract the polygons of the islands and the polygons of the peninsula.
canarias <- sapply(espPols@polygons, function(x)substr(x@ID, 1, 2) %in% c("35", "38"))
peninsulaPols <- espPols[!canarias]
islandPols <- espPols[canarias]
Then we shift the coordinates of the islands:
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)
and finally construct a new object binding the shifted islands with the peninsula:
espPols <- rbind(peninsulaPols, islandPols2)
The last step before drawing the map is to link the data with the polygons:
IDs <- sapply(espPols@polygons, function(x)x@ID)
idx <- match(IDs, votos2011$PROVMUN)
##Places without information
idxNA <- which(is.na(idx))
##Information to be added to the SpatialPolygons object
dat2add <- data.frame(prov = votos2011$PROV,
poblacion = votos2011$Nombre.de.Municipio,
Max = Max, pcMax = pcMax, who = whichMax)[idx, ]
row.names(dat2add) <- IDs
espMapVotes <- SpatialPolygonsDataFrame(espPols, dat2add)
## Drop those places without information
espMapVotes <- espMapVotes[-idxNA, ]
So let’s draw the map. I will produce a list of plots, one for each group. The «+.trellis» method of the latticeExtra package with Reduce superposes the elements of this list and produce a trellis object. I will use a set of sequential palettes from the colorspace package with a different hue for each group.
classes <- levels(factor(whichMax))
nClasses <- length(classes)
pList <- lapply(1:nClasses, function(i){
mapClass <- espMapVotes[espMapVotes$who==classes[i],]
step <- 360/nClasses ## distance between hues
pal <- rev(sequential_hcl(16, h = (30 + step*(i-1))%%360)) ## hues equally spaced
pClass <- spplot(mapClass['pcMax'], col.regions=pal, lwd=0.1,
at = c(0, 20, 40, 60, 80, 100))
})
p <- Reduce('+', pList)
However, the legend of this trellis object is not valid.
First, a title for the legend of each element pList will be
useful. Unfortunately, the levelplot function (the engine under the
spplot method) does not allow for a title with its colorkey
argument. The frameGrob and packGrob of the grid package will do the work.
addTitle <- function(legend, title){
titleGrob <- textGrob(title, gp=gpar(fontsize=8), hjust=1, vjust=1)
legendGrob <- eval(as.call(c(as.symbol(legend$fun), legend$args)))
ly <- grid.layout(ncol=1, nrow=2, widths=unit(0.9, 'grobwidth', data=legendGrob))
fg <- frameGrob(ly, name=paste('legendTitle', title, sep='_'))
pg <- packGrob(fg, titleGrob, row=2)
pg <- packGrob(pg, legendGrob, row=1)
}
for (i in seq_along(classes)){
lg <- pList[[i]]$legend$right
lg$args$key$labels$cex=ifelse(i==nClasses, 0.8, 0) ##only the last legend needs labels
pList[[i]]$legend$right <- list(fun='addTitle',
args=list(legend=lg, title=classes[i]))
}
Now, every component of pList includes a legend with a title below. The last step is to modify the legend of the p trellis object in order to merge the legends from every component of pList.
## 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))
Here is the result with the provinces boundaries superposed (only for the peninsula due to a problem with the definition of boundaries the Canarian islands in the file) and a rectangle to separate the Canarian islands from the
rest of the map (click on the image to get a SVG file):
provinces <- readShapePoly(fn="mapas_completo_municipal/spain_provinces_ag_2")
canarias <- provinces$PROV %in% c(35, 38)
peninsulaLines <- provinces[!canarias,]
p +
layer(sp.polygons(peninsulaLines, lwd = 0.5)) +
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')))
Me gusta esto:
Me gusta Cargando...
Omnia sunt Communia! 
