Maps with R (II)

[This article was first published on Omnia sunt Communia! » R-english, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

In my my last post 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.

First let’s load the packages.


Now, I define the geographical extent to be analyzed (approximately India and China).

ext <- extent(65, 135, 5, 55)

The first raster file is the population density in our planet, available at this NEO-NASA webpage (choose the Geo-TIFF floating option, ~25Mb). After reading the data with raster I subset the geographical extent and replace the 99999 with NA.

pop <- raster('875430rgb-167772161.0.FLOAT.TIFF')
pop <- crop(pop, ext)
pop[pop==99999] <- NA
pTotal <- levelplot(pop, zscaleLog=10, par.settings=BTCTheme)

The second raster file is the land cover classification (available at this NEO-NASA webpage)

landClass <- raster('241243rgb-167772161.0.TIFF')
landClass <- crop(landClass, ext)

The codes of the classification are described here. 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) 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!)

landClass[landClass %in% c(0, 254)] <- NA
landClass <- cut(landClass, c(0, 5, 11, 14, 16))
classes <- c('FOR', 'LAND', 'URB', 'SNOW')
nClasses <- length(classes)

EDIT: Following a question from a user of rasterVis I include some lines of code to display this qualitative variable in the map.

rng <- c(minValue(landClass), maxValue(landClass))
## define the breaks of the color key <- seq(rng[1]-1, rng[2])
## the labels will be placed vertically centered <- seq(rng[1], rng[2])-0.5

levelplot(landClass,, margin=FALSE, col.regions=terrain_hcl(4),
                        labels=classes, ## classes names as labels

This histogram shows the distribution of the population density in each land class.

s <- stack(pop, landClass)
layerNames(s) <- c('pop', 'landClass')
histogram(~log10(pop)|landClass, data=s,

Everything is ready for the map. I will create a list of trellis 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 argument of levelplot receives the correspondent at values from the global map)

at <- pTotal$legend$bottom$args$key$at

pList <- lapply(1:nClasses, function(i){
  landSub <- landClass
  landSub[!(landClass==i)] <- NA
  popSub <- mask(pop, landSub)
  step <- 360/nClasses
  pal <- rev(sequential_hcl(16, h = (30 + step*(i-1))%%360))
  pClass <- levelplot(popSub, zscaleLog=10, at=at,
                      col.regions=pal, margin=FALSE)

And that’s all. The rest of the code is exactly the same as in the previous post. If you execute it you will get this image (click on it for higher resolution).

To leave a comment for the author, please follow the link and comment on their blog: Omnia sunt Communia! » R-english. offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)