Below, I present a brief tutorial applying the methodology of Vásquez-Restrepo y Daza (2025) for assessing extinction risk according to IUCN Criterion B, using forest cover change as a proxy for habitat quality reduction. The intention is to facilitate the application of these analyses, as although they are simple, not everyone knows how to implement them programmatically. The goal is to make our categorizations a bit more objective, given that the application of subcriteria a, b, and c for geographic categorization often lacks evidence and is, in many cases, based solely on assumptions about habitat degradation.

All the files needed to replicate this tutorial can be downloaded here, except for the forest loss data, which are not included due to their size. However, these data can be downloaded via the links provided later in the tutorial when needed.

Before getting started

The first thing we need to do is load all the necessary libraries:

library(dplyr)
library(gdalUtilities)
library(geosphere)
library(raster)
library(rgdal)
library(sf)
library(terra)

Additionally, we’ll define our working directory:

setwd("D:/Daniel/Desktop/IUCN")

Next, we’ll load a couple of custom functions to perform some of the calculations and make our work easier. These functions should be run as is, without needing to modify them (unless you know what you’re doing). To keep this tutorial simple, I won’t go into detail about how each line works; for that, I’ve included some comments in the code. In short, just copy and paste them into your script or console.

source("scripts/functions.R")

Loading the data

The next step is to load the occurrence data. For this example, I’ll use the dataset from Vásquez-Restrepo y Daza (2025), where we describe a new species of the genus Echinosaura and use the methodology explained in this tutorial to assess the threat category of all species in the genus.

# Read the data
data <- read.csv("Echinosaura.csv")

# Filter unique occurrence records
data <- data %>% distinct(Species, Latitude, Longitude, .keep_all = TRUE)

The data format should look like the one shown below:

##            Species Latitude Longitude
## 1 E. brachycephala -0.43333 -78.96667
## 2 E. brachycephala -0.41667 -78.80000
## 3 E. brachycephala  0.11778 -78.60756
## 4 E. brachycephala  0.09230 -78.62006
## 5 E. brachycephala -0.43333 -78.95000
## 6 E. brachycephala -0.09300 -78.98959

The table above only shows one species for ease of visualization, but it contains data for nine species, so there is no limit on the number of species that can be processed simultaneously.

Generating a mosaic

For this methodology, we will use forest cover and loss data from Hansen et al. (2013) v1.10, which can be downloaded from the Global Forest Watch website. This dataset has a resolution of 30 m, covering the period 2001–2022, using the year 2000 as a baseline. Additionally, it’s important to know that in this context, “forest” is defined as any vegetation taller than 5 meters.

Because this dataset is very large due to its resolution, it must be downloaded in tiles. That’s why we need to create a mosaic and reassemble the pieces. For this example, I’m using the tiles highlighted in red, covering from 20°N, 90°W to 0°N, 80°W. The files we need are those labeled as treecover2000 and lossyear, meaning both for each tile.

The 2000 cover data indicate the proportion of vegetation per pixel for that year, while the forest loss data indicate the year in which a pixel transitioned from forest to non-forest, via a numeric value from 0 (no loss) to a value between 1 and 22 (year).

Once the rasters have been downloaded, we’ll merge them using the first custom function we loaded at the beginning. For this, it’s important to have each raster group (treecover2000 and lossyear) in separate subfolders.

# Merge all Hansen tiles into a single virtual raster mosaic
vrt_treecover2000 <- virtual_raster(path = "treecover2000", output = "mosaic_2000.vrt")
## Virtual raster created at: mosaic_2000.vrt
vrt_loss <- virtual_raster(path = "lossyear", output = "mosaic_loss.vrt")
## Virtual raster created at: mosaic_loss.vrt

The reason for creating a virtual raster instead of a true mosaic is computation time. Creating and saving a full mosaic could take hours depending on your computer’s capacity and the data resolution. A virtual raster is just a file that references the existing rasters and lets us perform calculations on them as if they were a single file, without physically merging them.

Minimum convex polygons

We will now generate minimum convex polygons for each species, which are equivalent to the extent of occurrence (EOO) used by the IUCN. These polygons will be used to crop our mosaic to reduce processing time, and we will also calculate the EOO area from them.

# Convert data to an sf object
# --- 4326 is the EPSG code for WGS84
data_sf <- st_as_sf(data, coords = c("Longitude", "Latitude"), crs = 4326)

# Group by species and create convex hulls
hulls <- data_sf %>%
  group_by(Species) %>%
  summarise(geometry = st_combine(geometry)) %>%
  mutate(geometry = st_convex_hull(geometry)) %>%
  st_as_sf()

Calculating forest loss

Because the forest loss data indicate a year rather than a proportion, the next step is to find a way to transform them. To do this, we can aggregate the original raster, that is, group pixels to reduce its resolution and use the cells within each new pixel to generate a metric of forest cover change.

To do this, we’ll first define an aggregation factor. In this example, we’ll go from 30-meter to 1-km pixels:

# Define the aggregation factor
# --- output resolution = aggregation factor × input resolution
# --- aggregation factor = output resolution / input resolution
# --- At the equator, 1 degree ≈ 111 km, so for a 1 km output we use 1 / 111
fact <- (1/111) / res(vrt_loss)[1]

Next, we’ll convert the previously generated convex hulls to a spatial object (necessary for some package operations) and also load a continental outline. The latter will be used for a second masking step, since some polygons may extend into the ocean, which would lead to an overestimation in this example, as we’re dealing with terrestrial species.

# Convert sf hulls to a Spatial object
hulls_sp <- as(hulls, "Spatial")

# Load continental contour
contour <- readOGR(dsn = "shapes", layer = "America_MERGED")

# Some people may experience issues when reading the contour.
# This happens because readOGR is a function from the rgdal package,
# which is no longer available in the official R repositories and can be
# somewhat difficult to install manually.
# Therefore, this function can be replaced with vect from the terra package,
# followed by a transformation to the appropriate class needed.
contour <- vect("shapes/America_MERGED.shp")
contour <- as(contour, "Spatial")

Finally, we’ll create two lists to store the results of the aggregation operations.

# Lists to store the clipped rasters and their aggregated values
aggregated_list_loss <- list()
aggregated_list_2000 <- list()

Then, we’ll loop through the list of species and for each one generate two new aggregated rasters. For treecover2000, we’ll assign the new cell the average value of all pixels within the aggregation. For lossyear, we’ll calculate the proportion of pixels with values greater than or equal to 1 over the total number of aggregated pixels. In the second case, we do this because we only have information about the year when forest turned into non-forest; in other words, we can treat the values as binary (0/1) and calculate how many cells lost forest cover during a given period. Put differently, we calculate the average forest cover in the year 2000 (FC00), forest loss between 2001 and 2022 (FL), and the current remaining forest cover (FC22) within each species’ EOO as follows:

\[ \text{FC}_{22} = \text{FC}_{00} - (\text{FC}_{00} \times \text{FL}) \] We also define a very conservative threshold of 95% and a more relaxed one of 80% to assess the reduction in habitat quality. Specifically, if current forest cover exceeds 95% of the 2000 value, we consider the impact of forest loss on habitat quality to be low; if it’s between 80% and 95%, we consider it moderate (concerning); and if it’s below 80%, we consider it high. We use the year 2000 as a baseline because landscapes naturally include non-forested areas that are not anthropogenically altered, so it doesn’t make sense to expect the entire area to be forested.

# Loop over each species polygon
for(i in 1:length(hulls_sp)){
  
  species_name <- hulls_sp@data$Species[i]
  cat("Processing:", species_name, "\n")
  
  # Get the polygon for species i
  poly <- hulls_sp[i, ]
  
  # Crop the VRT mosaics to the species polygon
  clipped_loss <- crop(vrt_loss, poly)
  clipped_2000 <- crop(vrt_treecover2000, poly)
  
  # Aggregate the clipped rasters
  agg_loss <- aggregate(clipped_loss, fact = fact, fun = prop, cores = parallel::detectCores())
  agg_2000 <- aggregate(clipped_2000, fact = fact, fun = mean, cores = parallel::detectCores())
  
  # Mask using the species polygon and the continental contour
  agg_loss <- mask(agg_loss, poly) %>% mask(., contour)
  agg_2000 <- mask(agg_2000, poly) %>% mask(., contour)
  
  # Store results in the corresponding lists
  aggregated_list_loss[[species_name]] <- agg_loss
  aggregated_list_2000[[species_name]] <- agg_2000
}
## Processing: E. brachycephala 
## Processing: E. centralis 
## Processing: E. embera 
## Processing: E. fischerorum 
## Processing: E. horrida 
## Processing: E. keyi 
## Processing: E. orcesi 
## Processing: E. palmeri 
## Processing: E. panamensis
# Create a new list by applying the forest cover calculation
forest_cover <- Map(forestCover, aggregated_list_2000, aggregated_list_loss)

Summary table

Finally, we just need to generate a summary table with all the information we’ve obtained. For this, we’ll use another custom function loaded at the beginning. By default, this function assumes corrected = FALSE, meaning it calculates the EOO area as is, without considering any additional contours. In this case, as previously mentioned, since the species are terrestrial, we use the continental outline to reduce overprediction. Additionally, the number of localities was calculated using a minimum distance of 18.5 km (approximately 10 arcmin).

generate_summary(hulls = hulls_sp,
                 data = data,
                 aggregated_list_2000 = aggregated_list_2000,
                 aggregated_list_loss = aggregated_list_loss,
                 forest_cover = forest_cover,
                 corrected = TRUE,
                 contour = contour,
                 min.dist.km = 18.5)
##            Species   EOO Localities  FC.2000 Threshold    prop.FL FC.current
## 1 E. brachycephala  2835          5 83.07404  78.92034 0.02449480   81.09095
## 2     E. centralis 58148         21 66.51405  63.18835 0.10244575   58.52673
## 3        E. embera 46671         17 89.72440  85.23818 0.04637628   85.69461
## 4   E. fischerorum   320          2 92.73193  88.09533 0.02229771   90.69385
## 5       E. horrida 29469         21 85.18532  80.92606 0.06656784   79.47579
## 6          E. keyi  6552          8 84.74488  80.50764 0.03100361   82.33109
## 7        E. orcesi  7883          5 96.86037  92.01736 0.02540984   94.47404
## 8       E. palmeri 28907         13 90.61501  86.08426 0.01679612   89.20893
## 9    E. panamensis 12731          7 75.71753  71.93165 0.08721947   68.67376
# EOO: extent of occurrence in km²  
# Localities: number of unique localities within an 18.5 km radius  
# FC.2000: average proportion of forest cover in the year 2000  
# Threshold: forest cover proportion used as a threshold for habitat loss  
# prop.FL: average proportion of forest loss during the time period  
# FC.current: average proportion of forest cover in the final year of data

# The data in this table may differ slightly from those published in our
# article, as this process was originally carried out using a combination
# of R and QGIS, as well as a thinned dataset, where the points forming the
# minimum convex polygons may vary

With this table, we can now refer to the IUCN guidelines to perform our assessment according to Criterion B, since we have the EOO, number of localities, and a forest cover threshold as a proxy for an estimated, inferred, or projected decline in habitat quality.

For example, for our new species Echinosaura embera, we have an EOO greater than 20,000 km², so Criterion B1 is not met, and the species does not qualify to be considered threatened. Likewise, it has more than 10 known localities, and the forest cover within its current EOO is above 95% of what was available in 2000. Therefore, its category would be Least Concern (LC).

On the other hand, Echinosaura keyi, currently listed as Vulnerable (VU), has an EOO greater than 5,000 km² but less than 20,000 km², thus meeting the first vulnerability condition. It also has 8 known localities, placing it again in the Vulnerable category. Now, for the species to be considered as such, it must meet two of the three conditions (a, b, and c). It already meets one in terms of number of localities, so we would need to justify some type of observed, estimated, inferred, or projected decline in one of the parameters indicated by the guidelines. In this case, we’ll use habitat quality, since that’s why we calculated forest loss. Referring to our table, we see that Echinosaura keyi still inhabits an area that retains, on average, more than 95% of the cover it had in 2000, so we consider the impact of deforestation to be low, and the second condition not applicable. Therefore, instead of Vulnerable, we consider it Near Threatened (NT), as if habitat conditions worsen in the near future, the species could become threatened.

Lastly, if we look at Echinosaura panamensis, this species is currently considered Least Concern (LC); however, our data show it has an EOO less than 20,000 km², fewer than 10 known localities, and an average forest cover below 95% of the proportion in the year 2000. Therefore, this species should be considered Vulnerable (VU).