Geoestadistica

Author

Marco Mora Pablo Pedraza Pedro Straub

Taller 3: Análisis de Autocorrelación espacial

1 SETUP

1.1 Cargar Librerías

library(tidyverse)
library (spdep)
library(terra)
library(sf)
library(gstat)
library(dplyr)
library(ggplot2)
library(tmap)

2 Febrero

2.1 Indice de Moran Global

2.1.1 Carga de datos

# Cargar datos
data_t_feb <- read_rds('C:/Re6/Geoestadistica/Talleres/data/procesada/data_temp_febrero.rds')

# Cálculo de la matriz de distancia

w <- 1/as.matrix(st_distance(data_t_feb))
diag(w) <- 0
w <- units::drop_units(w)
w <- mat2listw(w)
moran(data_t_feb$temp,w,n=length(w$neighbours),S0=Szero(w))

#Volvemos a calcular moran
w <- 1/as.matrix(st_distance(data_t_feb))
diag(w) <- 0
w <- units::drop_units(w)
w <- mat2listw(w)
m <- moran(data_t_feb$temp,w,n=length(w$neighbours),S0=Szero(w))
mp <- moran.plot(data_t_feb$temp,w)

# Buscar las estaciones que se encuentran muy cercanas
d <- st_distance(data_t_feb) #calcula la distance entre estaciones
diag(d) <- NA #convierte la diagonal de la matriz de distancia en NA
d <- units::drop_units(d) # quita las unidades de la distancia
ind <- which(d < 100,arr.ind = TRUE) # busca que estaciones se encuentran a menos de 100m
data_t_febf <- data_t_feb[-ind[,1],] #elimina las estaciones que estan a una distancia de menos de 100m
w <- 1/as.matrix(st_distance(data_t_feb))
diag(w) <- 0
w <- units::drop_units(w)
w <- mat2listw(w)
moran(data_t_feb$temp,w,n=length(w$neighbours),S0=Szero(w))

#Test estadístico Montecarlo para probar la significancia de la autocorrelación espacial
moran.mc(data_t_feb$temp,w,nsim=10000)
mp <- moran.plot(data_t_feb$temp,w) #Con esto se descarta la hipotesis nula.

Al realizar la prueba de Montecarlo para el índice global de Moran se puede afirmar que existe una correlación espacial positiva.
Porque la distribución de los datos tienen un valor estadístico de 0.23206 con un p-value de 0.0001 tras 10000 simulaciones.

2.2 Índice de Moran local

lmoran <- localmoran(data_t_feb$temp,w,alternative = "greater")
head(lmoran)

2.3 Mapa de autocorrelacion

data_t_feb <- st_as_sf(data_t_feb[c('geometry','temp')],coords = c('x','y'),crs = 32719)
data_t_feb$Ii <- lmoran[,'Ii']
data_t_feb$Z <- lmoran[,'Z.Ii']
data_t_feb$Pr <- lmoran[, "Pr(z > E(Ii))"]

p1 <- tm_shape(data_t_feb) +
  tm_dots(col = "temp", title = "Temperatura", style = "quantile", palette = "viridis") +
  tm_layout(legend.position = c("0.9", "0.1"), legend.outside = FALSE, frame = FALSE) +
   tm_scale_bar(position = c("right", "bottom"), text.size = 0.5) +
  tm_grid(proj = 4326, n.x = 8, n.y = 8, lines = TRUE, col = "black", lwd = 0.5, alpha = 0.5)  

# Repite para los otros mapas
p2 <- tm_shape(data_t_feb) +
  tm_dots(col = "Ii", title = "Índice de Moran Local", style = "quantile", palette = "viridis") +
  tm_layout(legend.position = c("0.9", "0.1"), legend.outside = FALSE, frame = FALSE) +
   tm_scale_bar(position = c("right", "bottom"), text.size = 0.5) +
  tm_grid(proj = 4326, n.x = 8, n.y = 8, lines = TRUE, col = "black", lwd = 0.5, alpha = 0.5)  

p3 <- tm_shape(data_t_feb) +
  tm_dots(col = "Z", title = "Z-score", style = "fixed", breaks = c(-Inf, 1.65, Inf), palette = "viridis") +
  tm_layout(legend.position = c("0.9", "0.1"), legend.outside = FALSE, frame = FALSE) +
    tm_scale_bar(position = c("right", "bottom"), text.size = 0.5) +
  tm_grid(proj = 4326, n.x = 8, n.y = 8, lines = TRUE, col = "black", lwd = 0.5, alpha = 0.5)  

p4 <- tm_shape(data_t_feb) +
  tm_dots(col = "Pr", title = "p-value", style = "fixed", breaks = c(-Inf, 0.05, Inf), palette = "viridis") +
  tm_layout(legend.position = c("0.9", "0.1"), legend.outside = FALSE, frame = FALSE) +
    tm_scale_bar(position = c("right", "bottom"), text.size = 0.5) +
  tm_grid(proj = 4326, n.x = 8, n.y = 8, lines = TRUE, col = "black", lwd = 0.5, alpha = 0.5)

tmap_mode('view')
tmap_arrange(p1,p2,p3,p4)

Mapa Interactivo

Mapa Estático

#Quad

quad <- attr(lmoran, 'quadr') 
head(quad)

data_t_feb <- data_t_feb |> 
  bind_cols(quad = quad$median) |>  # Agregando columna quad
  mutate(quad = case_when(
    Pr < 0.05 ~ quad,              # Mantiene valores significativos
    .default = 'No Significativo')) # Asigna "No Significativo" cuando Pr >= 0.05

# Crear el mapa con la columna "quad"
moran <- tm_shape(data_t_feb) +
  tm_dots(col = "quad", title = "Quad", style = "cat", palette = "Set1") +  # Cambiar a 'quad' y estilo 'cat' para categorías
  tm_layout(legend.outside = TRUE)

tmap_mode('view')
moran

El mapa del índice moran Quad nos muestra la ubicación y concentración de puntos de temperatura en autocorrelación espacial a sus vecinos.
Los puntos en rojo (High-High) son aquellos con valores altos de temperatura entorno a la media, rodeados por puntos vecinos con temperaturas similares.Son puntos con autocorrelación espacial positva.
Los puntos en azul (High-Low) son aquellos con valores altos de temperatura entorno a la media, rodeados por puntos vecinos con temperaturas más bajas. Son puntos con autocorrelación espacial negativa.
Los puntos en rojo (Low-High) son aquellos con valores bajos de temperatura entorno a la media, rodeados por puntos vecinos con temperaturas más altas.Son puntos con autocorrelación espacial negativa.
Los puntos en azul (Low-Low) son aquellos con valores bajos de temperatura entorno a la media, rodeados por puntos vecinos con temperaturas similares. Son puntos con autocorrelación espacial positiva.
Los puntos en naranjo (No Significativo) son aquellos cuyos valores de temperatura tienen un comportamiento aleatorio en relación a sus vecionos. Son puntos sin autocorrelación espacial.
Al comparar el índice de moran global y el índice de moran local, se pueden afirmar dos cosas principalmente:
La autocorrelación espacial general es positiva. Eso dice el índice de moral global (0.23206), lo que se puede corroborar visualmente con la presencia mayoritaria de puntos rojos (High-High) y morados (Low-Low).
Que los mayores clusters de autocorrelación se encuentran en las cordilleras de la Costa (High-High) y de los Andes (Low-Low).

2.4 Análisis del Variograma

2.4.1 Variograma Experimental

options(scipen = 10)

# Realizamos un variograma exploratorio para identificar qué variograma se ajusta más a nuestros datos
v <- variogram(temp~1,data_t_feb, cloud = TRUE) 
plot(v)

# calculamos la distancia de la diagonal para tener una referencia de cuanto debiese ser nuestro Cutoff
distance_matrix <- st_distance(data_t_feb)

max_dist <- max(dist_matrix, na.rm = TRUE)

diagonal <- max_dist/3
print(diagonal)

#Definimos el punto de corte.
v_feb <- variogram(temp~1,data_t_feb,cutoff = 90000) 
plot(v_feb)
show.vgms() #Verificar tipos de variogramas teoricos para ver cuál se ajusta más

# Definimos el sill (máximo valor de los datos), el tipo de variograma (Wav), y el nugget (la variabilidad inicial de los datos).
m_var_feb <- vgm(5.3,'Wav',82000,0.1)

# Seguimos reiterando hasta encontrar un variograma que se ajuste a nuestros datos e indique la distancia donde la semivarianza se estabiliza.

plot(v_feb,m_var_feb)

# El variograma experimental de Febrero nos indica que la máxima semiviaranza es de alredor de 5.3 y que se estabiliza a los 80000 m, y se asemaja un variograma de onda (Wav).

# Guardamos
png( "C:/Re6/Geoestadistica/Talleres/Outputs/Grafico/var_febrero.png",
       width = 25, height = 14, units = "cm", res = 300)
plot(v_feb,m_var_feb, 
     col = "blue", lwd = 1,
     main = "Variograma de Temperatura para Febrero 2023",
     xlab = "Distancia (m)",
     ylab = "Semivarianza")
dev.off()

write_rds(m_var_feb,'C:/Re6/Geoestadistica/Talleres/data/procesada/var_feb.rds')

2.5 Análisis de Anisotropía

# Definir variogramas en diferentes direcciones
v_febrero_dir <- variogram(temp ~ 1, data_t_feb, alpha = c(0, 45, 90, 135), cutoff = 120000, width=6000)
plot(v_febrero_dir)

El gráfico 0° presenta una curva plana y una casi nula variabilidad.
El gráfico 45° presenta un patrón un poco más definido, con una máxima semivarianza de 5 y que deja aumentar al rededor los 58000 m.
EL gráfico 90° presenta un patrón muy defenido y ajustado a un variograma Gausiano, con una semivarianza máxima cercana a 11 y que se estabiliza rodeando los 90000 m.
EL gráfico 135° presenta una máxima variabilidad cercana a 5 y que se estabiliza al rededor de los 78000 m.

# Se optó por la anisotropía de 90° porque es la dirección que presenta el mejor ajuste.
v_anis_feb <- vgm(13.5,'Gau',80000,0.2,anis = c(90,1))
plot(v_febrero_dir,v_anis_feb)

Anisotropía espacial del variograma es Gausiano

# Mapa Variográfico
plot(variogram(temp~1,data_t_feb,cutoff = 120000, width= 8000, map = TRUE))

Del mapa se interpreta una mayor variabilidad de los datos en las cordilleras de la Costa y de los Andes.

2.5.1 Ajuste del variograma

# Guardar en .rds
fit_v_feb <- fit.variogram(v_feb,m_var_feb) # Ajuste del variograma teórico (mv) con el variograma experimental (v).

write_rds(fit_v_feb,'C:/Re6/Geoestadistica/Talleres/data/procesada/var_feb_ajustado.rds')

plot(v_feb, fit_v_feb, scales = list(x = list(log = FALSE), y = list(log = FALSE)),
     col = "blue", lwd = 1,
     main = "Variograma ajustado de Temperatura para Febrero 2023",
     xlab = "Distancia (m)",
     ylab = "Semivarianza")

png( "C:/Re6/Geoestadistica/Talleres/Outputs/Grafico/var_febrero_ajustado.png", # Guardamos el gráfico
     width = 25, height = 14, units = "cm", res = 300)  
plot(v_feb, fit_v_feb, scales = list(x = list(log = FALSE), y = list(log = FALSE)),
     col = "blue", lwd = 1,
     main = "Variograma ajustado de Temperatura para Febrero 2023",
     xlab = "Distancia (m)",
     ylab = "Semivarianza")
dev.off()

3 Julio

3.1 Indice de Moran Global

3.1.1 Carga de datos

# Cargar datos
data_t_jul <- read_rds('C:/Re6/Geoestadistica/Talleres/data/procesada/data_temp_julio.rds')

# Cálculo de la matriz de distancia

w <- 1/as.matrix(st_distance(data_t_jul))
diag(w) <- 0
w <- units::drop_units(w)
w <- mat2listw(w)
moran(data_t_jul$temp,w,n=length(w$neighbours),S0=Szero(w))

# filtrar datos por distancia entre las estaciones
data_t_jul <- data_t_jul|> 
    distinct(geometry,.keep_all = TRUE) #deja sólo las geometrias únicas

#Volvemos a calcular moran
w <- 1/as.matrix(st_distance(data_t_jul))
diag(w) <- 0
w <- units::drop_units(w)
w <- mat2listw(w)
m <- moran(data_t_jul$temp,w,n=length(w$neighbours),S0=Szero(w))
mp <- moran.plot(data_t_jul$temp,w)

# Buscar las estaciones que se encuentran muy cercanas
d <- st_distance(data_t_jul) #calcula la distance entre estaciones
diag(d) <- NA #convierte la diagonal de la matriz de distancia en NA
d <- units::drop_units(d) # quita las unidades de la distancia
ind <- which(d < 100,arr.ind = TRUE) # busca que estaciones se encuentran a menos de 100m
data_t_julf <- data_t_jul[-ind[,1],] #elimina las estaciones que estan a una distancia de menos de 100m
w <- 1/as.matrix(st_distance(data_t_jul))
diag(w) <- 0
w <- units::drop_units(w)
w <- mat2listw(w)
moran(data_t_jul$temp,w,n=length(w$neighbours),S0=Szero(w))

#Test estadístico Montecarlo para probar la significancia de la autocorrelación espacial

moran.mc(data_t_jul$temp,w,nsim=10000)
mp <- moran.plot(data_t_jul$temp,w) #Con esto se descarta la hipotesis nula.

Al realizar la prueba de Montecarlo para el índice global de Moran se puede afirmar que existe una correlación espacial positiva.
Porque la distribución de los datos tiene un valor estadístico de 0.26426 con un p-value de 0.0001 tras 10000 simulaciones.

3.2 Indice de Moran Local

lmoran <- localmoran(data_t_jul$temp,w,alternative = "greater")
head(lmoran)

3.3 Mapa de autocorrelacion

data_t_jul <- st_as_sf(data_t_jul[c('geometry','temp')],coords = c('x','y'),crs = 32719)
data_t_jul$Ii <- lmoran[,'Ii']
data_t_jul$Z <- lmoran[,'Z.Ii']
data_t_jul$Pr <- lmoran[, "Pr(z > E(Ii))"]

# Crear los mapas
p1 <- tm_shape(data_t_jul) +
  tm_dots(col = "temp", title = "Temperatura", style = "quantile", palette = "viridis") +
  tm_layout(legend.position = c("0.9", "0.1"), legend.outside = FALSE, frame = FALSE) +
  tm_scale_bar(position = c("right", "bottom"), text.size = 0.5) +
  tm_grid(proj = 4326, n.x = 8, n.y = 8, lines = TRUE, col = "black", lwd = 0.5, alpha = 0.5)  

# Repite para los otros mapas
p2 <- tm_shape(data_t_jul) +
  tm_dots(col = "Ii", title = "Índice de Moran Local", style = "quantile", palette = "viridis") +
  tm_layout(legend.position = c("0.9", "0.1"), legend.outside = FALSE, frame = FALSE) +
  tm_scale_bar(position = c("right", "bottom"), text.size = 0.5) +
  tm_grid(proj = 4326, n.x = 8, n.y = 8, lines = TRUE, col = "black", lwd = 0.5, alpha = 0.5)  

p3 <- tm_shape(data_t_jul) +
  tm_dots(col = "Z", title = "Z-score", style = "fixed", breaks = c(-Inf, 1.65, Inf), palette = "viridis") +
  tm_layout(legend.position = c("0.9", "0.1"), legend.outside = FALSE, frame = FALSE) +
  tm_scale_bar(position = c("right", "bottom"), text.size = 0.5) +
  tm_grid(proj = 4326, n.x = 8, n.y = 8, lines = TRUE, col = "black", lwd = 0.5, alpha = 0.5)  

p4 <- tm_shape(data_t_jul) +
  tm_dots(col = "Pr", title = "p-value", style = "fixed", breaks = c(-Inf, 0.05, Inf), palette = "viridis") +
  tm_layout(legend.position = c("0.9", "0.1"), legend.outside = FALSE, frame = FALSE) +
  tm_scale_bar(position = c("right", "bottom"), text.size = 0.5) +
  tm_grid(proj = 4326, n.x = 8, n.y = 8, lines = TRUE, col = "black", lwd = 0.5, alpha = 0.5)

tmap_mode('view')
tmap_arrange(p1,p2,p3,p4)

Mapa Interactivo

Mapa Estático

#Quad

quad <- attr(lmoran, 'quadr') 
head(quad)

data_t_jul <- data_t_jul |> 
  bind_cols(quad = quad$median) |>  # Agregando columna quad
  mutate(quad = case_when(
    Pr < 0.05 ~ quad,              # Mantiene valores significativos
    .default = 'No Significativo')) # Asigna "No Significativo" cuando Pr >= 0.05

# Crear el mapa con la columna "quad"
moran <- tm_shape(data_t_jul) +
  tm_dots(col = "quad", title = "Quad", style = "cat", palette = "Set1") +  # Cambiar a 'quad' y estilo 'cat' para categorías
  tm_layout(legend.outside = TRUE)

tmap_mode('view')
moran

El mapa del índice moran Quad nos muestra la ubicación y concentración de puntos de temperatura en autocorrelación espacial a sus vecinos.
Los puntos en rojo (High-High) son aquellos con valores altos de temperatura entorno a la media, rodeados por puntos vecinos con temperaturas similares.Son puntos con autocorrelación espacial positva
Los puntos en azul (High-Low) son aquellos con valores altos de temperatura entorno a la media, rodeados por puntos vecinos con temperaturas más bajas. Son puntos con autocorrelación espacial negativa.
Los puntos en rojo (Low-High) son aquellos con valores bajos de temperatura entorno a la media, rodeados por puntos vecinos con temperaturas más altas.Son puntos con autocorrelación espacial negativa.
Los puntos en azul (Low-Low) son aquellos con valores bajos de temperatura entorno a la media, rodeados por puntos vecinos con temperaturas similares. Son puntos con autocorrelación espacial positiva.
Los puntos en naranjo (No Significativo) son aquellos cuyos valores de temperatura tienen un comportamiento aleatorio en relación a sus vecionos. Son puntos sin autocorrelación espacial.
Al comparar el índice de moran global y el índice de moran local, se pueden afirmar dos cosas principalmente:
La autocorrelación espacial general es positiva. Eso dice el índice de moral global (0.26426), lo que se puede corroborar visualmente con la presencia mayoritaria de puntos rojos (High-High) y morados (Low-Low).
Que los mayores clusters de autocorrelación se encuentran en las cordilleras de la Costa (High-High) y de los Andes (Low-Low).

3.4 Análisis del variograma

3.4.1 Variograma Experimental

options(scipen = 10)

# Realizamos un variograma exploratorio para identificar qué variograma se ajusta más a nuestros datos
v <- variogram(temp~1,data_t_jul) 
plot(v)
show.vgms() #Verificar tipos de variogramas teoricos para ver cuál se ajusta más
v_jul <- variogram(temp~1,data_t_jul,cutoff = 70000) #Definimos el punto de corte.
plot(v_jul)

# Definimos el sill (máximo valor de los datos), el tipo de variograma (Wav), y el nugget (la variabilidad inicial de los datos).
m_var_jul <- vgm(6,'Wav',60000, 0) #

plot(v_jul,m_var_jul)

El variograma experimental de Julio nos indica que la máxima semiviaranza es de alredor de 7.8 y que se estabiliza a los 68000 m, y se asemaja un variograma de onda (Wav). Con el punto de corte 70000 m se elimaron los puntos que evidencian la estabilización de la semivarianza.

# Guardamos

png( "C:/Re6/Geoestadistica/Talleres/Outputs/Grafico/var_julio.png",
     width = 25, height = 14, units = "cm", res = 300)
plot(v_jul,m_var_jul, 
     col = "blue", lwd = 1,
     main = "Variograma de Temperatura para julio 2023",
     xlab = "Distancia (m)",
     ylab = "Semivarianza")
dev.off()

write_rds(m_var_jul,'C:/Re6/Geoestadistica/Talleres/data/procesada/var_jul.rds')

3.5 Análisis de Anisotropía

# Definir variogramas en diferentes direcciones
v_julio_dir <- variogram(temp ~ 1, data_t_jul, alpha = c(0, 45, 90, 135), cutoff = 60000, width=4000)
plot(v_julio_dir)

El gráfico 0° presenta una curva plana, y una casi nula variabilidad.
El gráfico 45° presenta una máxima semivarianza máxima cercana a 1.8 que se estabiliza al rededor de los 45000 m.
EL gráfico 90° presenta un patrón más defenido y ajustado a un variograma Wav, con una semivarianza máxima cercana a 7 y que se estabiliza cerca de los 58000 m, aunque presenta cierta varia

# Se optó por la anisotropía de 90° porque es la dirección que presenta el mejor ajuste.
v_anis_jul <- vgm(8,'Wav',58000,anis = c(90,1))
plot(v_julio_dir,v_anis_jul)

Anisotropía espacial del variograma ondular (Wav)

# Mapa Variográfico
plot(variogram(temp~1,data_t_jul,cutoff = 120000, width= 8000, map = TRUE))

Del mapa se interpreta una mayor variabilidad de los datos en las cordilleras de la Costa y de los Andes.

3.5.1 Ajuste del variograma

# Guardar en .rds
fit_v_jul <- fit.variogram(v_jul,m_var_jul) # Ajuste del variograma teórico (mv) con el variograma experimental (v).

write_rds(fit_v_jul,'C:/Re6/Geoestadistica/Talleres/data/procesada/var_jul_ajustado.rds')

plot(v_jul, fit_v_jul, scales = list(x = list(log = FALSE), y = list(log = FALSE)),
     col = "blue", lwd = 1,
     main = "Variograma ajustado de Temperatura para julio 2023",
     xlab = "Distancia (m)",
     ylab = "Semivarianza")

png( "C:/Re6/Geoestadistica/Talleres/Outputs/Grafico/var_julio_ajustado.png", # Guardamos el gráfico
     width = 25, height = 14, units = "cm", res = 300)  
plot(v_jul, fit_v_jul, scales = list(x = list(log = FALSE), y = list(log = FALSE)),
     col = "blue", lwd = 1,
     main = "Variograma ajustado de Temperatura para julio 2023",
     xlab = "Distancia (m)",
     ylab = "Semivarianza")
dev.off()

4 Septiembre

4.1 Indice de Moran Global

4.1.1 Carga de datos

# Cargar datos
data_t_sept <- read_rds('C:/Re6/Geoestadistica/Talleres/data/procesada/data_temp_septiembre.rds')

# Cálculo de la matriz de distancia

w <- 1/as.matrix(st_distance(data_t_sept))
diag(w) <- 0
w <- units::drop_units(w)
w <- mat2listw(w)
moran(data_t_sept$temp,w,n=length(w$neighbours),S0=Szero(w))

# filtrar datos por distancia entre las estaciones

data_t_sept <- data_t_sept|> 
    distinct(geometry,.keep_all = TRUE) #deja sólo las geometrias únicas

#Volvemos a calcular moran

w <- 1/as.matrix(st_distance(data_t_sept))
diag(w) <- 0
w <- units::drop_units(w)
w <- mat2listw(w)
m <- moran(data_t_sept$temp,w,n=length(w$neighbours),S0=Szero(w))
mp <- moran.plot(data_t_sept$temp,w)

# Buscar las estaciones que se encuentran muy cercanas

d <- st_distance(data_t_sept) #calcula la distance entre estaciones
diag(d) <- NA #convierte la diagonal de la matriz de distancia en NA
d <- units::drop_units(d) # quita las unidades de la distancia
ind <- which(d < 100,arr.ind = TRUE) # busca que estaciones se encuentran a menos de 100m
data_t_septf <- data_t_sept[-ind[,1],] #elimina las estaciones que estan a una distancia de menos de 100m
w <- 1/as.matrix(st_distance(data_t_sept))
diag(w) <- 0
w <- units::drop_units(w)
w <- mat2listw(w)
moran(data_t_sept$temp,w,n=length(w$neighbours),S0=Szero(w))

#Test estadístico Montecarlo para probar la significancia de la autocorrelación espacial

moran.mc(data_t_sept$temp,w,nsim=10000)
mp <- moran.plot(data_t_sept$temp,w) #Con esto se descarta la hipotesis nula.

Al realizar la prueba de Montecarlo para el índice global de Moran se puede afirmar que existe una correlación espacial positiva.
Porque la distribución de los datos tiene un valor estadístico de 0.27261 con un p-value de 0.0001 tras 10000 simulaciones.

4.2 Indice de Moran Local

lmoran <- localmoran(data_t_sept$temp,w,alternative = "greater")
head(lmoran)

4.2.1 Mapa de autocorrelacion

data_t_sept <- st_as_sf(data_t_sept[c('geometry','temp')],coords = c('x','y'),crs = 32719)
data_t_sept$Ii <- lmoran[,'Ii']
data_t_sept$Z <- lmoran[,'Z.Ii']
data_t_sept$Pr <- lmoran[, "Pr(z > E(Ii))"]

# Crear los mapas
p1 <- tm_shape(data_t_sept) +
  tm_dots(col = "temp", title = "Temperatura", style = "quantile", palette = "viridis") +
  tm_layout(legend.position = c("0.9", "0.1"), legend.outside = FALSE, frame = FALSE) +
  tm_scale_bar(position = c("right", "bottom"), text.size = 0.5) +
  tm_grid(proj = 4326, n.x = 8, n.y = 8, lines = TRUE, col = "black", lwd = 0.5, alpha = 0.5)  

# Repite para los otros mapas
p2 <- tm_shape(data_t_sept) +
  tm_dots(col = "Ii", title = "Índice de Moran Local", style = "quantile", palette = "viridis") +
  tm_layout(legend.position = c("0.9", "0.1"), legend.outside = FALSE, frame = FALSE) +
  tm_scale_bar(position = c("right", "bottom"), text.size = 0.5) +
  tm_grid(proj = 4326, n.x = 8, n.y = 8, lines = TRUE, col = "black", lwd = 0.5, alpha = 0.5)  

p3 <- tm_shape(data_t_sept) +
  tm_dots(col = "Z", title = "Z-score", style = "fixed", breaks = c(-Inf, 1.65, Inf), palette = "viridis") +
  tm_layout(legend.position = c("0.9", "0.1"), legend.outside = FALSE, frame = FALSE) +
  tm_scale_bar(position = c("right", "bottom"), text.size = 0.5) +
  tm_grid(proj = 4326, n.x = 8, n.y = 8, lines = TRUE, col = "black", lwd = 0.5, alpha = 0.5)  

p4 <- tm_shape(data_t_sept) +
  tm_dots(col = "Pr", title = "p-value", style = "fixed", breaks = c(-Inf, 0.05, Inf), palette = "viridis") +
  tm_layout(legend.position = c("0.9", "0.1"), legend.outside = FALSE, frame = FALSE) +
  tm_scale_bar(position = c("right", "bottom"), text.size = 0.5) +
  tm_grid(proj = 4326, n.x = 8, n.y = 8, lines = TRUE, col = "black", lwd = 0.5, alpha = 0.5)

tmap_mode('view')
tmap_arrange(p1,p2,p3,p4)

Mapa Interactivo

Mapa Estático

#Quad

quad <- attr(lmoran, 'quadr') 
head(quad)

data_t_sept <- data_t_sept |> 
  bind_cols(quad = quad$median) |>  # Agregando columna quad
  mutate(quad = case_when(
    Pr < 0.05 ~ quad,              # Mantiene valores significativos
    .default = 'No Significativo')) # Asigna "No Significativo" cuando Pr >= 0.05

# Crear el mapa con la columna "quad"
moran <- tm_shape(data_t_sept) +
  tm_dots(col = "quad", title = "Quad", style = "cat", palette = "Set1") +  # Cambiar a 'quad' y estilo 'cat' para categorías
  tm_layout(legend.outside = TRUE)

tmap_mode('view')
moran

El mapa del índice moran Quad nos muestra la ubicación y concentración de puntos de temperatura en autocorrelación espacial a sus vecinos.
Los puntos en rojo (High-High) son aquellos con valores altos de temperatura entorno a la media, rodeados por puntos vecinos con temperaturas similares.Son puntos con autocorrelación espacial positva
Los puntos en azul (High-Low) son aquellos con valores altos de temperatura en torno a la media, rodeados por puntos vecinos con temperaturas más bajas. Son puntos con autocorrelación espacial negativa.
Los puntos en rojo (Low-High) son aquellos con valores bajos de temperatura en torno a la media, rodeados por puntos vecinos con temperaturas más altas.Son puntos con autocorrelación espacial negativa.
Los puntos en azul (Low-Low) son aquellos con valores bajos de temperatura en torno a la media, rodeados por puntos vecinos con temperaturas similares. Son puntos con autocorrelación espacial positiva.
Los puntos en naranjo (No Significativo) son aquellos cuyos valores de temperatura tienen un comportamiento aleatorio en relación a sus vecionos. Son puntos sin autocorrelación espacial.
Al comparar el índice de moran global y el índice de moran local, se pueden afirmar dos cosas principalmente:
La autocorrelación espacial general es positiva. Eso dice el índice de moral global (0.27261), lo que se puede corroborar visualmente con la presencia mayoritaria de puntos rojos (High-High) y morados (Low-Low).
Que los mayores clusters de autocorrelación se encuentran en las cordilleras de la Costa (High-High) y de los Andes (Low-Low).

4.3 Análisis del variograma

4.3.1 Variograma experimental

options(scipen = 10)

# Realizamos un variograma exploratorio para identificar qué variograma se ajusta más a nuestros datos
v <- variogram(temp~1,data_t_sept,cutoff = 90000, width = 6000) 
plot(v)

show.vgms() #Verificar tipos de variogramas teoricos para ver cuál se ajusta más

v_sept <- variogram(temp~1,data_t_sept,cutoff = 30000, width = 2000)  #Definimos el punto de corte.

plot(v_sept)

# Definimos el sill (máximo valor de los datos), el tipo de variograma (Wav), y el nugget (la variabilidad inicial de los datos).
m_var_sept <- vgm(1.35,'Wav',28000,0) # 

plot(v_sept,m_var_sept)

El variograma experimental de Septiembre nos indica que la máxima semiviaranza es de alredor de 1.5 y que se estabiliza a los 28000 m, y se asemaja un variograma de onda (Wav).Con el punto de corte 30000 m se elimaron los puntos que evidencian la estabilización de la semivarianza.

# Guardamos

png( "C:/Re6/Geoestadistica/Talleres/Outputs/Grafico/var_septiembre.png",
     width = 25, height = 14, units = "cm", res = 300)
plot(v_sept,m_var_sept, 
     col = "blue", lwd = 1,
     main = "Variograma de Temperatura para septiembre 2023",
     xlab = "Distancia (m)",
     ylab = "Semivarianza")
dev.off()

write_rds(m_var_sept,'C:/Re6/Geoestadistica/Talleres/data/procesada/var_sept.rds')

4.4 Análisis de Anisotropía

# Definir variogramas en diferentes direcciones
v_septiembre_dir <- variogram(temp ~ 1, data_t_sept, alpha = c(0, 45, 90, 135), cutoff = 150000, width = 10000)
plot(v_septiembre_dir)

El gráfico 0° presenta una curva plana, y una casi nula variabilidad.
El gráfico 45° presenta una máxima semivarianza máxima cercana a 4 que se estabiliza alrededor de los 55000 m.
El gráfico 90° presenta un patrón más defenido, con una semivarianza máxima cercana a 39 y que se estabiliza cerca de los 140000 m.
El gráfico 135° presenta una máxima variabilidad cercana a 10.5 y que se estabiliza alrededor de los 90000 m.

# Se optó por la anisotropía de 90° porque es la dirección que presenta el mejor ajuste.
v_anis_sept <- vgm(39,'Wav',137000,0.2,anis = c(90,1))
plot(v_septiembre_dir,v_anis_sept)

La anisotropía espacial del variograma ondular (Wav)

plot(variogram(temp~1,data_t_sept,cutoff = 120000, width= 8000, map = TRUE))

Del mapa se interpreta una mayor autocorrelación de los datos en las cordilleras de la Costa y de los Andes.

4.4.1 Ajuste del variograma

# Guardar en .rds
fit_v_sept <- fit.variogram(v_sept,m_var_sept) # Ajuste del variograma teórico (mv) con el variograma experimental (v).

write_rds(fit_v_sept,'C:/Re6/Geoestadistica/Talleres/data/procesada/var_sept_ajustado.rds')

plot(v_sept, fit_v_sept, scales = list(x = list(log = FALSE), y = list(log = FALSE)),
     col = "blue", lwd = 1,
     main = "Variograma ajustado de Temperatura para septiembre 2023",
     xlab = "Distancia (m)",
     ylab = "Semivarianza")

png( "C:/Re6/Geoestadistica/Talleres/Outputs/Grafico/var_septiembre_ajustado.png", # Guardamos el gráfico
     width = 25, height = 14, units = "cm", res = 300)  
plot(v_sept, fit_v_sept, scales = list(x = list(log = FALSE), y = list(log = FALSE)),
     col = "blue", lwd = 1,
     main = "Variograma ajustado de Temperatura para septiembre 2023",
     xlab = "Distancia (m)",
     ylab = "Semivarianza")
dev.off()

5 Noviembre

5.1 Indice de Moran Global

5.1.1 Carga de datos

# Cargar datos
data_t_nov <- read_rds('C:/Re6/Geoestadistica/Talleres/data/procesada/data_temp_noviembre.rds')

# Cálculo de la matriz de distancia

w <- 1/as.matrix(st_distance(data_t_nov))
diag(w) <- 0
w <- units::drop_units(w)
w <- mat2listw(w)
moran(data_t_nov$temp,w,n=length(w$neighbours),S0=Szero(w))

# filtrar datos por distancia entre las estaciones

data_t_nov <- data_t_nov|> 
    distinct(geometry,.keep_all = TRUE) #deja sólo las geometrias únicas

#Volvemos a calcular moran

w <- 1/as.matrix(st_distance(data_t_nov))
diag(w) <- 0
w <- units::drop_units(w)
w <- mat2listw(w)
m <- moran(data_t_nov$temp,w,n=length(w$neighbours),S0=Szero(w))
mp <- moran.plot(data_t_nov$temp,w)

#Test estadístico Montecarlo para probar la significancia de la autocorrelación espacial

moran.mc(data_t_nov$temp,w,nsim=10000)

mp <- moran.plot(data_t_nov$temp,w) 
# Con esto se descarta la hipotesis nula. Ya que índice de Moran global tras 10000 simulaciones genero una autocorrelación espacial positiva de 0.26516. 
# El p-value (0.0001) indica que esto ocurre en el 0.9999 de los casos.

5.2 Indice de Moran Local

lmoran <- localmoran(data_t_nov$temp,w,alternative = "greater")
head(lmoran)

5.2.1 Mapa de autocorrelacion

data_t_nov <- st_as_sf(data_t_nov[c('geometry','temp')],coords = c('x','y'),crs = 32719)
data_t_nov$Ii <- lmoran[,'Ii']
data_t_nov$Z <- lmoran[,'Z.Ii']
data_t_nov$Pr <- lmoran[, "Pr(z > E(Ii))"]

# Crear los mapas
p1 <- tm_shape(data_t_nov) +
  tm_dots(col = "temp", title = "Temperatura", style = "quantile", palette = "viridis") +
  tm_layout(legend.position = c("0.9", "0.1"), legend.outside = FALSE, frame = FALSE) +
  tm_scale_bar(position = c("right", "bottom"), text.size = 0.5) +
  tm_grid(proj = 4326, n.x = 8, n.y = 8, lines = TRUE, col = "black", lwd = 0.5, alpha = 0.5)  

# Repite para los otros mapas
p2 <- tm_shape(data_t_nov) +
  tm_dots(col = "Ii", title = "Índice de Moran Local", style = "quantile", palette = "viridis") +
  tm_layout(legend.position = c("0.9", "0.1"), legend.outside = FALSE, frame = FALSE) +
  tm_scale_bar(position = c("right", "bottom"), text.size = 0.5) +
  tm_grid(proj = 4326, n.x = 8, n.y = 8, lines = TRUE, col = "black", lwd = 0.5, alpha = 0.5)  

p3 <- tm_shape(data_t_nov) +
  tm_dots(col = "Z", title = "Z-score", style = "fixed", breaks = c(-Inf, 1.65, Inf), palette = "viridis") +
  tm_layout(legend.position = c("0.9", "0.1"), legend.outside = FALSE, frame = FALSE) +
  tm_scale_bar(position = c("right", "bottom"), text.size = 0.5) +
  tm_grid(proj = 4326, n.x = 8, n.y = 8, lines = TRUE, col = "black", lwd = 0.5, alpha = 0.5)  

p4 <- tm_shape(data_t_nov) +
  tm_dots(col = "Pr", title = "p-value", style = "fixed", breaks = c(-Inf, 0.05, Inf), palette = "viridis") +
  tm_layout(legend.position = c("0.9", "0.1"), legend.outside = FALSE, frame = FALSE) +
  tm_scale_bar(position = c("right", "bottom"), text.size = 0.5) +
  tm_grid(proj = 4326, n.x = 8, n.y = 8, lines = TRUE, col = "black", lwd = 0.5, alpha = 0.5)

tmap_mode('view')
tmap_arrange(p1,p2,p3,p4)

Mapa Interactivo

Mapa Estático

#Quad

quad <- attr(lmoran, 'quadr') 
head(quad)

data_t_nov <- data_t_nov |> 
  bind_cols(quad = quad$median) |>  # Agregando columna quad
  mutate(quad = case_when(
    Pr < 0.05 ~ quad,              # Mantiene valores significativos
    .default = 'No Significativo')) # Asigna "No Significativo" cuando Pr >= 0.05

# Crear el mapa con la columna "quad"
moran <- tm_shape(data_t_nov) +
  tm_dots(col = "quad", title = "Quad", style = "cat", palette = "Set1") +  # Cambiar a 'quad' y estilo 'cat' para categorías
  tm_layout(legend.outside = TRUE)

tmap_mode('view')
moran

El mapa del índice moran Quad nos muestra la ubicación y concentración de puntos de temperatura en autocorrelación espacial a sus vecinos.
Los puntos en rojo (High-High) son aquellos con valores altos de temperatura entorno a la media, rodeados por puntos vecinos con temperaturas similares.Son puntos con autocorrelación espacial positva
Los puntos en azul (High-Low) son aquellos con valores altos de temperatura entorno a la media, rodeados por puntos vecinos con temperaturas más bajas. Son puntos con autocorrelación espacial negativa.
Los puntos en rojo (Low-High) son aquellos con valores bajos de temperatura entorno a la media, rodeados por puntos vecinos con temperaturas más altas.Son puntos con autocorrelación espacial negativa.
Los puntos en azul (Low-Low) son aquellos con valores bajos de temperatura entorno a la media, rodeados por puntos vecinos con temperaturas similares. Son puntos con autocorrelación espacial positiva.
Los puntos en naranjo (No Significativo) son aquellos cuyos valores de temperatura tienen un comportamiento aleatorio en relación a sus vecionos. Son puntos sin autocorrelación espacial.
Al comparar el índice de moran global y el índice de moran local, se pueden afirmar dos cosas principalmente:
La autocorrelación espacial general es positiva. Eso dice el índice de moral global (0.26516), lo que se puede corroborar visualmente con la presencia mayoritaria de puntos rojos (High-High) y morados (Low-Low).
Que los mayores clusters de autocorrelación se encuentran en las cordilleras de la Costa (High-High) y de los Andes (Low-Low).

5.3 Análisis del Variograma

5.3.1 Variograma Experimental

options(scipen = 10)

# Realizamos un variograma exploratorio para identificar qué variograma se ajusta más a nuestros datos
v <- variogram(temp~1,data_t_nov, plot.numbers = TRUE) 
plot(v)


show.vgms() #Verificar tipos de variogramas teoricos para ver cuál se ajusta más

v_nov <- variogram(temp~1,data_t_nov,cutoff = 80000, width = 5500) #Definimos el punto de corte.

plot(v_nov)

# Definimos el sill (máximo valor de los datos), el tipo de variograma (Wav), y el nugget (la variabilidad inicial de los datos).
m_var_nov <- vgm(7.8,'Wav',70000,0.) #

plot(v_nov,m_var_nov)

El variograma experimental de Noviembre nos indica que la máxima semiviaranza es de alredor de 7.8 y que se estabiliza a los 68000 m, y se asemaja un variograma de onda (Wav).Con el punto de corte 80000 m se elimaron los puntos que evidencian la estabilización de la semivarianza.

# Guardamos

png( "C:/Re6/Geoestadistica/Talleres/Outputs/Grafico/var_noviembre.png",
     width = 25, height = 14, units = "cm", res = 300)
plot(v_nov,m_var_nov, 
     col = "blue", lwd = 1,
     main = "Variograma de Temperatura para noviembre 2023",
     xlab = "Distancia (m)",
     ylab = "Semivarianza")
dev.off()

write_rds(m_var_nov,'C:/Re6/Geoestadistica/Talleres/data/procesada/var_nov.rds')

# Definir variogramas en diferentes direcciones
v_noviembre_dir <- variogram(temp ~ 1, data_t_nov, alpha = c(0, 45, 90, 135), cutoff = 120000, width=6000)
plot(v_noviembre_dir)

El gráfico 0° presenta una curva plana, y una casi nula variabilidad.
El gráfico 45° presenta una máxima semivarianza máxima cercana a 3 que se estabiliza al rededor de los 50000 m.
EL gráfico 90° presenta un patrón más defenido y ajustado a un variograma Gausiano, con una semivarianza máxima cercana a 22.5 y que se estabiliza cerca de los 108000 m.
EL gráfico 135° presenta una máxima variabilidad cercana a 7.5 y que se estabiliza al rededor de los 78000 m.

5.4 Análisis de Anisotropía

# Se optó por la anisotropía de 90° porque es la dirección que presenta el mejor ajuste.
v_anis_nov <- vgm(22.5,'Wav',108000,0.2,anis = c(90,1))
plot(v_noviembre_dir,v_anis_nov)

La anisotropía espacial del variograma ondular (Wav)

plot(variogram(temp~1,data_t_nov,cutoff = 120000, width= 8000, map = TRUE))

Del mapa se interpreta una mayor variabilidad de los datos en las cordilleras de la Costa y de los Andes.

5.4.1 Ajuste del variograma

# Guardar en .rds
fit_v_nov <- fit.variogram(v_nov,m_var_nov) # Ajuste del variograma teórico (mv) con el variograma experimental (v).

write_rds(fit_v_nov,'C:/Re6/Geoestadistica/Talleres/data/procesada/var_nov_ajustado.rds')

plot(v_nov, fit_v_nov, scales = list(x = list(log = FALSE), y = list(log = FALSE)),
     col = "blue", lwd = 1,
     main = "Variograma ajustado de Temperatura para noviembre 2023",
     xlab = "Distancia (m)",
     ylab = "Semivarianza")

png( "C:/Re6/Geoestadistica/Talleres/Outputs/Grafico/var_noviembre_ajustado.png", # Guardamos el gráfico
     width = 25, height = 14, units = "cm", res = 300)  
plot(v_nov, fit_v_nov, scales = list(x = list(log = FALSE), y = list(log = FALSE)),
     col = "blue", lwd = 1,
     main = "Variograma ajustado de Temperatura para noviembre 2023",
     xlab = "Distancia (m)",
     ylab = "Semivarianza")
dev.off()

6 Variogramas de residuos de los modelos de regresión

6.1 Febrero

Para el mes de Febrero, el modelo con mejor relación RMSE-Rcuadrado fue el que solamente utilizó al DEM como predictor

vr6_feb <- variogram(temp ~ sqrt(dem), data_unida_feb, cutoff = 40000, width = 6000) # Definiendo los puntos de corte del variograma variograma.

plot(vr6_feb)


var_r6_feb <- vgm(0.8,'Pen',33500, 0.05)

plot(vr6_feb,var_r6_feb)

png( "C:/Re6/Geoestadistica/Talleres/Outputs/Grafico/residual_var_febrero.png",
     width = 25, height = 14, units = "cm", res = 300)
plot(vr6_feb,var_r6_feb, 
     col = "blue", lwd = 1,
     main = "Variograma residual de modelo de Temperatura para Febrero 2023",
     xlab = "Distancia (m)",
     ylab = "Semivarianza")
dev.off()

6.1.1 Variograma ajustado

fvar_r6_feb <- fit.variogram(vr6_feb,var_r6_feb)

write_rds(fvar_r6_feb,'C:/Re6/Geoestadistica/Talleres/data/procesada/residual_var_febrero_ajustado.rds')

plot(vr6_feb, fvar_r6_feb, scales = list(x = list(log = FALSE), y = list(log = FALSE)),
     col = "blue", lwd = 1,
     main = "Variograma residual ajustado de regresión de Temperatura para Febrero 2023",
     xlab = "Distancia (m)",
     ylab = "Semivarianza")

png( "C:/Re6/Geoestadistica/Talleres/Outputs/Grafico/residual_var_febrero_ajustado.png",
     width = 25, height = 14, units = "cm", res = 300)
plot(vr6_feb,var_r6_feb, 
     col = "blue", lwd = 1,
     main = "Variograma residual ajustado de regresión de Temperatura para Febrero 2023",
     xlab = "Distancia (m)",
     ylab = "Semivarianza")
dev.off()

El variograma residual experimental del mes de Febrero nos indica que la máxima semiviaranza es de al rededor de 0.85 y que se estabiliza a los 3900 m, y se asemeja un variograma tipo Pen. Con el punto de corte 40000 m se eliminaron los puntos que evidencian la estabilización de la semivarianza.
El variograma residual ajustado del mes de Febrero nos indica que la máxima semiviaranza es de al rededor de 0.85 y que se estabiliza a los 3900 m, y se asemeja un variograma tipo Pen. Con el punto de corte 40000 m se eliminaron los puntos que evidencian la estabilización de la semivarianza.

6.2 Julio

Para el mes de Julio, el modelo con mejor relación RMSE-Rcuadrado fue el que solamente utilizó al DEM como predictor

vr6_jul <- variogram(temp ~ sqrt(dem), data_unida_jul, cutoff = 40000, width = 5500) # Definiendo los puntos de corte del variograma variograma.

plot(vr6_jul)

var_r6_jul <- vgm(0.95,'Ste',39500, 0.15)

plot(vr6_jul,var_r6_jul)

png( "C:/Re6/Geoestadistica/Talleres/Outputs/Grafico/residual_var_julio.png",
     width = 25, height = 14, units = "cm", res = 300)
plot(vr6_jul,var_r6_jul, 
     col = "blue", lwd = 1,
     main = "Variograma residual de modelo de regresión de Temperatura para Julio 2023",
     xlab = "Distancia (m)",
     ylab = "Semivarianza")
dev.off()

6.2.1 Variograma ajustado

fvar_r6_jul <- fit.variogram(vr6_jul,var_r6_jul)

write_rds(fvar_r6_jul,'C:/Re6/Geoestadistica/Talleres/data/procesada/residual_var_julio_ajustado.rds')

plot(vr6_jul, fvar_r6_jul, scales = list(x = list(log = FALSE), y = list(log = FALSE)),
     col = "blue", lwd = 1,
     main = "Variograma residual ajustado de modelo de regresión de Temperatura para Julio 2023",
     xlab = "Distancia (m)",
     ylab = "Semivarianza")

png( "C:/Re6/Geoestadistica/Talleres/Outputs/Grafico/residual_var_julio_ajustado.png",
     width = 25, height = 14, units = "cm", res = 300)
plot(vr6_jul, fvar_r6_jul, 
     col = "blue", lwd = 1,
     main = "Variograma residual ajustado de modelo de regresión de Temperatura para Julio 2023",
     xlab = "Distancia (m)",
     ylab = "Semivarianza")
dev.off()

El variograma residual experimental del mes de Julio nos indica que la máxima semiviaranza es de al rededor de 0.85 y que se estabiliza a los 3900 m, y se asemeja un variograma tipo Ste. Con el punto de corte 40000 m se eliminaron los puntos que evidencian la estabilización de la semivarianza.
El variograma residual ajustado del mes de Julio nos indica que la máxima semiviaranza es de al rededor de 0.8 y que se estabiliza a los 3900 m, y se asemeja un variograma tipo Ste. Con el punto de corte 40000 m se eliminaron los puntos que evidencian la estabilización de la semivarianza.

6.3 Septiembre

El modelo de regresión con mejor relación RMSE-Rcuadrado para el mes de Septiembre fue el que utilizó como predictores al DEM, la distancia a la costa, el NDVI y el LST.

vr6_sept <- variogram(temp ~ sqrt(dem + dist_costa + ndvi_sept + lst_sept), data_unida_sept, cutoff = 40000, width = 2500) # Definiendo los puntos de corte del variograma variograma.

plot(vr6_sept)

var_r6_sept <- vgm(3,'Wav',39500, 0)

plot(vr6_sept,var_r6_sept)

png( "C:/Re6/Geoestadistica/Talleres/Outputs/Grafico/residual_var_septiembre.png",
     width = 25, height = 14, units = "cm", res = 300)
plot(vr6_sept, var_r6_sept, 
     col = "blue", lwd = 1,
     main = "Variograma residual de modelo de regresión de Temperatura para Septiembre 2023",
     xlab = "Distancia (m)",
     ylab = "Semivarianza")
dev.off()

6.3.1 Variograma ajustado

fvar_r6_sept <- fit.variogram(vr6_sept,var_r6_sept)

write_rds(fvar_r6_sept,'C:/Re6/Geoestadistica/Talleres/data/procesada/residual_var_septiembre_ajustado.rds')

plot(vr6_sept, fvar_r6_sept, scales = list(x = list(log = FALSE), y = list(log = FALSE)),
     col = "blue", lwd = 1,
     main = "Variograma residual ajustado de modelo de regresión de Temperatura para Septiembre 2023",
     xlab = "Distancia (m)",
     ylab = "Semivarianza")

png( "C:/Re6/Geoestadistica/Talleres/Outputs/Grafico/residual_var_septiembre_ajustado.png",
     width = 25, height = 14, units = "cm", res = 300)
plot(vr6_sept, fvar_r6_sept, 
     col = "blue", lwd = 1,
     main = "Variograma residual ajustado de modelo de regresión de Temperatura para Septiembre 2023",
     xlab = "Distancia (m)",
     ylab = "Semivarianza")
dev.off()

El variograma residual experimental del mes de Septiembre nos indica que la máxima semiviaranza es de al rededor de 3 y que se estabiliza a los 3900 m, y se asemeja un variograma tipo Wav. Con el punto de corte 40000 m se eliminaron los puntos que evidencian la estabilización de la semivarianza.
El variograma residual ajustado del mes de Septiembre nos indica que la máxima semiviaranza es de al rededor de 2.9 y que se estabiliza a los 3900 m, y se asemeja un variograma tipo Wav. Con el punto de corte 40000 m se eliminaron los puntos que evidencian la estabilización de la semivarianza.

6.4 Noviembre

El modelo de regresión con mejor relación RMSE-Rcuadrado para el mes de Noviembre fue el que utilizó al DEM y al LST como predictores.

vr6_nov <- variogram(temp ~ sqrt(dem + lst_nov), data_unida_nov, cutoff = 40000, width = 5000) # Definiendo los puntos de corte del variograma variograma.

plot(vr6_nov)

var_r6_nov <- vgm(0.62,'Ste',38500, 0.16)

plot(vr6_nov,var_r6_nov)

png( "C:/Re6/Geoestadistica/Talleres/Outputs/Grafico/residual_var_noviembre.png",
     width = 25, height = 14, units = "cm", res = 300)
plot(vr6_nov,var_r6_nov, 
     col = "blue", lwd = 1,
     main = "Variograma residual de modelo de regresión de Temperatura para Noviembre 2023",
     xlab = "Distancia (m)",
     ylab = "Semivarianza")
dev.off()

6.4.1 Variograma ajustado

fvar_r6_nov <- fit.variogram(vr6_nov,var_r6_nov)

write_rds(fvar_r6_nov,'C:/Re6/Geoestadistica/Talleres/data/procesada/residual_var_noviembre_ajustado.rds')

plot(vr6_nov, fvar_r6_nov, scales = list(x = list(log = FALSE), y = list(log = FALSE)),
     col = "blue", lwd = 1,
     main = "Variograma residual ajustado de modelo de regresión de Temperatura para Noviembre 2023",
     xlab = "Distancia (m)",
     ylab = "Semivarianza")

png( "C:/Re6/Geoestadistica/Talleres/Outputs/Grafico/residual_var_septiembre_ajustado.png",
     width = 25, height = 14, units = "cm", res = 300)
plot(vr6_nov, fvar_r6_nov, 
     col = "blue", lwd = 1,
     main = "Variograma residual ajustado de modelo de regresión de Temperatura para Noviembre 2023",
     xlab = "Distancia (m)",
     ylab = "Semivarianza")
dev.off()

El variograma residual experimental del mes de Noviembre nos indica que la máxima semiviaranza es de al rededor de 0.62 y que se estabiliza a los 3900 m, y se asemeja un variograma tipo Ste. Con el punto de corte 40000 m se eliminaron los puntos que evidencian la estabilización de la semivarianza.
El variograma residual ajustado del mes de Noviembre nos indica que la máxima semiviaranza es de al rededor de 0.58 y que se estabiliza a los 3900 m, y se asemeja un variograma tipo Ste. Con el punto de corte 40000 m se eliminaron los puntos que evidencian la estabilización de la semivarianza.