~ Sobre software, documentación y ciencia libres

Maps with R (II)

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.

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.

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

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

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)