A map of deforestation in Africa using R

Deforestation rates in Africa are in general high, but differences between countries are large. I want to illustrate these differences in a series of deforestation maps (change in forest and woodlands cover, including and excluding forest plantations, etc.). Normally I would use GRASS GIS or QGIS for such a task, but this should be easy to do in R too. Below I give a step-by-step description of my first attempt to make a Choropleth map in R (with the deforestation rates for 1990-2000 and 2000-2005).

Deforestation data

The Global Forest Assessment 2005 by the Food and Agriculture Organization of the United Nation (FAO) gives probably the the most comprehensive assessment of forests and forestry to date. On the site you can download a spreadsheet with various tables with information on e.g., forest cover, changes and characteristics per country. Needless to say that you should read the report before using this data. Also, have a look at Table 2 in the spreadsheet, it gives important information about how the presented numbers were estimated / measured.

The change in extent of forest and other wooded land for 1990-2000 and 2000-2005 are given in Table 4. After exporting this table as a text file (tab-delimited in the example below), it can be imported in R using read.table() .

Because the original tables in the excel file has a multiple row header, we import the table without column headers. This means that the imported data frame has the default column names V1…V27. We need the columns V5 (Annual change rate in forest cover between 1990-2000 in 1000 ha/yr) and V6 (Annual change rate in forest cover between 1990-2000 in %).

>	wd <- “/home/pvanb/Temp/”
>	FRA <- read.table(paste(wd,”/FRA.csv”, sep=””), header= FALSE, sep = “\t”, quote =””, dec =”.”,na.strings = “-“, nrows=245,skip = 5)

The data need some cleaning. More specifically, there are n.s. entries in the table. These represent not significant or very small values and can therefore be replaced by zero.

>	COUNTRY <- FRA[,1]
>	FRA <- as.data.frame(apply(FRA[,-1],2,function(x){gsub(“n.s.”,0,x)}))
>	FRA <- as.data.frame(apply(FRA,2,function(x){as.numeric(x)}))
>	FRA <- cbind(COUNTRY,FRA); rm(COUNTRY)

GADM database

The packages maps and mapdata come with global national boundary maps, the latter at an higher resolution than the former. More detailed / accurate national boundary maps are provided by the GADM database of Global Administrative Areas. The data can be downloaded per country as a shapefile, an ESRI geodatabase file, a Google Earth .kmz file or as an Rdata file (containing a SpatialPolygonsDataFrame). You can also download the data for the whole world in once, but only as shapefile or geodatabase.

Because we’ll need a vector layer for all countries in Africa, it is easier to download the shapefile for the whole world (I don’t think you can read in a geodatabase in R). In case the shapefile for the whole world is too big to import in R (it is on my computer), you first need to extract a shapefile with African countries, e.g., in QGIS (but any GIS program will do). In QGIS you need to: (1) open the GADM shapefile in QGIS file → , (2) open the attribute table and use the search function to select all records with UNREGION2 = ‘Africa’, (3) Close the table, right-click the data layer in the Layers panel and use ‘save selection as shapefile…’ to save the selection as a new shapefile ( Africa_admin.shp).

Update: it seems the gadm shapefile does not contain a column UNREGION2 (anymore). So you’ll need to find another way to extract the countries you need. For Africa it will be easy enough to select all countries manually, using the freehand selection tool in QGIS, and then save the selected polygons as a new shapefile. You can also get a shapefile with a polygon covering your region of interest and use the geoprocessing tools in QGIS to extract the countries within that region.

Next, we import the newly created shapefile in R as a SpatialPolygonsDataFrame.

>	require(maptools)
>	Africa_admin <- readShapeSpatial(paste(wd,”/Africa_admin”,sep=””), proj4string=CRS(“+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs”)

Merging the data

Both FRA and the attribute table of Africa_admin ( Africa_admin@data ) contain a column with country names. The latter actually contains more than one column with country names; we need the column NAME_FAO . We use merge() to merge the two tables into a new data frame africadat by the country name columns. Next we replace the origal attribute table of the vector map Africa_admin with the newly created data frame africadat .

Step 1 – To ensure that the order of the rows in the new table africadat is the same as in Africa_admin@data we first add a column with the row order, which will be used in the next step to re-sort the africadat table back to its original order after the merge.

Note that the original attribute table of Africa_admin contains 73 columns, most of which we do not need. We therefore only keep four columns: the new ‘order’ column and the columns with the country names, population numbers and land surface areas.

>	Africa_admin@data <- cbind(c(1:dim(Africa_admin)[1]),Africa_admin@data[,c(“NAME_FAO”,”POP2000″,”SQKM”)])
>	names(Africa_admin@data) <- c(“order”,names(Africa_admin@data)[-1])

Step 2 – Now we can merge FRA with Africa_admin@data into africadat and replace the original attribute table of Africa_admin with this newly created table.

>	africadat <- merge(Africa_admin@data,FRA, by.x=”NAME_FAO”,by.y=”COUNTRY”, all.x=TRUE)
>	Africa_admin@data <- africadat[order(africadat$order),]

Step 3 – Before we proceed we need to check if all country names of in FRA and Africa_admin match. Running the lines below shows that the names of three countries (Tanzania, DR Congo and Runion) do not match. To correct this, adapt the names of these countries in either the table FRA or in Africa_admin@data and repeat the two steps above.

>	row.names(FRA) <- FRA[,1]
>	test <- merge(africadat[, c(“order”,”NAME_FAO”)],FRA, by.x=”NAME_FAO”, by.y=0, all.x=TRUE)
>	levels(test[is.na(test[,3]),][,1][,drop=TRUE])

Plotting the data

We use spplot to plot the deforestation map for Africa. The map is written to a png file:

>	png(file=”deforestation_1.png”, width=20, height=12, units=”cm”,bg=”white”, pointsize=20, res=300)
>	spplot(Africa_admin, c(“V6″,”V8”), names.attr=c(“1990-2000″,”2000-2005″),scales=list(draw=TRUE), main=”Annual change rate in forest cover (%)”)
>	dev.off()

Here is the result:

Click on image for larger version

There are many options to further fine-tune the map to suit your needs and taste. For example, if you are only interested in the mainland, you might as well limit the map extend. You can also show deforestation rates in classes rather than on a continuous scale. To do so you need to create new columns with deforestation classes (V6a and V8a in the example below). You can define your own color scheme, or more conveniently, use one of the color schemes (RdYlGn in the example below) available through RColorBrewer().

>	require(RColorBrewer);
>	xl <- c(-19,52.5); yl <- c(-36,39)
>	qp <- c(-8,-2,-1,-0.5,0,0.5,1,2,7)
>	Africa_admin@data$V6a <- cut(Africa_admin@data$V6, qp, include.lowest=TRUE)
>	Africa_admin@data$V8a <- cut(Africa_admin@data$V8, qp, include.lowest=TRUE)
>	png(file=”deforestation_v2a.png”, width=20, height=12, units=”cm”,bg=”white”, pointsize=20, res=300)
>	spplot(Africa_admin, c(“V6a”,”V8a”), names.attr=c(“1990-2000″,”2000-2005″),scales=list(draw=TRUE, cex=0.75, col=”darkgrey”, fontfamily=”arial”), xlim=xl, ylim=yl,lw=0.1, main=”Annual change rate in forest cover (%)”, col.regions=brewer.pal(length(qp)-1,”RdYlGn”))
>	dev.off()

See below for the result. Just as a reminder, these rates are based on the changes in cover of all forests, including forest plantations.

Click on image for larger version

You can also overlay the map with other layers (raster or vector), e.g, a layer with the main road, capitals, etc. See ?spplot for all the options.

Some more examples:

There are different ways to make a Choropleth map in R. See the following links for some more examples using different methods.

• http://blog.lib.umn.edu/moor0554/canoemoore/2008/04/mapping_with_r_2.html

• http://stackoverflow.com/questions/1260965/developing-geographic-thematic-maps-with-r

• http://blog.revolution-computing.com/2009/11/choropleth-challenge-result.html