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

# 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, 
               field = "binomial",
               fun = function (x, ...) length(unique(na.omit(x))))

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, ],
                                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)

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, 2019 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, 2019 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, 2019 at 14:32
  • Damn again, my bad. I corrected it.
    – Farewell
    Jan 7, 2019 at 14:38
  • Its all working and I can reproduce your differences now. Explanation might take a bit longer...
    – Spacedman
    Jan 7, 2019 at 16:49

1 Answer 1


Thank you very much for reporting this. I have looked at this in detail, and in this case, method 2 produced the correct answer, whereas method 1 did not. This was related to polygons with holes and the order in which they were rasterized. I have fixed this in raster version 2.8-14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.