# How much miombo is in each country

## 2018-11-15

For some quick and dirty statistics to quote in the introduction to a report, I wanted to know how much of the miombo biome (as defined by White’s vegetation map) was in Angola. Afterwards, I decided to try and apply the same methods to all the countries in southern Africa. I used R to do the analyses.

``````library(rgdal)
library(rgeos)
``````

Next, import data on White&rsquo;s veg map and African countries .

``````white_veg <- readOGR(dsn = "whitesveg", layer = "Whites vegetation")

``````

The Coordinate reference system (CRS) isn’t explicitly defined in either of the spatial objects, but it’s a good guess that they will be WGS84 long-lat, so let’s add that.

``````proj4string(white_veg) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
proj4string(countries) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
``````

I’m not interested in the other biomes defined in White’s veg. map, only those that make up the “miombo”. I use the miombo definition used in Ryan et al. (2016) .

``````miombo <- white_veg[which(
(white_veg\$DESCRIPTIO == "Moist-infertile savanna") |
(white_veg\$DESCRIPTIO == "Mosaics of forest") |
(white_veg\$DESCRIPTIO == "Mopane savanna") |
(white_veg\$DESCRIPTIO == "Montane Forest") |
(white_veg\$DESCRIPTIO == "Hydropmorphic grassland") |
(white_veg\$DESCRIPTIO == "Arid-fertile savanna") |
(white_veg\$DESCRIPTIO == "Sedge and reed swamp")),]
``````

Then, I only want to keep countries that contain miombo, which I can do using `gIntersects()` from rgeos.

``````country_list <- split(countries, countries\$COUNTRY)

intersections <- lapply(country_list,
function(x)
as.vector(unlist(gIntersects(x, miombo))))

country_list_miombo <- country_list[intersections == TRUE]
``````

Then, I define a function to find UTM zones, because countries vary in their UTM zone, so I can’t just use a flat value. Converting from lat-long to UTM is necessary so I can get km^2 area estimates, rather than square degrees. The function takes any spatial object and uses the bounding box to estimate the UTM zone.

``````# Define a function to find the UTM zone
utm.zone <- function(x){
num <- floor(((mean(x@bbox[1,]) + 180)) / 6) + 1
let <- ifelse(mean(x@bbox[2,]) > 0, "N", "S")

return(paste0(num, let))
}
``````

Then, time for a big lengthy, possibly overly messy function to return a list of miombo area stats for each country.

``````miombo.country.perc <- function(country, miombo){
country_fix <- gBuffer(
country,
byid = TRUE,
width = 0)

miombo_country <- gIntersection(
country_fix,
miombo,
byid = TRUE,
drop_lower_td = TRUE)

miombo_utm <- spTransform(
miombo,
CRS(paste0("+proj=utm +zone=", utm.zone(miombo), " ellps=WGS84")))

miombo_country_utm <- spTransform(
miombo_country,
CRS(paste0("+proj=utm +zone=", utm.zone(miombo_country), " ellps=WGS84")))

country_utm <- spTransform(
country_fix,
CRS(paste0("+proj=utm +zone=", utm.zone(country_fix), " ellps=WGS84")))

area_miombo_km2 <- gArea(miombo_utm) / 1e6

area_miombo_country_km2 <- gArea(miombo_country_utm) / 1e6

area_country_km2 <- gArea(country_utm) / 1e6

perc_miombo_country <- area_miombo_country_km2 / area_country_km2 * 100

perc_miombo_all <- area_miombo_country_km2 / area_miombo_km2 * 100

return(data.frame(area_miombo_country_km2,
area_country_km2, perc_miombo_country, perc_miombo_all))
}
``````

Then I need to run the function on each country in the list of countries, using `lapply()`.

``````country_miombo_stats <- lapply(country_list_miombo, miombo.country.perc, miombo = miombo)
``````

And finally there is a bit of cleaning up to get a tidy data frame.

``````# Collapse the resulting list
country_miombo_stats <- do.call("rbind", country_miombo_stats)