1

I was trying to calculate species richness from IUCN distribution data in R, and I decided to check my results with two different methods, both of which rely on rasterize function from the raster package.

  1. The simplest method consists in counting the number of species in each grid cell by using the length function in rasterize (it is actually described in rasterize help as an example to be used to calculate richness).

  2. A slightly more complex method consist in rasterizing each species polygon into presence-absence rasters, and then summing all rasters to obtain a species richness raster.

However, I ran into minor inconsistencies in the results, and I am not able to figure out why.

I paste below the code I used for these two methods. The dataset I used can be downloaded from the IUCN website here, download the Tailed amphibians dataset (25MB). The IUCN requires an account to download the data, I am sorry about that but I cannot find another way to make a reproducible example of this dataset (and I am not allowed to redistribute it).

Note that in this dataset, there can be multiple polygons for each species. Species names are in the column binomial of the IUCN dataset.

Data preparation

# Packages
library(raster)
library(sf)

# Reading the polygons
iucn_sf <- st_read("./data/IUCN data/CAUDATA.shp")
# Some data filtering
iucn_sf <- iucn_sf[which(iucn_sf$presence == 1), ]
iucn_sf <- iucn_sf[which(iucn_sf$origin == 1), ]
iucn_sf$binomial <- droplevels(iucn_sf$binomial)

# Creating a raster at the desired resolution
r <- raster(nrows=180, ncols=360, xmn=-180, xmx=180, ymn=-90, ymx=90, 
            crs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0", 
            resolution = 1, vals=NULL)

Method 1: rasterize species polygons into richness raster with length

This method is described in the rasterize documentation.

a <- rasterize(iucn_sf, 
               r,
               field = "binomial",
               fun = function (x, ...) length(unique(na.omit(x))))
plot(a)

Results of method 1

The maximum richness with method 1 is 26.

Method 2: rasterize each species into a presence-absence raster and then calculate the sum

pa.stack <- stack()
for(sp in levels(iucn_sf$binomial))
{
 # We rasterize the polygons of the current species into values of 1
 pa.stack <- addLayer(pa.stack,
                      rasterize(iucn_sf[iucn_sf$binomial == sp, ],
                                r, 
                                field = "binomial", 
                                fun = function(x, ...) 1))
}
# After the loop, we sum all presences to obtain the species richness
b <- sum(pa.stack, na.rm = T)
plot(b)

Results of method 2 The maximum richness with method 2 is 25.

Differences between methods 1 and 2:

plot(a - b)

Differences between method 1 and 2

You can see that most of the results are the same, except for a few pixels where there can be a difference of 1 or 2 species. I do not understand where this difference comes from (is this a mistake in my code or a difference in the procedure?), which is why I turn to the wider community to seek advice.

I used R version 3.5.2, packages sf version 0.7-2 and raster version 2.8-4.

  • I've got your first example to run, but the second one references cur.distrib which you don't define. Not sure if this is the problem (maybe cur.distrib isn't quite the same as iucn_sf) but am stuck with that at the moment. – Spacedman Jan 7 at 13:58
  • Oops sorry, I made a modification on my script while writing the message and forgot to report it here. It is corrected. – Farewell Jan 7 at 14:09
  • Now I get Error in .polygonsToRaster(x, y, field = field, fun = fun, background = background, : when rasterizing multiple fields you must use "fun=first" or "fun=last" in the rasterize for method 2... Same raster version. – Spacedman Jan 7 at 14:32
  • Damn again, my bad. I corrected it. – Farewell Jan 7 at 14:38
  • Its all working and I can reproduce your differences now. Explanation might take a bit longer... – Spacedman Jan 7 at 16:49

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.