Sistemas de Información Geográfica con R

Modelos con Raster

knitr::opts_chunk$set(fig.align = "center",
                      warning = FALSE,
                      message = FALSE,
                      fig.width = 9)

• Notas tomadas del curso “Sistemas de Información Geográfica en R (SIG en R)” con Derek Corcoran.

Bibliotecas

library(rworldxtra)
library(tidyverse)
library(sf)
library(raster)

Datos Solenopsis

• Ver anexos para obtener código de la base de datos que se carga a continuación.
• Estos datos están en este directorio de Github.

datos <- read_csv("https://raw.githubusercontent.com/Edimer/Spatial-Data-Science/main/SIG_R/data/solenopsis.csv")
datos

Transformación a sf

• En el argumento coords se incorpora un vector con la posición de las columnas longitud y latitud, respectivamente.
• El sistema de coordenadas elegido en el siguiente código es el más sencillo de todos, sin embargo, se puede cambiar (ver clase 01).

# Datos a sf
datos_sf <- datos %>% 
  st_as_sf(coords = c(5, 6), crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")

Área de estudio

Área en Colombia

colombia <- raster::getData(name = "GADM", country = "COL", level = 0)
colombia_sf <- st_as_sf(colombia)
colombia_sf %>% 
  ggplot() + 
  geom_sf() +
  geom_sf(data = datos_sf)

Área específica

• El siguiente mapa corresponde a una región específica de Colombia.

# Mapa del mundo
data("countriesHigh")
mapa <- countriesHigh %>% 
  st_as_sf() %>%
  st_crop(datos_sf)
ggplot() +
  geom_sf(data = mapa) +
  geom_sf(data = datos_sf)

Otros mapas

• Gráfico con color por especie y tamaño por abundancia: filtro los NA para que no aparezcan en la leyenda.

ggplot() +
  geom_sf(data = mapa) +
  geom_sf(data = datos_sf %>% filter(!is.na(Species)),
          aes(color = Species, size = Measurement))

• Gráfico por especie: filtro los NA para que no aparezcan en la leyenda.

ggplot() +
  geom_sf(data = mapa) +
  geom_sf(data = datos_sf %>% filter(!is.na(Species)),
          aes(color = Species, size = Measurement)) +
  facet_wrap(~Species)

Resumen descriptivo

datos %>% 
  group_by(Species) %>% 
  summarise(Total = n(),
            AbundanciaPromedio = mean(Measurement),
            MáximoAbundancia = max(Measurement),
            MínimoAbundancia = min(Measurement))

Modelación

Especie Geminata

• Para los ejemplos siguientes sólo uso la especie Geminata. Se observa mayor abundancia en el norte de Colombia.

data_geminata <- datos_sf %>% filter(Species == "geminata")
ggplot() +
  geom_sf(data = mapa) +
  geom_sf(data = data_geminata, aes(size = Measurement))

Datos ambientales

• Aquí puede encontrar información acerca de las variables bioclimáticas que se descargan con el siguiente código.
• Se realiza el corte del área de interés con la función crop() del paquete raster.
  • Cuando utilizamos formato shape usamos la función st_crop y con formatos raster utilizamos la función crop().

clima <- getData("worldclim", var = "bio", res = 2.5)

Error in getData("worldclim", var = "bio", res = 2.5) : 
  unused arguments (var = "bio", res = 2.5)

• Número de capas o variables bioclimáticas:

nlayers(clima2)

[1] 19

• Para el ejemplo sólo uso las capas 1, 7, 12, 15 y 19. La elección la hago de forma aleatoria para el ejemplo.

clima3 <- clima2[[c(1, 7, 12, 15, 19)]]

Extracción de datos

• Variables bioclimáticas: algunos puntos quedan con valores NA para las variables bioclimáticas.

clima_data <- extract(clima3, data_geminata) %>% 
  as.data.frame()
clima_data

• Uniendo datos de abundancia con variables bioclimáticas:

data_final <- data_geminata %>% 
  bind_cols(clima_data)
data_final

Simple feature collection with 52 features and 12 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -76.52944 ymin: 3.107725 xmax: -72.65786 ymax: 11.09833
CRS:            +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0
First 10 features:
                             Reference Study_common_taxon Site_name Sampling_effort  Species
1  Dominguez-Haydar and Armbrecht 2010         Formicidae      0yrs              20 geminata
2  Dominguez-Haydar and Armbrecht 2010         Formicidae      1yrs              20 geminata
3  Dominguez-Haydar and Armbrecht 2010         Formicidae      2yrs              20 geminata
4  Dominguez-Haydar and Armbrecht 2010         Formicidae      4yrs              20 geminata
5  Dominguez-Haydar and Armbrecht 2010         Formicidae      6yrs              20 geminata
6  Dominguez-Haydar and Armbrecht 2010         Formicidae      7yrs              20 geminata
7  Dominguez-Haydar and Armbrecht 2010         Formicidae      8yrs              20 geminata
8  Dominguez-Haydar and Armbrecht 2010         Formicidae     12yrs              20 geminata
9  Dominguez-Haydar and Armbrecht 2010         Formicidae     13yrs              20 geminata
10 Dominguez-Haydar and Armbrecht 2010         Formicidae     14yrs              20 geminata
   Measurement Effort_corrected_measurement bio1 bio7 bio12 bio15 bio19
1            3                            3   NA   NA    NA    NA    NA
2           90                           90  275  130  1370    71   204
3          128                          128  275  130  1370    71   204
4            5                            5   NA   NA    NA    NA    NA
5            3                            3   NA   NA    NA    NA    NA
6           13                           13   NA   NA    NA    NA    NA
7            0                            0   NA   NA    NA    NA    NA
8           12                           12   NA   NA    NA    NA    NA
9            0                            0   NA   NA    NA    NA    NA
10           0                            0   NA   NA    NA    NA    NA
                     geometry
1    POINT (-72.71417 11.085)
2  POINT (-72.71222 11.07944)
3  POINT (-72.71443 11.08131)
4   POINT (-72.7235 11.09015)
5    POINT (-72.71167 11.095)
6  POINT (-72.70503 11.09543)
7  POINT (-72.71333 11.09389)
8  POINT (-72.69667 11.09833)
9  POINT (-72.69019 11.09255)
10     POINT (-72.6975 11.09)

Relaciones

• Nota: para esta especie en particular el número de datos es bajo.

Lineales

data_final %>% 
  as.data.frame() %>% 
  dplyr::select(bio1:bio19, Measurement) %>% 
  gather(key = "key", value = "value", -Measurement)  %>% 
  ggplot(aes(x = value, y = Measurement)) +
  facet_wrap(~key, scales = "free", ncol = 3) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  theme_bw()

Lineales (Logaritmos de abundancia)

data_final %>% 
  as.data.frame() %>% 
  dplyr::select(bio1:bio19, Measurement) %>% 
  gather(key = "key", value = "value", -Measurement)  %>% 
  ggplot(aes(x = value, y = Measurement)) +
  facet_wrap(~key, scales = "free", ncol = 3) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  scale_y_log10() +
  theme_bw()

Modelos

• Se prueban tres modelos y basado en el criterio de información de Akaike se elige el mejor.
  • El primero modelo es una regresión poisson, es decir, un modelo lineal generalizado con distribución de errores Poisson.
  • El segundo es un modelo de regresión poisson con un polinomio de segundo grado para bio1 y bio15.
  • El tercero es un modelo aditivo generalizado (GAM).

GLM Poisson

modelo1 <- glm(Measurement ~ bio1 + bio7 + bio12 + bio15 + bio19,
               data = data_final, family = "poisson")
summary(modelo1)

Call:
glm(formula = Measurement ~ bio1 + bio7 + bio12 + bio15 + bio19, 
    family = "poisson", data = data_final)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-12.0620   -1.6910   -1.2246    0.5594    5.8436  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) 55.808756  23.859812   2.339 0.019334 *  
bio1        -0.080275   0.075201  -1.067 0.285763    
bio7        -0.326484   0.084470  -3.865 0.000111 ***
bio12        0.008604   0.002602   3.307 0.000942 ***
bio15        0.087892   0.065417   1.344 0.179087    
bio19       -0.024663   0.008683  -2.840 0.004506 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 1169.74  on 38  degrees of freedom
Residual deviance:  309.68  on 33  degrees of freedom
  (13 observations deleted due to missingness)
AIC: 378.12

Number of Fisher Scoring iterations: 6

• Residuales:

par(mfrow = c(2,2))
plot(modelo1)

GLM Poisson + Polinomio 2

modelo2 <- glm(Measurement ~ I(bio1^2) + bio7 + bio12 + I(bio15^2) + bio19,
               data = data_final, family = "poisson")
summary(modelo2)

Call:
glm(formula = Measurement ~ I(bio1^2) + bio7 + bio12 + I(bio15^2) + 
    bio19, family = "poisson", data = data_final)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-12.0616   -1.7216   -1.3000    0.5815    5.8441  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) 43.8725633 15.4543461   2.839 0.004528 ** 
I(bio1^2)   -0.0001434  0.0001564  -0.917 0.358979    
bio7        -0.3003849  0.0798584  -3.761 0.000169 ***
bio12        0.0076670  0.0027049   2.834 0.004590 ** 
I(bio15^2)   0.0008692  0.0007621   1.141 0.254019    
bio19       -0.0224261  0.0092793  -2.417 0.015658 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 1169.74  on 38  degrees of freedom
Residual deviance:  310.28  on 33  degrees of freedom
  (13 observations deleted due to missingness)
AIC: 378.72

Number of Fisher Scoring iterations: 6

• Residuales:

par(mfrow = c(2,2))
plot(modelo2)

GAM

library(mgcv)
modelo3 <- gam(Measurement ~ s(bio1, k = 4) + s(bio7, k = 5) + 
                 bio12 + bio15 + bio19 + bio1 + bio7,
               data = data_final, family = "poisson")
summary(modelo3)

Family: poisson 
Link function: log 

Formula:
Measurement ~ s(bio1, k = 4) + s(bio7, k = 5) + bio12 + bio15 + 
    bio19 + bio1 + bio7

Parametric coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   0.0000     0.0000      NA       NA
bio12         0.2318     1.2777   0.181    0.856
bio15         4.0208     5.0339   0.799    0.424
bio19        -0.6589     4.8372  -0.136    0.892
bio1         -4.5917    23.6431  -0.194    0.846
bio7          7.1399    46.0812   0.155    0.877

Approximate significance of smooth terms:
          edf Ref.df Chi.sq p-value
s(bio1) 0.889  1.019  0.534   0.505
s(bio7) 3.963  3.966  2.168   0.667

Rank: 11/13
R-sq.(adj) =  0.495   Deviance explained =   76%
UBRE = 6.6983  Scale est. = 1         n = 39

• Residuales:

par(mfrow = c(2, 2))
gam.check(modelo3)

Method: UBRE   Optimizer: outer newton
full convergence after 11 iterations.
Gradient range [2.235383e-06,9.044007e-06]
(score 6.698331 & scale 1).
Hessian positive definite, eigenvalue range [9.541386e-06,0.0005525887].
Model rank =  11 / 13 

Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.

           k'   edf k-index p-value
s(bio1) 3.000 0.889    0.99    0.49
s(bio7) 4.000 3.963    1.22    0.94

Mejor Modelo

• Es mejor el modelo GAM (AIC más bajo). Con este modelo se realizan las predicciones de abundancia.

AIC(modelo1, modelo2, modelo3)

Predicciones

• Nota: dado el bajo número de datos (51 registros) para la especie bajo análisis, no se logra construir un modelo que proporcione abundancias predichas correctamente. Cabe mencionar que es posible mejorar el análisis incluyendo más información de la especie o ajustando otros modelos.

prediccion <- predict(clima3, modelo2, type = "response")
plot(prediccion, colNA = "black")

• Mapa con ggplot2:

prediccion_df <- prediccion %>% 
  as("SpatialPixelsDataFrame") %>% 
  as.data.frame() %>% 
  rename(Abundancia = layer)

ggplot() +
  geom_tile(data = prediccion_df, aes(x = x, y = y, fill = Abundancia)) +
  geom_sf(data = mapa, alpha = 0) +
  scale_fill_viridis_c() +
  theme_bw()

Predicciones futuras

• Imágenes futuras: las mismas con las cuales se entrenó el modelo.

futuro <- raster::getData("CMIP5", var = "bio", res = 2.5, rcp = 85,
                          model = "HD", year = 70) %>% 
  crop(data_geminata)

futuro <- futuro[[c(1, 7, 12, 15, 19)]]
names(futuro) <- c("bio1", "bio7", "bio12", "bio15", "bio19")
plot(futuro)

• Mapa de predicción de abundancia:

prediccion_futura <- predict(futuro, modelo2, type = "response")

prediccion_df_futuro <- prediccion_futura %>% 
  as("SpatialPixelsDataFrame") %>% 
  as.data.frame() %>% 
  rename(Abundancia = layer)

ggplot() +
  geom_tile(data = prediccion_df_futuro, aes(x = x, y = y, fill = Abundancia)) +
  geom_sf(data = mapa, alpha = 0) +
  scale_fill_viridis_c() +
  theme_bw()

Anexos

• En el ejemplo de clase se trabaja con la especie Bombus affinis (abejorros), sin embargo, obtuve directamente la base de datos de PREDICTS database y voy a trabajar con el género de hormigas Solenopsis, también conocida como hormiga colorada o de fuego..
• En este caso filtro información sólo para Colombia.
• Se filtran datos sólo de abundancia de especies.
• Mantengo sólo las siguientes variables:
  • Reference: referencia del estudio
  • Study_common_taxon: taxón
  • Site_name: nombre de sitio
  • Sampling_effort
  • Longitude
  • Latitude
  • Species
  • Measurement: medida de abundancia
  • Effort_corrected_measurement

# Base de datos total con más de 3 millones de registros
datos_total <- fread("database.csv", encoding = "UTF-8")
# Filtro solenopsis para Colombia: 551 registros de abundancia
datos_colombia <- datos_total %>% 
  filter(Country == "Colombia") %>% 
  filter(Genus == "Solenopsis") %>%  
  filter(Diversity_metric_type == "Abundance") %>% 
  select(Reference, Study_common_taxon, Site_name, Sampling_effort,
         Longitude, Latitude, Species, Measurement, 
         Effort_corrected_measurement) 
# Exportando datos de solenopsis
write.csv(datos_colombia, file = "solenopsis.csv", row.names = FALSE,
          fileEncoding = "UTF-8")