{"id":845,"date":"2015-04-01T14:02:27","date_gmt":"2015-04-01T21:02:27","guid":{"rendered":"http:\/\/homepages.uc.edu\/~yaozo\/wordpress\/?p=845"},"modified":"2015-04-01T14:02:27","modified_gmt":"2015-04-01T21:02:27","slug":"maps-with-r","status":"publish","type":"post","link":"https:\/\/zhuoyao.net\/index.php\/2015\/04\/01\/maps-with-r\/","title":{"rendered":"Maps with R"},"content":{"rendered":"

This is the first post of a short series to show some code I have learnt to produce maps with R.<\/p>\n

NOTE: Although the procedure described in this post is valid, there is anewer code version<\/a> in one of the chapters of the book \u201cDisplaying time series, spatial and space-time data with R<\/a>\u201c<\/strong><\/p>\n

Some time ago I found this infographic from The New York Times<\/a> (via this page<\/a>) 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<\/a> in a similar fashion.<\/p>\n

Some packages are needed:<\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n
3<\/div>\n
4<\/div>\n
5<\/div>\n
6<\/div>\n
7<\/div>\n
8<\/div>\n<\/td>\n
\n
\n
library<\/code>(maps)<\/code><\/div>\n
library<\/code>(maptools)<\/code><\/div>\n
## EDITED: if you have rgeos installed you won't need<\/code><\/div>\n
## gpclibPermit() below.<\/code><\/div>\n
library<\/code>(sp)<\/code><\/div>\n
library<\/code>(lattice)<\/code><\/div>\n
library<\/code>(latticeExtra)<\/code><\/div>\n
library<\/code>(colorspace)<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

Let\u2019s start with the data, which is available here<\/a> (thanks to Emilio Torres<\/a>, who \u201cmassaged\u201d the original dataset, available from here<\/a>).<\/p>\n

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 (\u201cPP\u201d and \u201cPSOE\u201d), the abstention results (\u201cABS\u201d), and the rest of parties (\u201cOTH\u201d).<\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n
3<\/div>\n
4<\/div>\n
5<\/div>\n
6<\/div>\n
7<\/div>\n
8<\/div>\n
9<\/div>\n
10<\/div>\n
11<\/div>\n
12<\/div>\n
13<\/div>\n<\/td>\n
\n
\n
load<\/code>(<\/code>url<\/code>(<\/code>'http:\/\/dl.dropbox.com\/u\/40293713\/spainVotes\/votos2011.rda<\/a>'<\/code>))<\/code><\/div>\n
votesData <- votos2011[, 12:1023] <\/code>##I don't need all the columns<\/code><\/div>\n
<\/div>\n
votesData$ABS <- <\/code>with<\/code>(votos2011, Total.censo.electoral - Votos.validos) <\/code>##abstention<\/code><\/div>\n
<\/div>\n
Max <- <\/code>apply<\/code>(votesData, 1, max)<\/code><\/div>\n
whichMax <- <\/code>apply<\/code>(votesData,\u00a0 1, <\/code>function<\/code>(x)<\/code>names<\/code>(votesData)[<\/code>which.max<\/code>(x)])<\/code><\/div>\n
<\/div>\n
## OTH for everything but PP, PSOE and ABS<\/code><\/div>\n
whichMax[!(whichMax %<\/code>in<\/code>% <\/code>c<\/code>(<\/code>'PP'<\/code>,\u00a0 <\/code>'PSOE'<\/code>, <\/code>'ABS'<\/code>))] <- <\/code>'OTH'<\/code><\/div>\n
<\/div>\n
## Finally, I calculate the percentage of votes with the electoral census<\/code><\/div>\n
pcMax <- Max\/votos2011$Total.censo.electoral * 100<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

The Spanish administrative boundaries are available as shapefiles at the INE webpage<\/a> (~70Mb):<\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n
3<\/div>\n
4<\/div>\n
5<\/div>\n
6<\/div>\n
7<\/div>\n
8<\/div>\n
9<\/div>\n<\/td>\n
\n
\n
espMap <- <\/code>readShapePoly<\/code>(fn=<\/code>\"mapas_completo_municipal\/esp_muni_0109\"<\/code>)<\/code><\/div>\n
Encoding<\/code>(<\/code>levels<\/code>(espMap$NOMBRE)) <- <\/code>\"latin1\"<\/code><\/div>\n
<\/div>\n
##There are some repeated polygons which can be dissolved with<\/code><\/div>\n
## unionSpatialPolygons.<\/code><\/div>\n
## EDITED: gpclibPermit() is needed for unionSpatialPolygons to work<\/code><\/div>\n
## but can be ommited if you have rgeos installed<\/code><\/div>\n
## (recommended, see comment of Roger Bivand below).<\/code><\/div>\n
espPols <- <\/code>unionSpatialPolygons<\/code>(espMap, espMap$PROVMUN)<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

(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.<\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n
3<\/div>\n<\/td>\n
\n
\n
canarias <-\u00a0 <\/code>sapply<\/code>(espPols@polygons, <\/code>function<\/code>(x)<\/code>substr<\/code>(x@ID, 1, 2) %<\/code>in<\/code>% <\/code>c<\/code>(<\/code>\"35\"<\/code>,\u00a0 <\/code>\"38\"<\/code>))<\/code><\/div>\n
peninsulaPols <- espPols[!canarias]<\/code><\/div>\n
islandPols <- espPols[canarias]<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

Then we shift the coordinates of the islands:<\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n
3<\/div>\n
4<\/div>\n
5<\/div>\n<\/td>\n
\n
\n
dy <- <\/code>bbox<\/code>(peninsulaPols)[2,1] - <\/code>bbox<\/code>(islandPols)[2,1]<\/code><\/div>\n
dx <- <\/code>bbox<\/code>(peninsulaPols)[1,2] - <\/code>bbox<\/code>(islandPols)[1,2]<\/code><\/div>\n
<\/div>\n
islandPols2 <- <\/code>elide<\/code>(islandPols, shift=<\/code>c<\/code>(dx, dy))<\/code><\/div>\n
bbIslands <- <\/code>bbox<\/code>(islandPols2)<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

and finally construct a new object binding the shifted islands with the peninsula:<\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n<\/td>\n
\n
\n
espPols <- <\/code>rbind<\/code>(peninsulaPols, islandPols2)<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

The last step before drawing the map is to link the data with the polygons:<\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n
3<\/div>\n
4<\/div>\n
5<\/div>\n
6<\/div>\n
7<\/div>\n
8<\/div>\n
9<\/div>\n
10<\/div>\n
11<\/div>\n
12<\/div>\n
13<\/div>\n
14<\/div>\n
15<\/div>\n
16<\/div>\n<\/td>\n
\n
\n
IDs <- <\/code>sapply<\/code>(espPols@polygons, <\/code>function<\/code>(x)x@ID)<\/code><\/div>\n
idx <- <\/code>match<\/code>(IDs, votos2011$PROVMUN)<\/code><\/div>\n
<\/div>\n
##Places without information<\/code><\/div>\n
idxNA <- <\/code>which<\/code>(<\/code>is.na<\/code>(idx))<\/code><\/div>\n
<\/div>\n
##Information to be added to the SpatialPolygons object<\/code><\/div>\n
dat2add <- <\/code>data.frame<\/code>(prov = votos2011$PROV,<\/code><\/div>\n
poblacion = votos2011$Nombre.de.Municipio,<\/code><\/div>\n
Max = Max,\u00a0 pcMax = pcMax,\u00a0 who = whichMax)[idx, ]<\/code><\/div>\n
<\/div>\n
row.names<\/code>(dat2add) <- IDs<\/code><\/div>\n
espMapVotes <- <\/code>SpatialPolygonsDataFrame<\/code>(espPols, dat2add)<\/code><\/div>\n
<\/div>\n
## Drop those places without information<\/code><\/div>\n
espMapVotes <- espMapVotes[-idxNA, ]<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

So let\u2019s draw the map. I will produce a list of plots, one for each group. The\u201c+.trellis\u201d method<\/a> of the latticeExtra<\/a> package with Reduce<\/code> superposes the elements of this list and produce a trellis<\/code> object. I will use a set of sequential palettes from the colorspace<\/a> package with a different hue<\/a> for each group.<\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n
3<\/div>\n
4<\/div>\n
5<\/div>\n
6<\/div>\n
7<\/div>\n
8<\/div>\n
9<\/div>\n
10<\/div>\n
11<\/div>\n
12<\/div>\n<\/td>\n
\n
\n
classes <- <\/code>levels<\/code>(<\/code>factor<\/code>(whichMax))<\/code><\/div>\n
nClasses <- <\/code>length<\/code>(classes)<\/code><\/div>\n
<\/div>\n
pList <- <\/code>lapply<\/code>(1:nClasses, <\/code>function<\/code>(i){<\/code><\/div>\n
\u00a0\u00a0<\/code>mapClass <- espMapVotes[espMapVotes$who==classes[i],]<\/code><\/div>\n
\u00a0\u00a0<\/code>step <- 360\/nClasses <\/code>## distance between hues<\/code><\/div>\n
\u00a0\u00a0<\/code>pal <- <\/code>rev<\/code>(<\/code>sequential_hcl<\/code>(16, h = (30 + step*(i-1))%%360)) <\/code>## hues equally spaced<\/code><\/div>\n
\u00a0\u00a0<\/code>pClass <- <\/code>spplot<\/code>(mapClass[<\/code>'pcMax'<\/code>], col.regions=pal, lwd=0.1,<\/code><\/div>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0<\/code>at = <\/code>c<\/code>(0, 20, 40, 60, 80, 100))<\/code><\/div>\n
\u00a0\u00a0<\/code>})<\/code><\/div>\n
<\/div>\n
p <- <\/code>Reduce<\/code>(<\/code>'+'<\/code>, pList)<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

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

First, a title for the legend of each element pList<\/code> will be
\nuseful. Unfortunately, the levelplot function (the engine under the
\nspplot method) does not allow for a title with its colorkey<\/code>
\nargument. The frameGrob<\/code> and packGrob<\/code> of the grid<\/code> package will do the work.<\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n
3<\/div>\n
4<\/div>\n
5<\/div>\n
6<\/div>\n
7<\/div>\n
8<\/div>\n
9<\/div>\n
10<\/div>\n
11<\/div>\n
12<\/div>\n
13<\/div>\n
14<\/div>\n
15<\/div>\n<\/td>\n
\n
\n
addTitle <- <\/code>function<\/code>(legend, title){<\/code><\/div>\n
\u00a0\u00a0<\/code>titleGrob <- <\/code>textGrob<\/code>(title, gp=<\/code>gpar<\/code>(fontsize=8), hjust=1, vjust=1)<\/code><\/div>\n
\u00a0\u00a0<\/code>legendGrob <- <\/code>eval<\/code>(<\/code>as.call<\/code>(<\/code>c<\/code>(<\/code>as.symbol<\/code>(legend$fun), legend$args)))<\/code><\/div>\n
\u00a0\u00a0<\/code>ly <- <\/code>grid.layout<\/code>(ncol=1, nrow=2, widths=<\/code>unit<\/code>(0.9, <\/code>'grobwidth'<\/code>, data=legendGrob))<\/code><\/div>\n
\u00a0\u00a0<\/code>fg <- <\/code>frameGrob<\/code>(ly, name=<\/code>paste<\/code>(<\/code>'legendTitle'<\/code>, title, sep=<\/code>'_'<\/code>))<\/code><\/div>\n
\u00a0\u00a0<\/code>pg <- <\/code>packGrob<\/code>(fg, titleGrob, row=2)<\/code><\/div>\n
\u00a0\u00a0<\/code>pg <- <\/code>packGrob<\/code>(pg, legendGrob, row=1)<\/code><\/div>\n
\u00a0\u00a0<\/code>}<\/code><\/div>\n
<\/div>\n
for <\/code>(i <\/code>in<\/code> seq_along<\/code>(classes)){<\/code><\/div>\n
\u00a0\u00a0<\/code>lg <- pList[[i]]$legend$right<\/code><\/div>\n
\u00a0\u00a0<\/code>lg$args$key$labels$cex=<\/code>ifelse<\/code>(i==nClasses, 0.8, 0) <\/code>##only the last legend needs labels<\/code><\/div>\n
\u00a0\u00a0<\/code>pList[[i]]$legend$right <- <\/code>list<\/code>(fun=<\/code>'addTitle'<\/code>,<\/code><\/div>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0<\/code>args=<\/code>list<\/code>(legend=lg, title=classes[i]))<\/code><\/div>\n
}<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

Now, every component of pList<\/code> includes a legend with a title below. The last step is to modify the legend of the p<\/code> trellis object in order to merge the legends from every component of pList<\/code>.<\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n
3<\/div>\n
4<\/div>\n
5<\/div>\n
6<\/div>\n
7<\/div>\n
8<\/div>\n
9<\/div>\n
10<\/div>\n
11<\/div>\n
12<\/div>\n
13<\/div>\n
14<\/div>\n
15<\/div>\n
16<\/div>\n
17<\/div>\n
18<\/div>\n
19<\/div>\n<\/td>\n
\n
\n
## list of legends<\/code><\/div>\n
legendList <- <\/code>lapply<\/code>(pList, <\/code>function<\/code>(x){<\/code><\/div>\n
\u00a0\u00a0<\/code>lg <- x$legend$right<\/code><\/div>\n
\u00a0\u00a0<\/code>clKey <- <\/code>eval<\/code>(<\/code>as.call<\/code>(<\/code>c<\/code>(<\/code>as.symbol<\/code>(lg$fun), lg$args)))<\/code><\/div>\n
\u00a0\u00a0<\/code>clKey<\/code><\/div>\n
})<\/code><\/div>\n
<\/div>\n
##function to pack the list of legends in a unique legend<\/code><\/div>\n
##adapted from latticeExtra::: mergedTrellisLegendGrob<\/code><\/div>\n
packLegend <- <\/code>function<\/code>(legendList){<\/code><\/div>\n
\u00a0\u00a0<\/code>N <- <\/code>length<\/code>(legendList)<\/code><\/div>\n
\u00a0\u00a0<\/code>ly <- <\/code>grid.layout<\/code>(nrow = 1,\u00a0 ncol = N)<\/code><\/div>\n
\u00a0\u00a0<\/code>g <- <\/code>frameGrob<\/code>(layout = ly, name = <\/code>\"mergedLegend\"<\/code>)<\/code><\/div>\n
\u00a0\u00a0<\/code>for <\/code>(i <\/code>in<\/code> 1:N) g <- <\/code>packGrob<\/code>(g, legendList[[i]], col = i)<\/code><\/div>\n
\u00a0\u00a0<\/code>g<\/code><\/div>\n
}<\/code><\/div>\n
<\/div>\n
## The legend of p will include all the legends<\/code><\/div>\n
p$legend$right <- <\/code>list<\/code>(fun = <\/code>'packLegend'<\/code>,\u00a0 args = <\/code>list<\/code>(legendList = legendList))<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

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
\nrest of the map (click on the image to get a SVG file):<\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n
3<\/div>\n
4<\/div>\n
5<\/div>\n
6<\/div>\n
7<\/div>\n
8<\/div>\n
9<\/div>\n
10<\/div>\n
11<\/div>\n
12<\/div>\n<\/td>\n
\n
\n
provinces <- <\/code>readShapePoly<\/code>(fn=<\/code>\"mapas_completo_municipal\/spain_provinces_ag_2\"<\/code>)<\/code><\/div>\n
<\/div>\n
canarias <- provinces$PROV %<\/code>in<\/code>% <\/code>c<\/code>(35, 38)<\/code><\/div>\n
peninsulaLines <- provinces[!canarias,]<\/code><\/div>\n
<\/div>\n
p +<\/code><\/div>\n
\u00a0\u00a0<\/code>layer<\/code>(<\/code>sp.polygons<\/code>(peninsulaLines,\u00a0 lwd = 0.5)) +<\/code><\/div>\n
\u00a0\u00a0<\/code>layer<\/code>(<\/code>grid.rect<\/code>(x=bbIslands[1,1], y=bbIslands[2,1],<\/code><\/div>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0<\/code>width=<\/code>diff<\/code>(bbIslands[1,]),<\/code><\/div>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0<\/code>height=<\/code>diff<\/code>(bbIslands[2,]),<\/code><\/div>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0<\/code>default.units=<\/code>'native'<\/code>, just=<\/code>c<\/code>(<\/code>'left'<\/code>, <\/code>'bottom'<\/code>),<\/code><\/div>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0<\/code>gp=<\/code>gpar<\/code>(lwd=0.5, fill=<\/code>'transparent'<\/code>)))<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

\"https:\/\/procomun.files.wordpress.com\/2012\/04\/spainvotes.jpeg\"<\/a><\/p>\n

In my my last post<\/a> I described how to produce a multivariate choropleth map with R. Now I will show how to create a map from raster files. One of them is a factor which will group the values of the other one. Thus, once again, I will superpose several groups in the same map.<\/p>\n

NOTE: Although the procedure described in this post is valid, there is anewer code version<\/a> in one of the chapters of the book \u201cDisplaying time series, spatial and space-time data with R<\/a>\u201c<\/strong><\/p>\n

First let\u2019s load the packages.<\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n
3<\/div>\n<\/td>\n
\n
\n
library<\/code>(raster)<\/code><\/div>\n
library<\/code>(rasterVis)<\/code><\/div>\n
library<\/code>(colorspace)<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

Now, I define the geographical extent to be analyzed (approximately India and China).<\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n<\/td>\n
\n
\n
ext <- <\/code>extent<\/code>(65, 135, 5, 55)<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

The first raster file is the population density<\/a> in our planet, available at thisNEO-NASA webpage<\/a> (choose the Geo-TIFF<\/a> floating option, ~25Mb). After reading the data with raster<\/code> I subset the geographical extent and replace the 99999 with NA<\/code>.<\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n
3<\/div>\n
4<\/div>\n
5<\/div>\n<\/td>\n
\n
\n
pop <- <\/code>raster<\/code>(<\/code>'875430rgb-167772161.0.FLOAT.TIFF'<\/code>)<\/code><\/div>\n
pop <- <\/code>crop<\/code>(pop, ext)<\/code><\/div>\n
pop[pop==99999] <- <\/code>NA<\/code><\/div>\n
pTotal <- <\/code>levelplot<\/code>(pop, zscaleLog=10, par.settings=BTCTheme)<\/code><\/div>\n
pTotal<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

\"\"<\/a>The second raster file is the land cover classification (available at this NEO-NASA webpage<\/a>)<\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n<\/td>\n
\n
\n
landClass <- <\/code>raster<\/code>(<\/code>'241243rgb-167772161.0.TIFF'<\/code>)<\/code><\/div>\n
landClass <- <\/code>crop<\/code>(landClass, ext)<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

The codes of the classification are described here<\/a>. In summary, the sea is labeled with 0, forests with 1 to 5, shrublands, grasslands and wetlands with 6 to 11, agriculture and urban lands with 12 to 14, and snow and barren with 15 and 16. This four groups (sea is replaced NA<\/code>) will be the levels of the factor. (I am not sure if these sets of different land covers is sensible: comments from experts are welcome!)<\/p>\n

EDIT: Following a question from a user of rasterVis<\/a> I include some lines of code to display this qualitative variable in the map.
\nEDIT2: raster<\/code> and rasterVis<\/code> are able to work with categorical data using ratify.
\n<\/code><\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n
3<\/div>\n
4<\/div>\n
5<\/div>\n
6<\/div>\n
7<\/div>\n
8<\/div>\n
9<\/div>\n
10<\/div>\n
11<\/div>\n
12<\/div>\n
13<\/div>\n<\/td>\n
\n
\n
landClass[landClass %<\/code>in<\/code>% <\/code>c<\/code>(0, 254)] <- <\/code>NA<\/code><\/div>\n
landClass <- <\/code>cut<\/code>(landClass, <\/code>c<\/code>(0, 5, 11, 14, 16))<\/code><\/div>\n
## Add a Raster Atribute Table and define the raster as categorical data<\/code><\/div>\n
landClass <- <\/code>ratify<\/code>(landClass)<\/code><\/div>\n
## Configure the RAT: first create a RAT data.frame using the<\/code><\/div>\n
## levels method; second, set the values for each class (to be<\/code><\/div>\n
## used by levelplot); third, assign this RAT to the raster<\/code><\/div>\n
## using again levels<\/code><\/div>\n
rat <- <\/code>levels<\/code>(landClass)[[1]]<\/code><\/div>\n
rat$classes <- <\/code>c<\/code>(<\/code>'Forest'<\/code>, <\/code>'Land'<\/code>, <\/code>'Urban'<\/code>, <\/code>'Snow'<\/code>)<\/code><\/div>\n
levels<\/code>(landClass) <- rat<\/code><\/div>\n
<\/div>\n
levelplot<\/code>(landClass, col.regions=<\/code>terrain_hcl<\/code>(4))<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

\"\"<\/a><\/p>\n

This histogram shows the distribution of the population density in each land class.<\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n
3<\/div>\n
4<\/div>\n<\/td>\n
\n
\n
s <- <\/code>stack<\/code>(pop, landClass)<\/code><\/div>\n
layerNames<\/code>(s) <- <\/code>c<\/code>(<\/code>'pop'<\/code>, <\/code>'landClass'<\/code>)<\/code><\/div>\n
histogram<\/code>(~<\/code>log10<\/code>(pop)|landClass, data=s,<\/code><\/div>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0<\/code>scales=<\/code>list<\/code>(relation=<\/code>'free'<\/code>))<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

\"\"<\/a>Everything is ready for the map. I will create a list of trellis<\/code> objects with four elements (one for each level of the factor). Each of these objects is the representation of the population density in a particular land class. I use the same scale for all of them to allow for comparisons (the at<\/code> argument of levelplot receives the correspondent at<\/code> values from the global<\/em> map)<\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n
3<\/div>\n
4<\/div>\n
5<\/div>\n
6<\/div>\n
7<\/div>\n
8<\/div>\n
9<\/div>\n
10<\/div>\n
11<\/div>\n<\/td>\n
\n
\n
at <- pTotal$legend$bottom$args$key$at<\/code><\/div>\n
<\/div>\n
pList <- <\/code>lapply<\/code>(1:nClasses, <\/code>function<\/code>(i){<\/code><\/div>\n
\u00a0\u00a0<\/code>landSub <- landClass<\/code><\/div>\n
\u00a0\u00a0<\/code>landSub[!(landClass==i)] <- <\/code>NA<\/code><\/div>\n
\u00a0\u00a0<\/code>popSub <- <\/code>mask<\/code>(pop, landSub)<\/code><\/div>\n
\u00a0\u00a0<\/code>step <- 360\/nClasses<\/code><\/div>\n
\u00a0\u00a0<\/code>pal <- <\/code>rev<\/code>(<\/code>sequential_hcl<\/code>(16, h = (30 + step*(i-1))%%360))<\/code><\/div>\n
\u00a0\u00a0<\/code>pClass <- <\/code>levelplot<\/code>(popSub, zscaleLog=10, at=at,<\/code><\/div>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0<\/code>col.regions=pal, margin=<\/code>FALSE<\/code>)<\/code><\/div>\n
})<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

And that\u2019s all. The rest of the code is exactly the same as in the previous post<\/a>. If you execute it you will get this image (click on it for higher resolution).<\/p>\n

\"\"<\/a><\/p>\n

In my previous posts (1<\/a> and 2<\/a>) I wrote about maps with complex legends but without any kind of interactivity. In this post I show how to produce an SVG<\/a>file with interactive functionalities with the gridSVG<\/a> package.<\/p>\n

As an example, I use a dataset about population from the Spanish Instituto Nacional de Estad\u00edstica<\/a>. This organism publishes information using the PC_Axis<\/a>format which can be imported into R with the pxR<\/a> package. This dataset is available at the INE webpage<\/a> or directly here<\/a>.<\/p>\n

Let\u2019s start loading packages.<\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n
3<\/div>\n
4<\/div>\n
5<\/div>\n
6<\/div>\n
7<\/div>\n
8<\/div>\n<\/td>\n
\n
\n
library<\/code>(gridSVG)<\/code><\/div>\n
library<\/code>(pxR)<\/code><\/div>\n
library<\/code>(sp)<\/code><\/div>\n
library<\/code>(lattice)<\/code><\/div>\n
library<\/code>(latticeExtra)<\/code><\/div>\n
library<\/code>(maptools)<\/code><\/div>\n
library<\/code>(classInt)<\/code><\/div>\n
library<\/code>(colorspace)<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

Then I read the px<\/code> file and make some changes to get the datWide<\/code> data.frame<\/code>(which is available here<\/a> if you are not interested in the PC-Axis details<\/code>).<\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n
3<\/div>\n
4<\/div>\n
5<\/div>\n
6<\/div>\n
7<\/div>\n
8<\/div>\n
9<\/div>\n<\/td>\n
\n
\n
datPX <- <\/code>read.px<\/code>(<\/code>'pcaxis-676270323.px'<\/code>, encoding=<\/code>'latin1'<\/code>)<\/code><\/div>\n
<\/div>\n
datWide <- <\/code>as.data.frame<\/code>(datPX, direction = <\/code>'wide'<\/code>, use.codes=<\/code>FALSE<\/code>)<\/code><\/div>\n
<\/div>\n
provID <- <\/code>strsplit<\/code>(<\/code>as.character<\/code>(datWide$provincias), <\/code>'\u00a0 '<\/code>)<\/code><\/div>\n
provID <- <\/code>as.data.frame<\/code>(<\/code>do.call<\/code>(rbind, provID))<\/code><\/div>\n
names<\/code>(provID) <- <\/code>c<\/code>(<\/code>'code'<\/code>, <\/code>'prov'<\/code>)<\/code><\/div>\n
<\/div>\n
datWide <- <\/code>cbind<\/code>(datWide, provID)<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

Now it\u2019s time to read a suitable shapefile<\/code> (read the first post of this series<\/a> for information about it).<\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n
3<\/div>\n
4<\/div>\n<\/td>\n
\n
\n
mapSHP <-\u00a0 <\/code>readShapePoly<\/code>(fn = <\/code>'mapas_completo_municipal\/spain_provinces_ind_2'<\/code>)<\/code><\/div>\n
Encoding<\/code>(<\/code>levels<\/code>(mapSHP$NOMBRE99)) <- <\/code>\"latin1\"<\/code><\/div>\n
## The encoding must be UTF8 to be correctly displayed in the SVG file<\/code><\/div>\n
levels<\/code>(mapSHP$NOMBRE99) <- <\/code>enc2utf8<\/code>(<\/code>levels<\/code>(mapSHP$NOMBRE99))<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

Both the shapefile<\/code> and the data.frame<\/code> have to be combined using the matches between the PROV<\/code> variable of the shapefile<\/code> and code<\/code> from the data.frame<\/code>. The numeric values of the row names (mapaIDs<\/code>) will be useful in the last step.<\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n
3<\/div>\n
4<\/div>\n
5<\/div>\n<\/td>\n
\n
\n
idx <- <\/code>match<\/code>(mapSHP$PROV, datWide$code)<\/code><\/div>\n
Total <- datWide[idx, <\/code>\"Total\"<\/code>]<\/code><\/div>\n
mapSHP@data <- <\/code>cbind<\/code>(mapSHP@data, Total)<\/code><\/div>\n
mapaDat <- <\/code>as.data.frame<\/code>(mapSHP)<\/code><\/div>\n
mapaIDs <- <\/code>as.numeric<\/code>(<\/code>rownames<\/code>(mapaDat))<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

A final step is needed before calling spplot<\/code>. I will use the functions of gridSVG<\/code>to include information in the SVG file according to the characteristics of each polygon. Since gridSVG<\/code> works after the plot has been created, I need a key to identify each polygon and match it with the correspondent element of the original dataset. Unfortunately, the panel.polygonsplot<\/code> (which is the function used by spplot<\/code> when drawing polygons) does not assign a name<\/a> to each polygon. Let\u2019s create a new function (panel.polygonNames<\/code>) adding a small change in panel.polygonsplot<\/code> (you should read the full code ofpanel.polygonsplot<\/code> to understand what is happening here).<\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n
3<\/div>\n
4<\/div>\n
5<\/div>\n
6<\/div>\n
7<\/div>\n
8<\/div>\n<\/td>\n
\n
\n
panel.str <- <\/code>deparse<\/code>(panel.polygonsplot, width=500)<\/code><\/div>\n
<\/div>\n
panel.str <- <\/code>sub<\/code>(<\/code>\"grid.polygon\\\\((.*)\\\\)\"<\/code>,<\/code><\/div>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0<\/code>\"grid.polygon(\\\\1, name=paste('ID', slot(pls\\\\[\\\\[i\\\\]\\\\], 'ID'\\\\), sep=':'))\"<\/code>,<\/code><\/div>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0<\/code>panel.str)<\/code><\/div>\n
<\/div>\n
panel.polygonNames <- <\/code>eval<\/code>(<\/code>parse<\/code>(text=panel.str),<\/code><\/div>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0<\/code>envir=<\/code>environment<\/code>(panel.polygonsplot))<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

EDITED: Thanks to Andre\u2019s comment (see below), I have found thatpanel.polygonsplot<\/code> has been changed to add hole-handling. This hack will only work if this new behaviour is disabled with set_Polypath(FALSE)<\/code> before the next code chunk.
\nNow everything is ready for drawing. I use the jenks<\/a> style of theclassIntervals<\/code> function from the classInt<\/a> package to set the breaks of the color key.<\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n
3<\/div>\n
4<\/div>\n
5<\/div>\n
6<\/div>\n
7<\/div>\n<\/td>\n
\n
\n
n=7<\/code><\/div>\n
int <- <\/code>classIntervals<\/code>(Total, n, style=<\/code>'jenks'<\/code>)<\/code><\/div>\n
pal <- <\/code>brewer.pal<\/code>(n, <\/code>'Blues'<\/code>)<\/code><\/div>\n
<\/div>\n
p <- <\/code>spplot<\/code>(mapSHP[<\/code>\"Total\"<\/code>], panel=panel.polygonNames,<\/code><\/div>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0<\/code>col.regions=pal, at=<\/code>signif<\/code>(int$brks, digits=2))<\/code><\/div>\n
p<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

EDITED: You can recover the default hole-handling setting withset_Polypath(TRUE)<\/code> after the p<\/code> object has been printed.<\/p>\n

Once the plot has been created (do not close the graphic window!) thegrid.garnish<\/code> attaches SVG attributes (in this example onmouseover<\/code> andonmouseout<\/code>) to each polygon, which is identified with its name thanks to thepanel.polygonNames function<\/code>.<\/p>\n

These attributes are related to javascript<\/code> functions (showTooltip<\/code> andhideTooltip<\/code>) included in this javascript file<\/a> (this file is only a minor modification of the original file available at the webpage of the creator of grid<\/code>and gridSVG<\/code><\/a>). The grid.script<\/code> function attaches the javascript<\/code> file to the grob<\/code>and gridToSVG<\/code> produces the SVG file with a simple HTML page. EDITED: The javascript file (tooltip.js) must be saved in the same folder as the SVG and HTML files.<\/p>\n

\n
\n\n\n\n
\n
1<\/div>\n
2<\/div>\n
3<\/div>\n
4<\/div>\n
5<\/div>\n
6<\/div>\n
7<\/div>\n
8<\/div>\n
9<\/div>\n
10<\/div>\n
11<\/div>\n
12<\/div>\n
13<\/div>\n
14<\/div>\n
15<\/div>\n
16<\/div>\n
17<\/div>\n
18<\/div>\n
19<\/div>\n
20<\/div>\n
21<\/div>\n
22<\/div>\n
23<\/div>\n
24<\/div>\n
25<\/div>\n<\/td>\n
\n
\n
## grobs in the graphical output<\/code><\/div>\n
grobs <- <\/code>grid.ls<\/code>()<\/code><\/div>\n
## only interested in those with \"ID:\" in the name<\/code><\/div>\n
nms <- grobs$name[grobs$type == <\/code>\"grobListing\"<\/code>]<\/code><\/div>\n
idxNames <- <\/code>grep<\/code>(<\/code>'ID:'<\/code>, nms)<\/code><\/div>\n
IDs <- nms[idxNames]<\/code><\/div>\n
<\/div>\n
for <\/code>(id <\/code>in<\/code> unique<\/code>(IDs)){<\/code><\/div>\n
\u00a0\u00a0<\/code>## extract information from the data<\/code><\/div>\n
\u00a0\u00a0<\/code>## according to the ID value<\/code><\/div>\n
\u00a0\u00a0<\/code>i <- <\/code>strsplit<\/code>(id, <\/code>'ID:'<\/code>)<\/code><\/div>\n
\u00a0\u00a0<\/code>i <- <\/code>sapply<\/code>(i, <\/code>function<\/code>(x)<\/code>as.numeric<\/code>(x[2]))<\/code><\/div>\n
\u00a0\u00a0<\/code>dat <- mapaDat[<\/code>which<\/code>(mapaIDs==i),]<\/code><\/div>\n
\u00a0\u00a0<\/code>## Information to be attached to each polygon<\/code><\/div>\n
\u00a0\u00a0<\/code>info <- <\/code>paste<\/code>(dat$NOMBRE99, dat$Total, sep=<\/code>':'<\/code>)<\/code><\/div>\n
\u00a0\u00a0<\/code>g <- <\/code>grid.get<\/code>(id)<\/code><\/div>\n
\u00a0\u00a0<\/code>## attach SVG attributes<\/code><\/div>\n
\u00a0\u00a0<\/code>grid.garnish<\/code>(id,<\/code><\/div>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0<\/code>onmouseover=<\/code>paste<\/code>(<\/code>\"showTooltip(evt, '\"<\/code>, info, <\/code>\"')\"<\/code>),<\/code><\/div>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0<\/code>onmouseout=<\/code>\"hideTooltip()\"<\/code>)<\/code><\/div>\n
}<\/code><\/div>\n
<\/div>\n
grid.script<\/code>(filename=<\/code>\"tooltip.js\"<\/code>)<\/code><\/div>\n
<\/div>\n
gridToSVG<\/code>(<\/code>'map_with_annotations.svg'<\/code>)<\/code><\/div>\n<\/div>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<\/div>\n<\/div>\n

(Click on the image to show the SVG graphic. Move the mouse over it to display the information)<\/p>\n

\"https:\/\/procomun.files.wordpress.com\/2012\/02\/map.jpg\"<\/a><\/p>\n

<\/h6>\n","protected":false},"excerpt":{"rendered":"

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… <\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"closed","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[20],"tags":[],"class_list":["post-845","post","type-post","status-publish","format-standard","hentry","category-r"],"_links":{"self":[{"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/posts\/845","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/comments?post=845"}],"version-history":[{"count":0,"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/posts\/845\/revisions"}],"wp:attachment":[{"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/media?parent=845"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/categories?post=845"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/zhuoyao.net\/index.php\/wp-json\/wp\/v2\/tags?post=845"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}