Maps with R (I)

February 18, 2012
By

(This article was first published on Omnia sunt Communia! » R-english, and kindly contributed to R-bloggers)

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')))

http://procomun.files.wordpress.com/2012/04/spainvotes.jpeg?w=1280&h=1064


To leave a comment for the author, please follow the link and comment on his blog: Omnia sunt Communia! » R-english.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.