A while back, a colleague from Cuba contacted me seeking help with making species richness maps for plants. I had written about species richness maps in R before, but only when working with point occurrence data or species range polygons. In their case, the task was to reclassify and sum a bunch of MaxEnt models to create a species richness layer.
I tried to point them to existing tutorials, but didn’t find any recent ones with good exposition, so here we are. I haven’t run my own SDMs in years, but as far as I know, this should still be the overall process.
Run a distribution model (repeat for n species)
Pick a threshold and reclassify each species’ raster layer with continuous values of suitability or probability of occurrence into binary presence/absence values.
Sum the reclassified layers to count the total species richness for each cell
Disclaimer: I don’t know how others are doing this recently (I suspect with QGIS or ESRI software).
Thanks to all the ongoing development happening with spatial methods in R, a stars + sf + tidyverse approach to raster calculations and plotting is now possible, and below is my take on it.
Proof of concept
Let’s start with a proof of concept. Before working with real data I wanted to see if the approach I cooked up works. For this, we’ll create toy rasters with random suitability values.
The function below creates a 4x5 matrix with random values, converted to a stars object with a silly name created by ids::adjective_animal.
We’ll run this 8 times with purrr::rerun, stack the objects with c(), then combine them with merge().
stars objects have attributes and dimensions, which are the cell values, and the metadata for the dimensions in the array, respectively. We can set these names by passing character vectors to setNames (for attributes) and st_set_dimensions (for the xy coordinates and the species, in this case).
Printing the object gives us this:
This set of layers is ready for plotting with ggplot2
To reclassify all of the arrays, we can use tidyverse methods. Here we consider every cell value >70 as present (1) and <70 as absent (0). It’s very cool how we can use dplyr verbs directly on the stars attributes.
We can also plot the reclassified arrays as a faceted ggplot2 plot.
To sum these array dimensions, we use st_apply for the raster calculation. We apply sum to the dimensions with the xy coordinates to get the total value for each pixel.
This is our richness raster, which we can plot now. The math checks out (i.e. the pixel values in this layer are the sum of all the 0/1 arrays)
We can plot the stacks in “3d” with this approach by soil scientist Lauren O’Brien. After converting the rasters (stars) to polygons (simple features) and reshaping them, we can warp the geometries directly with dplyr::mutate.
We can use patchwork to plot the original and reclassified stacks side by side