Mapping species richness in cells based on georeferenced species records

Example of one relatively simple way to map grid cell-based species richness starting with individual georeferenced species records. This is only a method for processing the data, and does not consider sampling biases, for example. There are multiple ways to work between data frames/point data and rasters in R, this is just one. You can also 'rasterize' point data or assign values to a pre-generated raster layer etc.

The basic premise below is to take a set of species records with their geographic (longlat) coordinates, truncate the coordinates so that the records fall within a regular grid pattern (here 0.1 degrees resolution) and then for each point on the truncated grid, count the number of unique species recorded and display as a grid. A moving or sliding (neighbourhood) window can be used to even out the scores a little.

##############################

require(raster)

n <- c(1:2000); for (i in 1:2000) {n[i] <- sample(as.factor(1:100),1)} #create a random set of 100 'species' (numnbers 1 to 100 as factor levels, each representing a species), i.e. pick one species 2000 times

spp_records <- data.frame(latitude=runif(2000, min=-33, max=-30), longitude=runif(2000, min=135, max=140), value=n) #create a data frame simulating 2000 spatial records of the 100 simulated species within a predetermined spatial area - numbers under 'value' are in this case factors representing unique species from the above pool of 100 

#in reality the above step would be replaced by the raw data - a data frame latitude, longitude and a column for the name of the species recorded at those coordinates.

coordinates(spp_records) <- c(1,2) #set to spatial points data frame object by highlighting the latitude and longitude columns

plot(spp_records, pch=20, cex=0.5) # plot as points to visualise where the records were made

spp_records <- data.frame(latitude=round(spp_records$latitude, digits=1), longitude=round(spp_records$longitude, digits=1), value=n) #now recreate, but truncate the coordinates to a regular grid at 0.1 degrees resolution

coordinates(spp_records) <- c(1,2)

rich_count <- function(x) length(unique(x)) # simple function to count the number of unique 'species' (factor levels)

richness <- aggregate(spp_records$value, list(spp_records$latitude, spp_records$longitude), rich_count) #count the number of unique species that have the same truncated coordinates

coordinates(richness) <- c(2,1) #set coordinates

richness2 <- rasterFromXYZ(richness) #create a gridded map (at 0.1 resolution) with values as the number of unique species recorded in that cell

projection(richness2) <- "+proj=longlat +ellps=WGS84" #set the projection for plotting

plot(richness2, col=rev(heat.colors(max(richness$x)))) # cell-based 'species richness' map

Cell-based count of species richness calculated from simulated set of individual georeferenced species records, using 'raster' package in R

plot(focal(richness2, w=matrix(1,3,3), fun=mean, na.rm=TRUE), col=rev(heat.colors(max(richness$x)))) #you can also play around with a moving window type add-on, in this case averaging the richness scores for a given cell using the surrounding set of cells, to even out idiosyncrasies of how records are assigned to one cell or another

The same map of cell-based (simulated) species richness as above with a sliding window average to even out scores among neighbouring cells