TITLE: How much miombo is in each country
DATE: 2018-11-15
AUTHOR: John L. Godlee
====================================================================


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.

First, load some packages:

    library(rgdal)
    library(rgeos)

Next, import data on White's veg map and African countries.

  [White's veg map]: 
http://omap.africanmarineatlas.org/BIOSPHERE/pages/3_terrestrial%20v
egetation.htm
  [African countries]: 
http://maplibrary.org/library/stacks/Africa/index.htm

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

    countries <- readOGR(dsn="africa", layer="Africa")

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

  [Ryan et al. (2016)]: 
http://rstb.royalsocietypublishing.org/content/royptb/371/1703/20150
312.full.pdf

    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)

    # Add country as column
    country_miombo_stats$country <- rownames(country_miombo_stats)

I could probably spend more time to just have the function and 
lapply call give me the tidy dataframe straight off, but I don't 
have the inclination.

I think the most useful thing to come out of this little exercise 
is actually the UTM zone function, I think it's pretty neat.