I was looking into a paper called “Determinants of woody cover in African savannas” by Sankaran et al. (2005). The paper looks at the large scale environmental factors that affect percentage woodland cover in African savanna landscapes. One figure in particular that got me interested was this one:
I haven’t fully worked out the implications of this figure yet, but what stands out to me the most is that many plots in high rainfall areas with low woody cover are classed as ‘arid fertile savanna’ by White’s veg. classification. Secondly moist infertile savanna seems to straddle the saturation point of the MAP limited woody cover.
It shows that savannas with Mean Annual Precipitation values less than ~650 mm have their upper woody cover potential limited by precipitation, but above that threshold an increase in MAP doesn’t increase the maximum potential woody cover. It also shows that lots of sites have woody cover below their MAP limited maximum, pointing to lots of other environmental factors, like fire, herbivory, soil characteristics.
Looking into this graph more, I wanted to see whether there was any biogeographic patterns that could be drawn from the data. Were all the sites with particularly low actual woody cover from a particular woodland cover biome, for example.
To do this I compared the data from Sankaran et al. (2005), which is publicly available as supplementary information, to White’s seminal vegetation classification map of 1983, which I accessed as as shapefiles from here.
I did this analysis in R, so all the code below is for R.
You can find the code and data that I used by cloning this github repo.
First, I loaded the packages and the data:
# Set working directory to the location of the source file ---- setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) # Packages ---- library(ggplot2) library(dplyr) library(rgeos) library(rgdal) # Import data ---- ## sankaran data cover <- read.csv("data/sankaran_2005_data.csv") str(cover) ## white 1983 veg data white_veg <- readOGR(dsn="data/whitesveg", layer="Whites vegetation") ## Country outline countries <- readOGR(dsn="data/africa", layer="Africa")
First I can have a go at plotting White’s map:
# Plot White's veg map data ---- ## Fortify country outline for ggplot countries@data countries_fort <- fortify(countries, region = "COUNTRY") ## Exploring whiteveg white_veg@data white_veg@polygons[] white_veg@proj4string white_veg@bbox white_veg@plotOrder ## Fortify white shape file for ggplot2 white_veg_fort <- fortify(white_veg, region = "DESCRIPTIO") names(white_veg_fort) length(unique(white_veg_fort$id)) ## Create colour palette for ggplot2 palette_veg_type_19 <- c("#FF4A46","#008941","#006FA6","#A30059","#FFDBE5", "#7A4900","#0000A6","#63FFAC","#B79762","#004D43", "#8FB0FF","#997D87","#5A0007","#809693","#FEFFE6", "#1B4400","#4FC601","#3B5DFF","#4A3B53") ## ggplot ggplot() + geom_polygon(aes(x = long, y = lat, group = group, fill = id), data = white_veg_fort) + geom_polygon(aes(x = long, y = lat, group = group, fill = NA), colour = "black", data = countries_fort) + theme_classic() + scale_fill_manual(values = palette_veg_type_19) + labs(fill = "Biome") + xlab("Longitude") + ylab("Latitude") + coord_map()
Then I had to convert the data from Sankaran et al. into a
SpatialPoints object so I can use it in future analyses.
# Create a data frame with only the latitude and longitude data cover_coords <- cover %>% select(lon_dec_deg, lat_dec_deg) # Important for later on to have lon then lat as columns # Convert to SpatialPoints object cover_spoints <- SpatialPoints(cover_coords,proj4string=CRS(proj4string(white_veg)))
The bit of this project that took me some time to work out was how to compute whether a data point from Sankaran et al. fell into a polygon of a certain type in White’s map. I ended up using
over() from the
sp package. Then I can add that data back into the original Sankaran dataset
# Add vegetation class column by referencing White map cover_veg_class <- cover %>% mutate(veg_class = over(cover_spoints, white_veg)$DESCRIPTIO)
Finally, I can use
cover_veg_class to create a ggplot of MAP vs woody cover, with the points coloured according to which of White’s vegetation classes the point is in:
# Plot with vegetation classification from White et al. 1983 ---- veg_class_plot <- ggplot(cover_veg_class, aes(x = map_mm, y = woody_cover_per, colour = veg_class)) + geom_point(size = 4) + guides(colour=guide_legend(title="Vegetation Class")) + xlab("MAP (mm)") + ylab("Woody Cover (%)")