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 %).

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.

BothFRAand the attribute table ofAfrica_admin(Africa_admin@data) contain a column with country names. The latter actually contains more than one column with country names; we need the columnNAME_FAO. We usemerge()to merge the two tables into a new data frameafricadatby the country name columns. Next we replace the origal attribute table of the vector mapAfrica_adminwith the newly created data frameafricadat.

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 theafricadattable 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.

Step 3 – Before we proceed we need to check if all country names of inFRAandAfrica_adminmatch. 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 tableFRAor inAfrica_admin@dataand repeat the two steps above.

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().

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.

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

On the FRA site, you can download the data in spreadsheet (excel) format, so no need to copy/paste from pdf files. Also note that you can now download the FRA 2010 from http://www.fao.org/forestry/fra/fra2010/en/, also in spreadsheet format.

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

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

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?

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.

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!