Macroecología
Instituto de Ecología, A.C. - INECOL
dismo (Hijmans & Elith 2011), código de Vijay Barve (https://vijaybarve.wordpress.com) y Luis Verde-Arregoitia (https://luisdva.github.io/)Ya vimos de dónde vienen los datos primarios de biodiversidad (e.g. datos de presencia; almacenados en museos, colecciones científicas). Muchos de esos datos ahora están disponibles en internet!!!
Cargar los paquetes (TODOS!) necesarios para este ejercicio (estos ya deben estar instalados en su computadora)
library(sf)
library(terra)
library(maps)
library(dismo)
library(rnaturalearth)
library(tidyterra)
library(ggplot2)
library(rgbif)
library(dplyr)
library(alphahull)
library(rangeBuilder)
Obtener datos de presencia de una especie cualquiera a partir de la base de datos GBIF (data.gbif.org) directamente desde R! NOTA: Necesitan estar conectados a internet. Este proceso puede tardar, dependiendo de la especie, si esta tiene (o no) muchos registros
sp_datos <- occ_search(scientificName = "Musonycteris harrisoni", limit = 500)
NOTA: Uds podrían colocar el nombre del género y la especie deseada.
El objeto sp_datos es una lista con datos sobre los resultados obtenidos en GBIF, para trabjar únicamente con la tabla de registros hay que seleccionar el objeto data dentro del mismo
sp_data <- sp_datos$data
Ver las dimensiones del objeto generado (i.e. cuántas líneas y cuántas columnas tiene)
dim(sp_data)
Checar el nombre de las columnas (para después buscar únicamente las de posición geográfica: lat/long)
names(sp_data)[1:30]
Crear otro objeto a partir del anterior para quedarse únicamente con long/lat. Quedarse únicamente con los puntos/registros individuales (i.e., excluir duplicados). Transformarlo en un objeto espacial
sp_p1<-sp_data %>%
select(decimalLongitude,decimalLatitude,species) %>%
mutate(lat=decimalLatitude,lon=decimalLongitude) %>%
distinct() %>%
na.omit() %>%
sf::st_as_sf(coords = c('decimalLongitude','decimalLatitude'),crs="EPSG: 4326")
NOTA: el nombre de la variable puede ser diferente (.e.g “LATITUDE”, “Latidude”, “lat”, etc. Siempre hay que checar antes)
Graficar (poner en un mapa) los puntos de presencia de nuestra especie
ggplot()+
geom_sf(data=sp_p1, col="blue",pch=19)
Agregar el mapa del mundo. Para eso, necesitamos cargar el mapa (que viene con el paquete “rnaturalearth”)
wrld <- ne_countries(scale = "small",returnclass = "sf")
ggplot()+
geom_sf(data=wrld)+geom_sf(data=sp_p1,col="blue",pch=19,size=1)+coord_sf(expand = F) +
labs(x="Longitud decimal ",
y="Latitud decimal",
title=expression(paste("Puntos reportados ", italic("Musonycteris harrisoni"))))+
theme(plot.title = element_text(hjust = 0.5))
¿Hay algún dato extraño? Los puntos/registros necesitan ser “curados” (limpiados)
Eliminar los puntos con mala georeferencia (en este caso, puntos obvios en el “viejo mundo”)
sp_p1<-sp_data%>%
select(decimalLongitude,decimalLatitude,species)%>%
mutate(lat=decimalLatitude,lon=decimalLongitude)%>%
distinct() %>% na.omit() %>%
sf::st_as_sf(coords = c('decimalLongitude','decimalLatitude'),crs="EPSG: 4326")%>%
filter(lat> 0.5) %>% filter(lat< 22)
Ahora sí, mapeamos de nuevo pero sólamente en la región de interés (México)
mex_map <- filter(wrld,name=="Mexico")
ggplot()+geom_sf(data=mex_map)+
geom_sf(data=sp_p1,col="blue",pch=19,size=1)+coord_sf(expand = F)
Y ¿Cómo eliminamos los registros que están en el mar?
sp_p1<-sp_data%>%
select(decimalLongitude,decimalLatitude,species)%>%
mutate(lat=decimalLatitude,lon=decimalLongitude)%>%
distinct() %>%
na.omit() %>%
sf::st_as_sf(coords = c('decimalLongitude','decimalLatitude'),crs="EPSG: 4326")%>%
filter(lat> 0.5) %>%
filter(lat< 22)%>%
filter(lon> -105.56611)
En esta parte, se usaran funciones para generar mapas “simples” basados en la geometría (e.g. distancias/ángulos entre puntos, etc.). NO consideran variables ambientales!
Este modelo crea un polígono convexo alrededor de los puntos de presencia.
sp1_mcp <- st_convex_hull(st_union(sp_p1))
¿Cómo se ve?
ggplot()+
geom_sf(data=mex_map)+
geom_sf(data=sp1_mcp,
fill="blue")
Usamos el paquete alphahull
NOTA: Esta función solo acepta tablas como entrada
sp_p2<-as.data.frame(st_coordinates(sp_p1))
sp1_alphahull <- ahull(sp_p2, alpha = 6)
Error: shull: duplicate points found
Falla porque encuentra puntos duplicados o, como en este caso, puntos en una línea recta (i.e, mismo X y/o mismo Y).
¿Cómo podemos identificar y solucionar este error?
sp_p3<-sp_p2 %>% select(X, Y)%>%
mutate(X = ifelse(duplicated(sp_p2$X), X + rnorm(1, mean = 0, sd = 0.0001), X))%>%
mutate(Y = ifelse(duplicated(sp_p2$Y), Y + rnorm(1, mean = 0, sd = 0.0001), Y))
Ahora si, podemos crear el Alpha Hull con un valor de alpha escogido (por la razón que crean relevante)
sp1_alphahull<- ahull(sp_p3, alpha = 6)
Para observar el alpha hull, necesitamos que el objeto sea de tipo espacial del paquete sf. Para eso usaremos una función independiente, disponible en su carpeta de trabajo
source(file ="ejercicios_datos/ah2sf.R")
sp1_alphahull.poly <- ah2sf(sp1_alphahull)
¿Cómo se ve?
ggplot()+
geom_sf(data=mex_map)+geom_sf(data=sp1_alphahull.poly,fill="blue")
Usamos el paquete rangeBuilder, el cual crea un polígono alpha hull con un valor de alpha “óptimo” basado en la distribución espacial de los puntos.
sp_range <- getDynamicAlphaHull(
sp_p3, #Tabla de puntos/registros de la especie
coordHeaders = c("decimalLongitude", "decimalLatitude"),# x y y
fraction = 0.95, # la fracción mínima de registros que debe incluir el polígono
partCount = 2, # el máximo de polígonos disyuntos permitidos
initialAlpha = 1, # Alpha inicial
alphaIncrement = 0.5,
alphaCap = 1000,
clipToCoast = "terrestrial" # solo la parte terrestre del polígono se mantendrá (se cortan las partes no-terrestres/acuáticas con base en un mapa descargado de naturalearth).
)
alpha <- sp_range[[2]]
alpha
Convertir el polígono alpha a un objeto sf
sp1_dynalpha <- st_make_valid(st_as_sf(sp_range[[1]]))
¿Cómo podemos visualizarlo?
ggplot()+ geom_sf(data=mex_map)+
geom_sf(data=sp1_dynalpha, fill="blue")
ggplot()+
geom_sf(data=mex_map)+ geom_sf(data=sp1_mcp,fill="red",alpha=0.1) +
geom_sf(data=sp1_alphahull.poly,fill="yellow",alpha=0.5) +
geom_sf(data=sp1_dynalpha, fill="blue",alpha=0.5)
Finalmente, podemos salvar esos polígonos como shapefiles, para usarlos en otros software (e.g. ArcGIS) y eventualmente juntar los de varias especies para otros análisis (ejercicio siguiente)
st_write(sp1_mcp, "sp1_min_convex.shp")
st_write(sp1_alphahull.poly, "sp1_alphahull.shp")
st_write(sp1_dynalpha, "sp1_dyn_alphahull.shp")
Ahora, vamos a incluir datos ambientales/climáticos
Bajar los datos de clima del site www.worldclim.org (esto puede tardarse un poco, un par de minutos)
world_bioclim<-getData("worldclim", var = "bio", res = 10)
world_bioclim <- getData(name='worldclim',res=10,var='bio',path="/Users/fabricio/Library/CloudStorage/Dropbox/Macroecologia_2023/",download = F)
Estas variables están todas “juntas”. Vamos a separarlas y dejar un objeto a parte que se puede usar para agarrar cada una de las variables separadas una por una
world_bioclim_inds<-unstack(world_bioclim)
BIO1 = Annual Mean Temperature
BIO2 = Mean Diurnal Range (Mean of monthly (max temp - min temp))
BIO3 = Isothermality (BIO2/BIO7) (* 100)
BIO4 = Temperature Seasonality (standard deviation *100)
BIO5 = Max Temperature of Warmest Month
BIO6 = Min Temperature of Coldest Month
BIO7 = Temperature Annual Range (BIO5-BIO6)
BIO8 = Mean Temperature of Wettest Quarter
BIO9 = Mean Temperature of Driest Quarter
BIO10 = Mean Temperature of Warmest Quarter
BIO11 = Mean Temperature of Coldest Quarter
BIO12 = Annual Precipitation
BIO13 = Precipitation of Wettest Month
BIO14 = Precipitation of Driest Month
BIO15 = Precipitation Seasonality (Coefficient of Variation)
BIO16 = Precipitation of Wettest Quarter
BIO17 = Precipitation of Driest Quarter
BIO18 = Precipitation of Warmest Quarter
BIO19 = Precipitation of Coldest Quarter
Graficar una de las variables
plot(world_bioclim_inds[[1]], zlim=c(-270,320), axes = F, box=F, col=rev(heat.colors(25)), ext = matrix(c(180,90,-180,-90), nrow = 2))
O un plot más “fancy”
bio1 <- as(world_bioclim_inds[[1]],"SpatRaster")
ggplot()+
geom_spatraster(data=bio1,aes(fill=bio1)) + scale_fill_gradientn(colours=rev(c("brown","red","yellow","darkgreen","green")))
Grafiquen algunas de esas variables para América y para México. Para eso, hay que modificar el “ext=matrix(c(XX,XX,XX,XX)))” de la función anterior, colocando las coordenadas correctas para América y para México (Por separado, o sea, un mapa para cada una)
##Ecological Niche Model (Bioclim)
Vamos a quedarnos únicamente con dos variables
bio1<- world_bioclim_inds[[1]]
bio12<- world_bioclim_inds[[12]]
env.variables<- stack (bio1, bio12)
Establecer el “extent” de América del Sur y cortar las variables para ese extent
mex.extent <- c(-117,-87,15,33)
mex.env <- crop(env.variables,mex.extent)
Generar el modelo Bioclim (Ecological Niche Model: ENM)
# Convertir los puntos de la especie en un objeto tipo `matrix`
sp1_points_mat <- as.matrix(sp_p3)
# Correr el ENM (Bioclim)
sp1_bioclim <- bioclim(mex.env, sp1_points_mat)
plot(sp1_bioclim, a=1, b=2, p=0.85)
# Generar una predicción del área adecuada ("suitable") basada en el ENM
sp1.map<- predict (mex.env, sp1_bioclim)
# Plotar o resultado da predição baseada no ENM
plot (sp1.map)
# Evaluar el modelo ENM
group <- kfold(sp_p3, 5)
pres_train <- sp_p3[group != 1, ]
pres_test <- sp_p3[group == 1, ]
backg <- randomPoints(mex.env, n=1000, ext=mex.extent, extf = 1.25)
colnames(backg) = c('lon', 'lat')
group <- kfold(backg, 5)
backg_train <- backg[group != 1, ]
backg_test <- backg[group == 1, ]
e <- evaluate(pres_test, backg_test, sp1_bioclim, mex.env)
# Establecer un "threshold" para cortar la predicción y generar un mapa binario
threshold <- e@t[which.max(e@TPR + e@TNR)]
plot(sp1.map > threshold, main='presence/absence')
# Agregar los puntos "reales" de la especie
points(sp1_points_mat[,1],sp1_points_mat[,2],cex=0.2)
¿Cómo se ve?
¿Hay regiones en donde la especie fue “predicha” (en realidad, las condiciones) pero en las que sabemos que no está?
NOTA: Salven el “workspace” de R, es decir, todo lo que se hizo en este ejercicio debe quedar salvado como un archivo “.Rdata”. Estos objetos serán usados en ejercicios posteriores. NO CIERREN R SIN HABER GUARDADO ESTE ARCHIVO!!!