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:

Deforestation rates in Africa (version 1)

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.

Deforestation rates in Africa (version 2)

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.

About these ads

11 thoughts on “A map of deforestation in Africa using R

  1. Luca De Benedictis

    Dear Paulo van Breugel,
    thanks for the nice post.

    I am an international trade economist and I am starting to use R to do some statial analysis on Africa.

    I was trying to replicate your example but I was unable to do it.

    I think there is a mistake in the secon line of the script
    FRA <- read.table(paste(wd,”/FRA.csv”, sep=”"), header= FALSE, sep = “\t”, quote =”", dec =”.”,na.strings = “-”, nrows=245,skip = 5)

    that should be
    FRA <- read.table(paste(wd,”/FRA.csv”, sep=”"), header= FALSE, sep = “;”, quote =”", dec =”.”,na.strings = “-”, nrows=245,skip = 5)

    But my main problem is with the shape file Africa_admin.

    I was unable to download it from QGIS (but I am a dummy on that!!!). Can you make the Africa_admin.shp file available (send it to me) or even write a post on how to use QGIS to extract data for R.

    Thanks a lot.
    All the best

    Luca

    Reply
    1. pvanb Post author

      It works for me. Be careful to use straight quotation and double quotation marks. As for your suggested change, the ‘sep’ parameter in ‘read.table’ is to define the field delimiter. In my post I am importing a tab-delimited text file. You would use your line to import a text file with a “;” as field delimiter. I think that in Excel you can choose between ‘csv’, which is comma delimited by default, or tab delimited. In some versions of Open Office, the default text delimited is “;”.

      The Africa_admin.shp was extracted from the GADM database and saved as a shapefile using the three steps given in the post. I might write a short post how to do this in more detail at a later stage, when time allows. As I mentioned, you can of course use any GIS software you are comfortable with to do the same. If you are more comfortable in R, you really should try the ‘maps’ or ‘mapdata’ packages. Redistribution of the GADM data is not allowed, but I guess sending you the file personally will not be a problem so I’ll do that.

      Reply
  2. Pingback: Ecostudies

  3. Clovis

    Hello,
    Thanks for your site, I find many useful information!
    I want to test your exercices using R/QGis/Grass on deforestation.
    However i’m struggling to prepare the FRA dataset (problems with pdf copie paste operation).
    It there a another way to get it? (send it via mail for instance)
    Thanks
    Clovis

    Reply
  4. Celeste

    Dear Mr. van Breugel,
    I am a Peace Corps volunteer serving in Southwestern Ethiopia doing Conservation and Natural Resource work. I am hoping to help the Woreda (district) Agricultural Office map the rate of deforestation in the area. I am new to mapping– although, I have a bit of experience programming in Python and Java. I’ve begun to read your blog, but it’s pretty high level. Do you have any recommendations on where to start?
    Thanks very much!
    Celeste

    Reply
    1. pvanb Post author

      Hi Celeste,

      That is a bit hard to answer without knowing a bit more about what you want to do. First question would be what data do you have or are planning to collect? Depending on the type of data, you can decide what methods and tools (software) to use. And are you planning to monitor / map the forest extend at regular intervals (e.g., once a year), or is this a one-time assessment of the rate of deforestation?

      In any case, the data presented in this post is not suitable at all for what you want, it only gives data at a national level.

      I have worked with botanists in Ethiopia and Denmark to make a potential vegetation map of Ethiopia. Potential in this case means that vegetation that would become established under current conditions if you ban all human interferences. In other words, it does not give the actual vegetation or land use cover. But it is very useful as a baseline to which you can compare current land use or vegetation distribution, which is for which we are using it. Our problem is probably the same in that respect, the lack of accurate up-to-date forest cover data (and changes therein). If you have data or can collect field data on deforestation in the region, that could be interesting to compare with our work.

      As for the software, for display and simple calculations of deforestation rates, QGIS could be very suitable and is relatively easy to use (and it is free, open source software). At the same time it offers the GRASS GIS toolbox with more advanced functions (including remote sensing capabilities, but that is not really my domain).

      Reply
  5. annoporci

    Hi, thanks for the tutorial! I’m experimenting, great stuff!

    I’m having problems with extracting the African shapefiles. I have downloaded the latest GADM (gadm_v2_shp.zip) and have QGIS 2.0.

    You write “open the attribute table and use the search function to select all records with UNREGION2 = ‘Africa’” The problem is that I don’t see the UNREGION2 variable. So I managed to pull the “attribute table” via the layer menu, but I can’t find a way to select “Africa”.

    (If I do search in name0, then I pull “Central African Republic” and “South Africa”, which is not the whole thing by a long way)

    Do you have any idea what the problem might be?

    Also, I’d appreciate if you’d give a hint on how to “look” inside of Africa_admin@data, to see what available variables there are, I know of head(), is there another way you use?

    I’m downloading QGIS 2.2 to see if the update helps at all, but I’d have thought the problem would be rather in the missing UNREGION2 in GADM 2, wouldn’t you think?

    It would be really cool if gadm allowed to download “by continent”…

    Reply
    1. pvanb Post author

      Hi, I checked and in the latest version there is indeed no column ‘UNREGION2′. I am not sure anymore what version I have been using. I’ll look into this when time allows (that will not be very soon). You can of course always select the countries manually. I would first use the ‘freehand selection’ option (in the map navigation toolbar) to make a first selection and then, if needed, get a map with the outlines of Africa (you’ll find plenty of links if you Google for it) and use the geoprocessing tools in QGIS to make the final selection.

      The Africa_admin@data is just a data frame, so to see what variables (columns) it contains, you can use names(Africa_admin@data), or str(Africa_admin@data)

      You can always send the people behind gdam a suggestion to make the files available per continent. There is a contact form on their web page.

      Reply
      1. annoporci

        brilliant, thanks a lot! (though it doesn’t help right away, I’m glad it wasn’t something totally stupid I was doing). I’ll try this: get a list of all African countries and pull them with a script, I think it’s feasible. Me too don’t have time to do it right now. Thanks again!

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s