A continuación, presento un pequeño tutorial aplicando la metodología de Vásquez-Restrepo y Daza (2025) para la evaluación del riesgo de extinción según los criterios B de la UICN, utilizando el cambio en la cobertura de bosque como un proxy de la reducción en la calidad del hábitat. La intención es facilitar la aplicación de estos análisis, ya que, si bien son simples, no todo el mundo sabe cómo implementarlos a nivel programático. El propósito es intentar que nuestras categorizaciones sean un poco más objetivas, dado que la aplicación de los subcriterios a, b y c de la categorización geográfica no siempre cuenta con evidencia y, en muchos casos, se basa únicamente en supuestos sobre un deterioro en la calidad del hábitat.

Todos los archivos necesarios para replicar este tutorial pueden descargarse aquí, con excepción de los datos de pérdida de bosque, los cuales, debido a su tamaño, no se incluyen. Sin embargo, pueden descargarse a través de los enlaces que se proporcionan más adelante, cuando sean requeridos.

Antes de iniciar

Lo primero que tenemos que hacer es cargar todas las librerías necesarias:

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

Adicionalmente, vamos a definir nuestro directorio de trabajo:

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

A continuación, vamos a cargar un par de funciones personalizadas para realizar algunos de los cálculos y facilitarnos el trabajo. Estas funciones deberían ejecutarse tal cual, sin necesidad de modificarlas (a menos que sepas lo que estás haciendo). Para no complicar demasiado este tutorial, no entraré en detalles sobre el funcionamiento de cada línea; para ello, incluí algunos comentarios en el código. En pocas palabras, solo copia y pega en tu script o consola:

source("scripts/functions.R")

Cargando los datos

Lo siguiente en nuestra lista es cargar los datos de ocurrencia. Para este ejemplo, voy a utilizar los datos de Vásquez-Restrepo y Daza (2025), en los que describimos una nueva especie del género Echinosaura y utilizamos la metodología que se explica en este tutorial para evaluar la categoría de amenaza de todas las especies del género.

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

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

El formato de los datos debe ser como el que se muestra a continuación:

##            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

En la tabla anterior se muestra solo una especie para facilitar su visualización, pero esta contiene información para nueve, por lo que no hay un límite en cuanto al número de especies que se pueden procesar simultáneamente.

Generando un mosaico

Para esta metodología vamos a utilizar los datos de cobertura y pérdida de bosque de Hansen et al. (2013) v.1.10, los cuales se pueden descargar desde el sitio web de Global Forest Watch. Este conjunto de datos tiene una resolución de 30 m, para el período 2001–2022, tomando el año 2000 como punto de referencia. Por otro lado, es importante saber que, en este contexto, se define como “bosque” a toda vegetación con una altura superior a 5 m.

Ya que este conjunto de datos es sumamente pesado debido a su resolución, es necesario descargar la información por bloques. Es por eso que tenemos que crear un mosaico y ensamblar las piezas de nuevo. Para este ejemplo, estoy usando los recuadros que aparecen en rojo, los cuales abarcan desde 20°N, 90°W hasta 0°N, 80°W. Los archivos que necesitamos son los que se indican como: treecover2000 y lossyear, es decir, ambos para cada bloque.

Los datos de cobertura en el 2000 indican la proporción de cobertura vegetal por píxel para ese año, mientras que los de pérdida de bosque indican el año en que se produjo el cambio de un estado de bosque a uno sin bosque, a través de un valor numérico que va de 0 (sin pérdida) hasta un valor entre 1 y 22 (año).

Una vez descargados los rásters, procederemos a juntarlos usando la primera de las funciones personalizadas que cargamos al comienzo. Para esto es importante tener cada grupo de rásters (treecover2000 y lossyear) en subcarpetas diferentes.

# 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

La razón para crear un ráster virtual y no un mosaico verdadero es el tiempo de cómputo. Crear un mosaico y guardarlo podría tomar horas, dependiendo de la capacidad de la computadora que estemos usando y de la resolución de los datos. En cambio, un ráster virtual no es más que un archivo que referencia los rásters que ya tenemos descargados, y nos permite realizar cálculos sobre ellos como si se tratara de un único archivo, sin necesidad de fusionarlos físicamente.

Mínimos polígonos convexos

Procederemos ahora a generar los mínimos polígonos convexos para cada una de nuestras especies, los cuales son equivalentes al extent of occurrence (EOO) utilizado por la UICN. Estos polígonos los usaremos para recortar nuestro mosaico y así reducir el tiempo de cálculo, además de calcular el área de EOO posteriormente.

# 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()

Calculando la pérdida de bosque

Debido a que los datos de pérdida de bosque corresponden a un año en lugar de una proporción, lo siguiente que tenemos que hacer es encontrar la manera de transformarlos. Para esto, podemos agregar el ráster original, es decir, agrupar píxeles para reducir su resolución y utilizar las celdas dentro de cada nuevo píxel para generar una métrica de cambio de cobertura.

Para ello, primero definiremos un factor de agregación. En este ejemplo, llevaremos los píxeles desde 30 metros a 1 km:

# 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]

Posteriormente, convertiremos los mínimos polígonos convexos que generamos anteriormente a un objeto espacial (necesario para operar con algunos paquetes) y cargaremos, además, el contorno de América. Este último lo vamos a utilizar para realizar un segundo enmascaramiento, ya que algunos polígonos pueden tener áreas en el océano, las cuales, para este ejemplo y al tratarse de especies terrestres, nos generarían una sobreestimación.

# 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")

Por último, vamos a generar dos listas, cada una pensada para almacenar el resultado de las operaciones de agregación que realizaremos.

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

Luego, recorreremos la lista de especies y para cada una generaremos dos nuevos rásters agregados. Para treecover2000, asignaremos a la nueva celda el valor promedio de todos los píxeles contenidos dentro de la agregación, mientras que para lossyear, calcularemos la proporción de píxeles con valores mayores o iguales a 1 sobre el total de píxeles de la nueva celda agrupada. En el segundo caso, hacemos esto porque solo tenemos información sobre el año en que el bosque pasó a un estado sin bosque; es decir, podemos asumir los años como un carácter binario (0/1) y calcular cuántas celdas ya no tienen cobertura para un determinado período de tiempo. En otras palabras, calculamos los valores promedio de la cobertura de bosque en el año 2000 (FC00) y la pérdida de bosque durante el período 2001–2022 (FL) dentro del EOO de cada especie, así como la proporción de bosque remanente en la actualidad (FC22), de la siguiente manera:

\[ \text{FC}_{22} = \text{FC}_{00} - (\text{FC}_{00} \times \text{FL}) \] Además, definimos un umbral muy conservador del 95% y uno más relajado del 80% para evaluar la reducción en la calidad del hábitat. Específicamente, si la cobertura de bosque actual supera el 95% de la proporción disponible en el año 2000, consideramos que el efecto de la pérdida de bosque sobre la calidad del hábitat es bajo; si se encuentra por encima del 80% pero por debajo del 95%, lo consideramos un efecto intermedio (preocupante); y si la proporción es inferior al 80%, consideramos que el efecto es alto. Usamos el año 2000 como punto de comparación porque los paisajes naturalmente incluyen áreas no boscosas no antropizadas, por lo tanto, no tiene sentido esperar que la totalidad de un área esté cubierta por bosque.

# 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)

Tabla resúmen

Por último, solo tenemos que generar la tabla resumen con toda la información que hemos obtenido. Para ello, usaremos otra de las funciones personalizadas que cargamos al comienzo. Por defecto, esta función asume corrected = FALSE, es decir, calcula el área del EOO tal cual, sin tener en cuenta ningún contorno adicional. En este caso, como se mencionó anteriormente, dado que las especies son terrestres, usamos el contorno del continente para reducir el área de sobrepredicción. Además, el número de localidades se calculó con base en una distancia mínima de 18.5 km (aproximadamente 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

# Los datos en esta tabla pueden diferir sutilmente de los publicados en nuestro
# artículo, ya que originalmente este proceso se realizó usando una combinación
# de R y QGIS, además de un conjunto de datos adelgazado, donde los puntos que
# forman los polígonos mínimos convexos pueden variar ligeramente

Con esta tabla, ya podemos remitirnos a las guías de la UICN para realizar nuestra evaluación según el criterio B, ya que contamos con el EOO, el número de localidades y un umbral de cobertura de bosque como proxy de un declive estimado, inferido o proyectado en la calidad del hábitat.

Por ejemplo, para nuestra nueva especie Echinosaura embera, tenemos un EOO mayor a 20,000 km², por lo que no se cumple el criterio B1 y, por lo tanto, la especie no clasifica para ser considerada dentro de una categoría de amenaza. De igual forma, tiene más de 10 localidades conocidas y la cobertura de bosque en su EOO actual está por encima del 95% de la cobertura que había en el año 2000. Por ende, su categorización sería Preocupación Menor (LC).

Por otro lado, Echinosaura keyi, actualmente clasificada como Vulnerable (VU), tiene un área mayor a 5,000 km² pero menor a 20,000 km², por lo que cumple la primera condición de vulnerabilidad. Además, tiene 8 localidades conocidas, lo que la ubica de nuevo en la categoría Vulnerable Ahora bien, para que la especie pueda ser considerada como tal, es necesario que cumpla dos de las tres condiciones (a, b y c). Por número de localidades ya cumple una, por lo que tendríamos que justificar algún tipo de declive observado, estimado, inferido o proyectado en alguno de los parámetros que la guía nos indica. En este caso, utilizaremos la calidad del hábitat, pues para eso calculamos la pérdida de bosque. Si nos referimos a nuestra tabla, veremos que Echinosaura keyi todavía habita en un área que conserva, en promedio, más del 95% de la cobertura que tenía en el año 2000, por lo que consideramos el efecto de la deforestación bajo, y la segunda condición no aplicable. Por ende, en lugar de ser Vulnerable, la consideramos Casi Amenazada (NT), pues si las condiciones de favorabilidad del hábitat cambian en un futuro cercano, la especie podría pasar a estar amenazada.

Finalmente, si miramos Echinosaura panamensis, es una especie que actualmente se considera Preocupación Menor (LC); sin embargo, nuestros datos muestran que tiene un EOO menor a 20,000 km², menos de 10 localidades conocidas y un promedio de cobertura de bosque que se encuentra por debajo del 95% de la proporción de cobertura en el año 2000. Por ende, esta especie debería considerarse como Vulnerable (VU).