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(cacc)
library(ggsankey)
library(gt)

Nota

Los paquetes cageo, cacc 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)	-84.58961	10.427559	gbif	2024-01-02	4510305838
Agalychnis callidryas (Cope, 1862)	-84.64558	10.455464	gbif	2024-01-03	4510361946
Agalychnis callidryas (Cope, 1862)	-84.12152	10.451470	gbif	2024-01-04	4510345248
Agalychnis callidryas (Cope, 1862)	-84.12165	10.451480	gbif	2024-01-04	4510349158
Agalychnis callidryas (Cope, 1862)	-82.14613	9.254863	gbif	2024-01-02	4510122523
Agalychnis callidryas (Cope, 1862)	-84.21459	9.576750	gbif	2024-01-07	4510302449

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 <- ca_historic_worldclim(var = "bio",
                      res = 2.5,
                      path = tempdir(),
                      return_stack = TRUE)

Error in download.file(url, temp) : 
  download from 'https://biogeo.ucdavis.edu/data/worldclim/v2.1/base/wc2.1_2.5m_bio.zip' failed

Lo convertimos en un objeto del paquete terra para usarlo en el modelaje

Show the code

ca_bioc_t <- rast(ca_bioc)

Realizamos un mapa de las 19 variables bioclimáticas.

Show the code

plot(ca_bioc_t, 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 objetos 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_t, datos_entrenamiento)

Warning in .local(x, p, ...): 37 (6.94%) of the presence points have NA
predictor values

Loading required namespace: rJava

Graficams 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_t)

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_t, 1000)

Warning in backgroundSample(ca_bioc_t, 1000): generated random points = 0.741
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_t)

Visualizamos diferentes métricas de evaluación

Show the code

evaluacion_01

@stats
   np  na prevalence   auc   cor pcor   ODP
1 659 741      0.471 0.927 0.753    0 0.529

@thresholds
  max_kappa max_spec_sens no_omission equal_prevalence equal_sens_spec
1     0.377         0.377       0.004            0.472           0.401

@tr_stats
    treshold kappa  CCR  TPR TNR FPR  FNR  PPP  NPP  MCR  OR
1          0     0 0.47    1   0   1    0 0.47  NaN 0.53 NaN
2          0     0 0.47    1   0   1    0 0.47  NaN 0.53 NaN
3          0     0 0.47    1   0   1    0 0.47    1 0.53 Inf
4        ...   ...  ...  ... ... ...  ...  ...  ...  ... ...
949     0.99  0.01 0.53 0.01   1   0 0.99    1 0.53 0.47 Inf
950        1     0 0.53    0   1   0    1  NaN 0.53 0.47 NaN
951        1     0 0.53    0   1   0    1  NaN 0.53 0.47 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, cacc 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_inm_cm4_8_ssp126_2041_2060 <- ca_future_worldclim(var = "bioc",
                    res = 2.5,
                    gcm = "INM-CM4-8",
                    ssp = "ssp585",
                    interval = "2041-2060",
                    path = tempdir(),
                    return_stack = TRUE)

Al raster obtenido lo convertimos en un objeto terra para usarlo en el modelaje para analizar el impacto del cambio climático.

Show the code

ca_bio_2041_2060_t <- rast(ca_inm_cm4_8_ssp126_2041_2060)

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_t) <- names(ca_bioc_t)

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

Show the code

mapa_maxent_2041_2060 <- predict(modelo_maxent, ca_bio_2041_2060_t)

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)

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ínola2023,
  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 = {2023-12-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. 2023. “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.” December 8, 2023. https://mspinola-ciencia-de-datos.netlify.app/posts/2023-12-06-actividad_2020_2022/actividad_2020_2022.html.