Muchas aplicaciones de sistemas de información geográfica en macroecología requieren de análisis estadísticos. La plataforma por excelencia para realizarlos es R, aún más con las excelentes capacidades de R por sí solo como SIG. En este tutorial veremos cómo:

  1. Importar datos geográficos: puntos, líneas, polígonos y raster

  2. Transformar entre tipos de datos

  3. Reproyectar datos

  4. Hacer un análisis de regresión con los datos geográficos

    • Interpretar el análisis

Importar datos

Al igual que en el tutorial de GBIF utilizaremos el paquete terra para importar y manipular todos los datos geográficos. De los datos usados, los puntos, líneas y polígonos pueden ser manejados como vectoriales. Las extensiones de archivo más comunes para estos son gpkg y shp. Los puntos, por carecer de más características geométricas son más sencillos y se les suele también encontrar como csv.

Comenzaremos cargando el paquete terra, y la capa vectorial de Yucatán, Campeche y Quitana Roo:

library(terra); library(tidyr)
py <- vect("Datos/Peninsula.gpkg")

el objeto py es el que contiene la capa, y podemos ver su tabla de atributos imprimiendo el objeto como un data.frame (usando la sintaxis de tidyr):

py |> as.data.frame()
##      GID_1 GID_0 COUNTRY       NAME_1 VARNAME_1 NL_NAME_1 TYPE_1 ENGTYPE_1 CC_1
## 1  MEX.4_1   MEX  México     Campeche        NA        NA Estado     State   NA
## 2 MEX.23_1   MEX  México Quintana Roo        NA        NA Estado     State   NA
## 3 MEX.31_1   MEX  México      Yucatán        NA        NA Estado     State   NA
##   HASC_1  ISO_1
## 1  MX.CM MX-CAM
## 2  MX.QR MX-ROO
## 3  MX.YU     NA

Para ver la proyección geográfica de la capa:

crs(py, describe = T)
##     name authority code  area             extent
## 1 WGS 84      EPSG 4326 World -180, 180, 90, -90

Para importar capas de líneas se usa la misma función vect. Para importar datos de puntos, podemos usar read.csv ó vect dependiendo de la extensión del archivo. Si los datos los obtuvimos de GBIF serán csv. Si los datos de GBIF lo hemos limpiado con anterioridad, es posible que los hayamos guarado como gpkg ó shp:

cd.csv <- read.csv("Datos/Cachryx.csv")
cd.gpkg <- vect("Datos/Cachryx.csv")

Al verificar la proyección de estos datos, podemos ver que está vacía:

crs(cd.gpkg, describe = T)
##      name authority code area         extent
## 1 unknown      <NA> <NA> <NA> NA, NA, NA, NA

y podemos expecificarla utilizando la proyección de la capa vectorial de la península, si de antemano sabemos que tienen la misma:

crs(cd.gpkg) <- crs(py)
crs(cd.gpkg, describe = T)
##     name authority code  area             extent
## 1 WGS 84      EPSG 4326 World -180, 180, 90, -90

Para especificar la proyección de puntos importados desde un csv, necesitamos transformarlo a un objeto reconocible por terra:

cd.v <- vect(cd.csv)
crs(cd.v) <- crs(py)
crs(cd.v, describe = T)
##     name authority code  area             extent
## 1 WGS 84      EPSG 4326 World -180, 180, 90, -90

Reproyectar datos

El sistema de coordenadas más común es WGS84, que en la lista EPSG tiene el número 4326. WGS84 tiene coordenadas polares, medidas en grados, lo que quiere decir que son coordenadas geográficas y se deben utilizar para ubicar objetos en la superficie terrestre. Cuando queremos medir distancias ó áreas, necesitamos coordenadas proyectadas, que representen la superficie de la tierra con la menor cantidad de distorsiones en una superficie plana como la pantalla de la computadora. Los sistemas de coordenadas están diseñados para regiones específicas, y las podemos encontrar utilizando el nombre de la región y los códigos EPSG. En QGIS podemos buscar los sistemas de coordenadas haciendo click en el botón “EPSG” de al esquina inferior derecha:

El sistema por defecto es precisamente EPSG, con cobertura global:

En esta ventana podemos buscar sistemas de coordenadas propias de una región geográfica. Para México por ejemplo, está ITRF2008, EPSG:6363.

Además, hay proyecciones para zonas específicas, como la zona 16N de ITRF2008 (EPSG:6371):

Una vez identificadas las proyecciones adecuadas para nuestros datos, haremos el cambio de proyección a ITRF2008-16N (EPSG:6371):

cd.v.it <- project(cd.v, "EPSG:6371")
py.it <- project(py, "EPSG:6371")
plot(py.it); points(cd.v.it)

El sistema de referencia de las coordenadas es importante cuando necesitamos calcular áreas ó distancias. En macroecología y biogeografía entonces, es de muy alta importancia poder estimar adecuadamente las áreas.

Análisis de regresión con datos geográficos

Importaremos aquí las distribuciones de mamíferos de PHYLACINE, para hacer un análisis de las variables ambientales a que está asociada la riqueza. Comenzaremos importando las distribuciones con la función rast de terra. A esta función le podemos dar uno ó cientos de nombres de archivos para generar una capa raster con múltiples bandas. La lista de archivos podemos crearla con list.files:

f <- list.files("../Current/Current/", pattern = ".tif", full.names = T)
length(f)
## [1] 5831

Como podemos ver, hay 5831 archivos raster con extensión tif, de ahí que a list.files primero le indicamos en qué carpeta buscar (../Current/Current/) y después la cadena de texto que debe identificar en los archivos, que es la extensión. Una vez generada la lista de archivos, la importamos con rast:

r <- rast(f)

Debido a que ahora el objeto r tiene 5831 capas, es poco factible graficar todas las capas, aunque podemos mostrar un par de las distribuciones que contiene:

plot(r[[400]])

Para ver la riqueza global de mamíferos, podemos hacer la suma de todas las capas, de modo que obtengamos el número de especies por píxel:

riq <- sum(r)
plot(riq)

Como podemos ver, estos datos contienen las especies marinas y terrestres. Podríamos extraer las spp. terrestres usando como máscara una capa vectorial de las regiones terrestres. De paso, la proyectaremos a el sistema de coordenadas proyectadas global para tratar de preservar en la mayor medida posible las áreas de los polos:

library(maptools)
data("wrld_simpl")
m <- vect(wrld_simpl)
m <- project(m, "EPSG:4326")
plot(m)

riq <- project(riq, "EPSG:4326")
riq.ter <- mask(riq, m)
plot(riq.ter)

Ahora, importaremos las variables bioclimáticas de Chelsa para usar como variables independientes de nuestro modelo de regresión. De Chelsa usaremos las variables bioclimáticas 1, 7 y 12 (temperatura anual promedio, rango anual de temperatura y precipitación total). Estas variables también están en formato tif y las podemos importar y proyectar con terra:

bio <- rast(paste0("../Current/bio",c(1, 7, 12), ".tif"))
bio.res <- resample(bio, riq, method = "lanczos")
bio.res <- mask(bio.res, m)
plot(bio.res)

Para hacer un análisis de regresión usando la riqueza como variable dependiente de un conjunto de variables climáticas debemos transformar la capa raster de riqueza a tabla. Antes de ello, necesitamos cambiar el sistema de coordenadas de las capas de riqueza y climáticas a una proyección que preserve mejor las áreas que el sistema de coordenadas geográficas WGS84:

m.88 <- project(m, "EPSG:8857")
riq.88 <- project(riq.ter, "EPSG:8857")
bio.88 <- project(bio.res, "EPSG:8857")

riq.88 <- mask(riq.88, m.88)
bio.88 <- mask(bio.88, m.88)

Transformando de raster a tabla:

todos <- c(bio.88, riq.88)
riq.df <- as.data.frame(todos, xy = T)
riq.df$sum <- round(riq.df$sum, 0)
knitr::kable(head(riq.df))
x y bio1 bio7 bio12 sum
69 -12085265 7313734 -5.073144 -243.8959 289.2390 32
70 -12010090 7313734 -6.002041 -241.0531 387.0372 34
71 -11934914 7313734 -6.714367 -237.8318 591.0377 34
72 -11859738 7313734 -5.802148 -241.6073 608.9218 35
77 -11483858 7313734 -3.872388 -243.9792 190.6438 43
78 -11408682 7313734 -4.499725 -239.2193 344.9317 44

Ahora, podemos hacer un análisis exploratorio, viendo por ejemplo los gráficos de dispersión de la riqueza como función de las variables bioclimáticas:

library(ggplot2)

ggplot(riq.df) + geom_point(aes(x = bio1, y = sum)) + 
  geom_smooth(aes(x = bio1 , y = sum), method = "gam") +
  labs(y = "Riqueza")

ggplot(riq.df) + geom_point(aes(x = bio7, y = sum)) + 
  geom_smooth(aes(x = bio7 , y = sum), method = "gam") + 
  labs(y = "Riqueza")

ggplot(riq.df) + geom_point(aes(x = bio12, y = sum)) + 
  geom_smooth(aes(x = bio12 , y = sum), method = "gam")+ 
  labs(y = "Riqueza")

Para analizar la riqueza utilizaremos un modelo lineal generalizado Poisson, utilizando los siguientes términos:

  1. Bio1: lineal
  2. Bio7: lineal y cuadrático
  3. Bio12: lineal y cuadrático
m1 <- glm(sum ~ bio1 + bio7 + I(bio7^2) + 
            bio12 + I(bio12^2), 
          data = riq.df, family = "poisson")
summary(m1)
## 
## Call:
## glm(formula = sum ~ bio1 + bio7 + I(bio7^2) + bio12 + I(bio12^2), 
##     family = "poisson", data = riq.df)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.556e+00  2.897e-01  -12.27   <2e-16 ***
## bio1         1.108e-02  1.194e-04   92.77   <2e-16 ***
## bio7        -5.345e-02  2.336e-03  -22.88   <2e-16 ***
## I(bio7^2)   -9.740e-05  4.703e-06  -20.71   <2e-16 ***
## bio12        7.041e-04  2.464e-06  285.71   <2e-16 ***
## I(bio12^2)  -1.383e-07  6.708e-10 -206.24   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 665170  on 24967  degrees of freedom
## Residual deviance: 326630  on 24962  degrees of freedom
## AIC: 474850
## 
## Number of Fisher Scoring iterations: 4

De acuerdo con este análisis todas las variables afectan significativamente la riqueza. Ahora veremos cómo representa este modelo los patrones de riqueza:

p <- predict(m1)
riq.df$pred.riq <- exp(p)
riq.pred <- rast(riq.df)
par(mfrow = c(2, 1))
plot(riq.pred[["pred.riq"]], main = "Modelo")
plot(riq.pred[["sum"]], main = "Observada")

Podemos ver que de manera general, las tres variables sí influencian la riqueza, aunque los mecanismos precisos no los hemos dilucidado aún.

Regresar al índice del curso