Sistemas de Información Geográfica con R

Autocorrelación Espacial

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

Bibliotecas

library(tidyverse)
library(sp)
library(sf)
library(gstat)
library(rgdal)
library(rworldxtra)
library(raster)

Ley de Tobler

“Todas las cosas están relacionadas entre sí, pero las cosas más próximas en el espacio tienen una relación mayor que las distantes.” Waldo Tobler

Bibliotecas

library(gstat)
library(raster)
library(rgdal)
library(rworldxtra)
library(sp)
library(sf)
library(tidyverse)

Datos meuse

• Datos de ejemplo contenidos en la biblioteca sp.
• El crs asignado es determinado por las unidades en las cuales está dada la longitud y latitud.
• En este caso se transforman los datos a sf y a SpatialPoints.

# Data originakl
data("meuse") 

# Data sf
meuse_sf <- meuse %>% st_as_sf(coords = c(1, 2), crs = "+init=epsg:28992")

# Data SpatialPoints
coordinates(meuse) <- ~ x + y

# Clase de objetos
class(meuse)

[1] "SpatialPointsDataFrame"
attr(,"package")
[1] "sp"

class(meuse_sf)

[1] "sf"         "data.frame"

• Concentración de Zinc:

meuse_sf %>% 
  ggplot() +
  geom_sf(aes(color = zinc)) +
  theme_bw() +
  scale_color_viridis_c()

Ubicación meuse

library(leaflet)
leaflet(as_Spatial(st_transform(meuse_sf, crs = 4326))) %>% addTiles() %>% addCircles()

Variograma

• El Variograma es una medida de variación entre valores a distintas distancias.
• A mayor distancia se espera mayor variación.
• El logaritmo se adiciona para garantizar la obtención de valores positivos.
• Es posible modelar la variación en función de diferentes fuentes o factores, dando lugar con ello a diferentes variogramas.

Modelo Nulo

#Modelo nulo
variograma_nulo <- variogram(object = log(zinc) ~ 1, data = meuse)
variograma_nulo

• Resultados de tabla anterior: la primera variable (np) proporciona información del número de puntos, es decir, que para la primera fila diremos que hay 57 puntos a una distancia (dist) de 79.29244 metros y la semi-varianza entre estos 57 puntos es 0.1234479 (gamma). El punto de mayor variación se da con 0.7033984 a una distancia de 1117.86 metros.

variograma_nulo %>% 
  ggplot(aes(x = dist, y = gamma)) +
  geom_point() +
  theme_bw() +
  labs(x = "Distancia (metros)", y = "Gamma (semi-varianza)")

Comparando variogramas

# Variogramas
var_nulo <- variogram(log(zinc) ~ 1, data = meuse) %>% mutate(Modelo = "Nulo")
var_spat <- variogram(log(zinc) ~ x + y, data = meuse) %>% mutate(Modelo = "Espacial")
var_dist <- variogram(log(zinc) ~ dist, data = meuse) %>% mutate(Modelo = "Distancia")
var_sqrt <- variogram(log(zinc) ~ sqrt(dist), data = meuse) %>% mutate(Modelo = "Raíz Distancia")

# Uniendo resultados
all_var <- list(var_nulo, var_spat, var_dist, var_sqrt) %>% 
  reduce(bind_rows)
all_var

• Gráfico de variogramas:

all_var %>%
  ggplot(aes(x = dist, y = gamma, color = Modelo)) +
  geom_point(size = 1.7) +
  theme_bw() +
  labs(x = "Distancia (metros)", y = "Gamma (semi-varianza)")

Ajuste de variograma

• psill: meseta parcial. Estabilización de la semivarianza.
• range: separación o distancia entre pares de puntos en la cual ya no hay dependencia espacial.
• nugget: semivarianza esperada a la separación o distancia de 0m
• Los valores de psill, range y nugget sirven como valores iniciales para el ajuste del modelo.

ajuste_var <- fit.variogram(object = var_spat,
                            model = vgm(psill = 1, model = "Sph", range =  700,
                                        nugget =  1))
ajuste_var

• Resultados: el range proporcionado por el modelo indica la distancia esperada donde se estabiliza la semivarianza. La suma de psill (0.08234213 + 0.38866509) dará como resultado el valor donde se evidencia o espera la “meseta” total (estabilización de la semivarianza). Nugget brinda información acerca de la semivarianza esperada a una distancia de 0 metros (intercepto). A continuación se muestra el gráfico del modelo ajustado:

ggplot(variogramLine(ajuste_var, 1600), aes(x = dist, y = gamma)) +
  geom_path() + 
  geom_point(data = var_spat, color = "red") + 
  geom_vline(xintercept = ajuste_var$range[2], lty = 2) + 
  geom_text(x = ajuste_var$range[2], y = ajuste_var$psill[2]/2, label = "range") + 
  theme_bw() +
  geom_hline(yintercept = ajuste_var$psill[2] + ajuste_var$psill[1] , lty = 2) + 
  geom_text(x = ajuste_var$range[2]/2, y = ajuste_var$psill[2] + ajuste_var$psill[1], label = "psill") +
  geom_text(x = ajuste_var$range[1], y = ajuste_var$psill[1], label = "Nugget") +
  ylim(c(0,max(var_spat$gamma)))

Otros modelos

# Modelo nulo
ajuste_nulo <- fit.variogram(object = var_nulo,
                             model = vgm(psill = 1, model = "Sph", range =  700,
                                         nugget =  1))

# Modelo con distancia
ajuste_dist <- fit.variogram(object = var_dist,
                            model = vgm(psill = 1, model = "Sph", range =  700,
                                        nugget =  1))

# Modelo con raíz cuadrada de la distancia
ajuste_sqrt <- fit.variogram(object = var_sqrt,
                            model = vgm(psill = 1, model = "Sph", range =  700,
                                        nugget =  1))

• Gráfico con cuatro modelos ajustados:

# Cuatro modelos
Abn_fit_null <- variogramLine(ajuste_nulo, 1600) %>% mutate(Modelo = "Nulo")
Abn_line_Spat <- variogramLine(ajuste_var, 1600) %>% mutate(Modelo = "Espacial")
Abn_line_Dist <- variogramLine(ajuste_dist, 1600) %>% mutate(Modelo = "Distancia")
Abn_line_Dist_sq <- variogramLine(ajuste_sqrt, 1600) %>% mutate(Modelo = "Raíz Distancia")

# Unión de datos
Abn_line <- list(Abn_fit_null, Abn_line_Spat, Abn_line_Dist, Abn_line_Dist_sq) %>%
  reduce(bind_rows)

ggplot(Abn_line, aes(x = dist, y = gamma)) + 
  geom_path(aes(color = Modelo)) +
  geom_point(data = all_var, aes(color = Modelo)) + 
  theme_bw() +
  scale_color_brewer(palette = "Set1")

Kriging

• El krigeaje o regresión en procesos gaussianos permite realizar interporlación a partir de un variograma previamente ajustado.
• En este caso se usa la base de datos meuse.grid que tiene información que podrá ser utilizada para la interpolación. Las distancias del siguiente gráfico indican proximidad a una fuente hídrica.

Datos

# Data inicial
data("meuse.grid")

# Conversión a sf
meusegrid_sf <- meuse.grid %>% 
  st_as_sf(coords = c(1, 2), crs = "+init=epsg:28992")

# Gráfico
meusegrid_sf %>% 
  ggplot() +
  geom_sf(aes(color = dist)) +
  scale_color_viridis_c() +
  theme_bw()

Predicciones (interpolación)

pred_spat <- krige(log(zinc) ~ 1, meuse_sf, meusegrid_sf, model = ajuste_var)

[using ordinary kriging]

class(pred_spat)

[1] "sf"         "data.frame"

pred_spat

Simple feature collection with 3103 features and 2 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 178460 ymin: 329620 xmax: 181540 ymax: 333740
projected CRS:  Amersfoort / RD New
First 10 features:
   var1.pred  var1.var              geometry
1   6.535255 0.2461491 POINT (181180 333740)
2   6.626923 0.2093807 POINT (181140 333700)
3   6.528438 0.2197185 POINT (181180 333700)
4   6.429889 0.2315686 POINT (181220 333700)
5   6.730607 0.1709532 POINT (181100 333660)
6   6.622317 0.1818457 POINT (181140 333660)
7   6.508323 0.1944079 POINT (181180 333660)
8   6.396911 0.2074667 POINT (181220 333660)
9   6.835826 0.1364733 POINT (181060 333620)
10  6.728034 0.1437661 POINT (181100 333620)

• Gráfico de predicciones: como los datos están en logaritmos se usa la función exp() para obtener valores en las mismas unidades

pred_spat %>% 
  ggplot() +
  geom_sf(aes(color = exp(var1.pred))) +
  scale_color_viridis_c(name = "Zinc (predicho)")

• Gráfico de variación en las predicciones: el gráfico permite evidencicar lugares de mayor incertidumbre para la predicción.

pred_spat %>% 
  ggplot() +
  geom_sf(aes(color = exp(var1.var))) +
  scale_color_viridis_c(name = "")

Todos los modelos

Null_pred <- krige(log(zinc) ~ 1,  meuse_sf, meusegrid_sf, model = ajuste_nulo) %>%
  mutate(Modelo = "Nulo")

[using ordinary kriging]

Spat_pred <- krige(log(zinc) ~ 1, meuse_sf, meusegrid_sf, model = ajuste_var) %>%
  mutate(Modelo = "Espacial")

[using ordinary kriging]

Dist_pred <- krige(log(zinc) ~ 1,  meuse_sf, meusegrid_sf, model = ajuste_dist) %>% 
  mutate(Modelo = "Distancia")

[using ordinary kriging]

Dist_sq_pred <- krige(log(zinc) ~ 1, meuse_sf, meusegrid_sf, model = ajuste_sqrt)  %>%
  mutate(Modelo = "Raíz Distancia")

[using ordinary kriging]

Pred <- list(Null_pred, Spat_pred, Dist_pred, Dist_sq_pred) %>%
  reduce(bind_rows)
Pred

Simple feature collection with 12412 features and 3 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 178460 ymin: 329620 xmax: 181540 ymax: 333740
projected CRS:  Amersfoort / RD New
First 10 features:
   var1.pred  var1.var Modelo              geometry
1   6.499611 0.3198081   Nulo POINT (181180 333740)
2   6.622347 0.2520183   Nulo POINT (181140 333700)
3   6.505156 0.2729842   Nulo POINT (181180 333700)
4   6.387580 0.2955285   Nulo POINT (181220 333700)
5   6.764490 0.1779374   Nulo POINT (181100 333660)
6   6.635508 0.2022009   Nulo POINT (181140 333660)
7   6.497544 0.2277371   Nulo POINT (181180 333660)
8   6.361471 0.2524902   Nulo POINT (181220 333660)
9   6.904628 0.1099817   Nulo POINT (181060 333620)
10  6.780233 0.1280782   Nulo POINT (181100 333620)

• Predicciones:

ggplot() +
  geom_sf(data = Pred, aes(color = exp(var1.pred))) +
  scale_color_viridis_c(name = "Zinc") +
  facet_wrap(~Modelo)  + 
  theme_bw()

• Incertidumbre en predicciones: el modelo de mayor incertidumbre es el modelo nulo.

ggplot() +
  geom_sf(data = Pred, aes(color = exp(var1.var))) +
  scale_color_viridis_c(name = "Zinc") +
  facet_wrap(~Modelo)  + 
  theme_bw()

Mejor modelo

• Una manera de seleccionar el modelo podría ser comparar la bondad de ajuste a través de una métrica que represente el error, por ejemplo, Raíz del Cuadrado Medio del Error (RSME - Root Square Mean Error).
• También es posible implementar validación cruzada.
• La función krige.cv permite ajustar modelos a través de validación cruzada.

Null_pred_cv <- krige.cv(log(zinc) ~ 1,  meuse_sf,
                         model = ajuste_nulo, nfold = 5) %>%
  mutate(Modelo = "Nulo")

Spat_pred_cv <- krige.cv(log(zinc) ~ 1, meuse_sf,model = ajuste_var,
                         nfold = 5) %>%
  mutate(Modelo = "Espacial")

Dist_pred_cv <- krige.cv(log(zinc) ~ 1,  meuse_sf, model = ajuste_dist,
                         nfold = 5) %>% 
  mutate(Modelo = "Distancia")

Dist_sq_pred_cv <- krige.cv(log(zinc) ~ 1, meuse_sf, model = ajuste_sqrt,
                         nfold = 5)  %>%
  mutate(Modelo = "Raíz Distancia")

Pred_cv <- list(Null_pred_cv, Spat_pred_cv, Dist_pred_cv, Dist_sq_pred_cv) %>%
  reduce(bind_rows)

• Predichos vs observados:

Pred_cv %>% 
  ggplot(aes(x = var1.pred, y = observed)) +
  geom_point() +
  theme_bw() +
  labs(x = "Predichos", y = "Observados")

• Comparando modelos a través de RMSE:

Pred_cv %>% 
  as.data.frame() %>% 
  group_by(Modelo) %>% 
  summarise(RMSE = sqrt(sum(residual^2)/length(residual)))