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[[1]]
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 (%)")
