Maps with R (I)

This is the first post of a short series to show some code I have learnt to produce maps 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('votos2011.rda')
votesData
votesData$ABS
Max whichMax
## 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 

The Spanish administrative boundaries are available as shapefiles at the INE webpage (~70Mb):

espMap 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 

(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 peninsulaPols islandPols 

Then we shift the coordinates of the islands:

dy dx
islandPols2 bbIslands 

and finally construct a new object binding the shifted islands with the peninsula:

espPols 

The last step before drawing the map is to link the data with the polygons:

IDs idx
##Places without information
idxNA
##Information to be added to the SpatialPolygons object
dat2add                       poblacion = votos2011$Nombre.de.Municipio,
                      Max = Max,  pcMax = pcMax,
                      who = whichMax)[idx, ]

row.names(dat2add) espMapVotes
## Drop those places without information
espMapVotes 

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 nClasses
pList  mapClass  step  pal  pClass  at = c(0, 20, 40, 60, 80, 100))
 })

p</pre>
However, the legend of this <code>trellis</code> object is not valid.

First, a title for the legend of each element <code>pList</code> will be useful. Unfortunately, the <code>levelplot</code> function (the engine under the <code>spplot</code> method) does not allow for a title with its <code>colorkey</code> argument. The <code>frameGrob</code> and <code>packGrob</code> of the <a href="http://cran.r-project.org/web/packages/grid">grid</a> package will do the work.

1
addTitle   titleGrob   legendGrob   ly   fg   pg   pg   }

for (i in seq_along(classes)){
  lg   lg$args$key$labels$cex=ifelse(i==nClasses, 0.8, 0) ##only the last legend needs labels
  pList[[i]]$legend$right                                   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   lg   clKey   clKey
})

##function to pack the list of legends in a unique legend
##adapted from latticeExtra::: mergedTrellisLegendGrob
packLegend   N   ly   g   for (i in 1:N) g   g
}

## The legend of p will include all the legends
p$legend$right</pre>
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):

1
provinces
canarias peninsulaLines
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')))

Related articles

• Maps with R (II)
• Maps with R (III)
• Spatial data visualization with R