Modelaje de la distribución actual de la rana verde de ojos rojos

Este es un flujo de trabajo de como realizar un análisis del impacto del cambio climático sobre la biodiversidad en América Central. En este caso utilizamos como ejemplo la rana verde de ojos rojos (Agalychnis callidryas).

Para poder realizar el análisis es necesario instalar y cargar los siguientes paquetes de R (R Core Team 2023).

Show the code

library(tidyverse)
library(here)
library(spocc)
library(sf)
library(stars)
library(terra)
library(tidyterra)
library(predicts)
library(cageo)
library(paisaje)
library(ggsankey)
library(gt)

Nota

Los paquetes cageo, paisaje y ggsankey se deben bajar desde el enlace que se provee al pinchar en los respectivos nombres.

1. Obtención de los límites del área de estudio

Para obtener los límites de América Central usamos el mapa ca_outline_d que se incluye en el paquete de R, cageo (Spínola 2023b)

Nota

El paquete de R, cageo contiene datos espaciales para America Central. El paquete se desarrolló para facilitar este flujo de trabajo.

Show the code

ca <- ca_outline_d

Show the code

ggplot() +
  theme_minimal() +
  geom_sf(data = ca)

2. Obtenención de los registros de la especie

Para obtener los registros de la especie en el portal de datos GBIF usamos el paquete spocc (Owens, Barve, y Chamberlain 2023)

Show the code

ac_gbif <- occ(query = "Agalychnis callidryas", from = "gbif", has_coords = TRUE, limit = 5000, geometry = st_bbox(ca))

Los registros obtenidos para la región de América Central los convertimos en una estructura de datos de 2 dimensiones (filas y columnas).

Show the code

ac_gbif_df <- occ2df(ac_gbif) |>
  drop_na() |>
  mutate_if(is.character, as.factor)

Show the code

head(ac_gbif_df) |> gt() |> tab_options(table.font.size = 12)

name	longitude	latitude	prov	date	key
Agalychnis callidryas (Cope, 1862)	-83.50432	10.53795	gbif	2025-01-01	5006812286
Agalychnis callidryas (Cope, 1862)	-84.06817	10.45032	gbif	2025-01-03	5006976124
Agalychnis callidryas (Cope, 1862)	-84.65591	10.45616	gbif	2025-01-04	5007621545
Agalychnis callidryas (Cope, 1862)	-84.66306	10.46189	gbif	2025-01-04	5007625747
Agalychnis callidryas (Cope, 1862)	-84.06824	10.45039	gbif	2025-01-03	5007746235
Agalychnis callidryas (Cope, 1862)	-83.50235	10.54248	gbif	2025-01-01	5007782176

Luego convertimos los registros en un objeto espacial, donde solo nos quedamos con los registros que se incluyen en la región continental de América Central

Show the code

ac_gbif_sp <- st_as_sf(ac_gbif_df, coords = c("longitude", "latitude"), crs = 4326) |>
  st_intersection(ca)

Se remueven los registros que puedan estar duplicados

Show the code

ac_gbif_sp <- distinct(ac_gbif_sp)

Elaboramos un mapa que muestre los registros en la región de América Central

Show the code

ggplot() +
  geom_sf(data = ca) +
  geom_sf(data = ac_gbif_sp)

3. Preparación de las variables ambientales

Las variables ambientales a utilizar en este caso son las 19 variables bioclimáticas elaboradas por WorldClim. Para obtener las variables para la región de América Central utilizamos el paquete de R, cacc (Spínola 2023a)

Show the code

ca_bioc <- get_worldclim_historic(var = "bio",
                      res = 2.5,
                      aoi = ca)

Realizamos un mapa de las 19 variables bioclimáticas.

Show the code

plot(ca_bioc, maxnl = 19)

4. Modelación

Ajuste del modelo

Para modelar la distribución de la rana verde de ojos rojos usamo el algoritmo Maxent del paquete de R, predicts (Hijmans 2023)

Primero transformamos el objeto ac_gbif_sp en un objeto terra

Luego dividimos los registros de ocurrencia en 5 subconjuntos y se dividen los datos en un subconjunto de entrenamiento y se mantiene un 20% para poner a prueba el modelo

Show the code

sc <- folds(ac_gbif_sp_t, k=5)
datos_prueba <- ac_gbif_sp_t[sc == 1, ]
datos_entrenamiento <- ac_gbif_sp_t[sc != 1, ]

Ajustamos el modelo de maxent

Show the code

modelo_maxent <- MaxEnt(ca_bioc, datos_entrenamiento)

Loading required namespace: rJava

Graficamos la importancia de las variables.

Show the code

plot(modelo_maxent)

Predicción

Predecimos la cálidad del hábitat para la especie para la región de América Central

Show the code

mapa_maxent <- predict(modelo_maxent, ca_bioc)

Elaboramos un mapa con la calidad del hábitat

Show the code

ggplot() +
  theme_minimal() +
  geom_spatraster(data = mapa_maxent) +
  scale_fill_whitebox_c(name = "Calidad \ndel hábitat",
                       palette = "Bl_Yl_Rd",
                       na.value = NA,
                       direction = 1)

Evaluación

Generamos datos de “background”, en este caso 1000 puntos

Show the code

bg <- backgroundSample(ca_bioc, 1000)

Warning in backgroundSample(ca_bioc, 1000): generated random points = 0.787
times requested number

Para evaluar el modelo usamos los puntos de “background” con los datos de prueba que habíamos asignado previamente

Show the code

evaluacion_01 <- pa_evaluate(modelo_maxent, p = datos_prueba, a = bg, x = ca_bioc)

Visualizamos diferentes métricas de evaluación

Show the code

evaluacion_01

@stats
   np  na prevalence   auc   cor pcor   ODP
1 979 787      0.554 0.966 0.851    0 0.446

@thresholds
  max_kappa max_spec_sens no_omission equal_prevalence equal_sens_spec
1     0.513         0.513       0.025            0.558           0.544

@tr_stats
     treshold kappa  CCR  TPR TNR FPR  FNR  PPP  NPP  MCR  OR
1           0     0 0.55    1   0   1    0 0.55  NaN 0.45 NaN
2           0     0 0.55    1   0   1    0 0.55    1 0.45 Inf
3           0     0 0.55    1   0   1    0 0.55    1 0.45 Inf
4         ...   ...  ...  ... ... ...  ...  ...  ...  ... ...
1044     0.99  0.01 0.45 0.01   1   0 0.99    1 0.45 0.55 Inf
1045     0.99  0.01 0.45 0.01   1   0 0.99    1 0.45 0.55 Inf
1046     0.99     0 0.45    0   1   0    1  NaN 0.45 0.55 NaN

Graficamos la evaluación del modelo.

Show the code

plot(evaluacion_01, "ROC")

Obtenemos el threshold para este modelo para usarlo en la evaluación del impacto del cambio climático.

5. Evaluación del impacto del cambio climático

Para evaluar el impacto del cambio climático sobre el hábitat de la rana verde de ojos rojos utilizamos un escenario de cambio climático.

En este caso usaremos el escenario INM-CM4-8 para SSP5-8.5 y el intervalo de tiempo 2041 - 2060.

Nota

El paquete de R, paisaje cuenta con funciones para bajar datos climatológicos para América Central. El paquete se desarrolló para facilitar este flujo de trabajo.

Show the code

ca_bio_2041_2060 <- get_worldclim_future(var = "bioc",
                    res = "2.5m",
                    gcm = "INM-CM4-8",
                    scenario = "585",
                    time_range = "2041-2060",
                    aoi = ca)

[1] "Download URL: https://geodata.ucdavis.edu/cmip6/2.5m/INM-CM4-8/ssp585/wc2.1_2.5m_bioc_INM-CM4-8_ssp585_2041-2060.tif"

Renombramos las 19 variables bioclimáticas con el mismo nombre de las variables usadas en el modelo para el tiempo presente.

Show the code

names(ca_bio_2041_2060) <- names(ca_bioc)

Realizamos la predicción para el escenario climático escogido.

Show the code

mapa_maxent_2041_2060 <- predict(modelo_maxent, ca_bio_2041_2060)

Elaboramos un mapa de la predicción.

Show the code

ggplot() +
  theme_minimal() +
  geom_spatraster(data = mapa_maxent_2041_2060) +
  scale_fill_whitebox_c(name = "Calidad \ndel hábitat",
                       palette = "Bl_Yl_Rd",
                       na.value = NA,
                       direction = 1)

Combinamos ambos mapas para poder compararlos.

Show the code

ambos <- c(mapa_maxent, mapa_maxent_2041_2060)

Le asignamos nombres a las capas.

Show the code

names(ambos) <- c("Actual", "Futuro")

Como se puede ver, se registran cambios en la calidad del hábitat para el escenario climático utilizado.

Show the code

ggplot() +
  theme_minimal() +
  geom_spatraster(data = ambos) +
  facet_wrap(~ lyr) +
  scale_fill_whitebox_c(name = "Calidad \ndel hábitat",
                       palette = "Bl_Yl_Rd",
                       na.value = NA,
                       direction = 1)

Una manera de cuantificar los cambios en un escenario futuro es transformar la calidad del hábitat, por ejemplo, en 5 categorías.

Lo hacemos para el mapa de calidad de hábitat actual.

Show the code

mapa_maxent_c <- classify(mapa_maxent, rclmat, include.lowest = TRUE)

Show the code

levels(mapa_maxent_c) <- clases

Lo hacemos para el mapa de calidad de hábitat en el futuro.

Show the code

mapa_maxent_2041_2060_c <- classify(mapa_maxent_2041_2060, rclmat, include.lowest = TRUE)

Show the code

levels(mapa_maxent_2041_2060_c) <- clases

Combinamos ambos mapas (rasters)

Show the code

ambos_c <- c(mapa_maxent_c, mapa_maxent_2041_2060_c)

Show the code

names(ambos_c) <- c("Actual", "Futuro")

Elaboramos un gráfico con ambos mapas para hacer más fácil la comparación.

Show the code

ggplot() +
  theme_minimal() +
  geom_spatraster(data = ambos_c, color = "black") +
  facet_wrap(~ lyr) +
  scale_fill_whitebox_d(name = "Calidad \ndel hábitat",
                       palette = "Bl_Yl_Rd",
                       na.value = NA,
                       direction = 1)

Luego procesamos ambos mapas para poder contabilizar las celdas (pixeles) que corresponden a cada categoría de calidad del hábitat en el mapa actual y en el mapa futuro.

Show the code

r_df <- as.data.frame(ambos_c) |>
  drop_na()

Show the code

r_df_sankey <- r_df %>%
  make_long(Actual, Futuro)

Show the code

r_df_sankey$node <- factor(r_df_sankey$node, levels = c("Muy bajo", "Bajo", "Medio", "Alto", "Muy alto"))

Elaboramos un sankey plot para visualizar la transformación de las celdas (pixeles) entre diferentes categorías de calidad del hábitat entre el mapa actual y el mapa futuro. .

Show the code

ggplot(r_df_sankey, aes(x = x, 
               next_x = next_x, 
               node = node, 
               next_node = next_node,
               fill = factor(node),
               label = node)) +
  theme_sankey(base_size = 18) +
  labs(x = NULL) +
  geom_sankey(flow.alpha = 0.5
                      , node.color = "black"
                      ,show.legend = FALSE) +
  geom_sankey_label(size = 3, color = "black", fill= "white", hjust = -0.5) +
  scale_fill_brewer(palette = "Set1", direction=-1)

Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
ℹ Please use the `linewidth` argument instead.
ℹ The deprecated feature was likely used in the ggsankey package.
  Please report the issue at <https://github.com/davidsjoberg/ggsankey/issues>.

Referencias

Hijmans, Robert J. 2023. predicts: Spatial Prediction Tools. https://CRAN.R-project.org/package=predicts.

Owens, Hannah, Vijay Barve, y Scott Chamberlain. 2023. spocc: Interface to Species Occurrence Data Sources. https://CRAN.R-project.org/package=spocc.

R Core Team. 2023. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Spínola, Manuel. 2023a. cacc: Central America climate change. https://github.com/ManuelSpinola/cacc.

———. 2023b. cageo: Geospatial data of Central America. https://github.com/ManuelSpinola/cageo.

Cómo citar

BibTeX

@online{spínola2024,
  author = {Spínola, Manuel and Sáenz, Joel and Retamosa, Mónica},
  title = {Flujo de trabajo para la evaluación del cambio climático en
    la biodiversidad de América Central mediante el uso de datos
    abiertos y software libre},
  date = {2024-01-08},
  url = {https://mspinola-ciencia-de-datos.netlify.app/posts/2023-12-06-actividad_2020_2022/actividad_2020_2022.html},
  langid = {es}
}

Por favor, cita este trabajo como:

Spínola, Manuel, Joel Sáenz, and Mónica Retamosa. 2024. “Flujo de trabajo para la evaluación del cambio climático en la biodiversidad de América Central mediante el uso de datos abiertos y software libre.” January 8, 2024. https://mspinola-ciencia-de-datos.netlify.app/posts/2023-12-06-actividad_2020_2022/actividad_2020_2022.html.