Taller 3 Región Metropolitana

Author

José Bórquez & Sergio Barrera

Published

September 27, 2023

Magister en Teledetección. Análisis de Datos Espaciales con R. Análisis de Autocorrelación Espacial

Introducción:

Para comprender bien este taller número tres hay que referirse en primera instancia al concepto de autocorrelación espacial. Este cuantifica la falta de aleatoriedad de los valores de una variable, debido a su estructura espacial. Bajo este concepto se espera que valores pequeños se encuentren cercanos de otros valores pequeños, mientras que valores altos lo hagan en proximidades a otros valores altos.

Para establecer la existencia y significancia de una posible autocorrelación espacial dentro de nuestro conjunto de datos de temperatura se utilizaron 2 medidas estadísticas diferentes. En primer lugar se calculó el índice de Moran de forma global y local, y luego se utilizó la herramienta del variograma para describir la dependencia espacial entre las diferentes estaciones con los datos de interés.

Por último, luego de realizar los ajustes pertinentes a los distintos variogramas experimentales generados se consideraron los resultados del taller anterior para obtener un modelo de residuos. Este modelo fue analizado, ajustado e interpretado. Cada una de estas operaciones fueron efectuadas los meses trabajados en el taller n°2.

00. Instalación de los paquetes requeridos:

#install.packages(c('tidyverse','sf','terra','stars','tmap','remotes','geodata','rstac','earthdatalogin','viridis','gt','quarto','tinytex', 'knitr', 'rmarkdown','gstat','gridExtra','spdep','automap'))

01. Carga de las librerías requeridas:

library(remotes)          # Para traer función para instalar paquetes desde github
library(tidyverse)        # Para manipulación y visualización de datos
library(sf)               # Para manejar datos espaciales
library(terra)            # Para manipulación de datos raster
library(tmap)             # Para crear mapas temáticos
library(geodata)          # Para manejar datos geográficos
library(rstac)            # Para trabajar con datos STAC
library(earthdatalogin)   # Para autenticación con Earthdata
library(viridis)          # Para paletas de colores
library(gt)               # Para crear tablas
library(quarto)           # Para crear documentos dinámicos
library(tinytex)          # Para la gestión de LaTeX
library(knitr)            # Para la generación de informes
library(rmarkdown)        # Para documentos de R Markdown
library(gstat)            # Para geostatística y modelado de variogramas
library(stars)            # Para manejar datos en formato stars
library(gridExtra)        # Para organizar gráficos en cuadrículas
library(spdep)            # Para análisis de dependencia espacial
library(automap)          # Para ajuste automático de variogramas

02. Preparación de los datos:

La preparación de los datos para realizar los distintos análisis de correlación espacial constaron de cargar los datos previamente generados, visualizar la ubicación de las estaciones de toma de temperatura reales con las ficticias y preparar una grilla para las interpolaciones.

# Cargar los datos necesarios para el análisis

region <- st_read('data3/region_RM_UTM.gpkg')  # Cargar el archivo geoespacial de la Región Metropolitana en formato GPKG
Reading layer `region_RM_UTM' from data source 
  `C:\METRO\02_metro\data3\region_RM_UTM.gpkg' using driver `GPKG'
Simple feature collection with 1 feature and 11 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 249135.3 ymin: 6205240 xmax: 428386.9 ymax: 6355827
Projected CRS: WGS 84 / UTM zone 19S
data_temp <- read_rds('data3/data_estas_temp_con_predictores.rds')  # Cargar los datos de temperatura con predictores
preds <- rast('data3/predictores.tif')  # Cargar los datos raster de predictores (MDE, orientación, pendiente, etc.)

# Agregar una columna que clasifica las estaciones en ficticias o reales según el station_id

data_temp <- data_temp |> 
  mutate(estacion_ficticia = ifelse(station_id >= 1000, "Ficticia", "Real"))

# Visualizar la ubicación de las estaciones de temperatura (reales y ficticias)

tmap_mode('view')  # Configurar el modo interactivo para visualización en tmap
tm_shape(region) +  # Mostrar la región
  tm_borders() +  # Dibujar los bordes de la región
  tm_shape(data_temp) +  # Añadir las estaciones de temperatura
  tm_dots(col = "estacion_ficticia", size = 0.05, 
          palette = c("blue", "red"),  # Colores para estaciones reales y ficticias
          title = "Tipo de Estación") +  # Leyenda que indica el tipo de estación
  tm_layout(legend.outside = TRUE,       # Leyenda colocada fuera del mapa
            legend.title.size = 0.8,     # Tamaño del título de la leyenda
            legend.text.size = 0.6)      # Tamaño del texto de la leyenda
# Preparar una grilla para la interpolación

bb <- st_bbox(region)  # Obtener los límites geográficos de la región
grilla <- rast(xmin = bb$xmin, xmax = bb$xmax, ymin = bb$ymin, ymax = bb$ymax, 
               res = 500, crs = "EPSG:32719")  # Crear una grilla con resolución de 500 metros y sistema de coordenadas UTM
grilla_st <- grilla |> st_as_stars()  # Convertir la grilla en formato 'stars'

03. Índice de Moran Global para analizar la temperatura de los meses trabajados:

El Índice de Moran, de manera general, determina si una variable se distribuye en el espacio de forma agrupada, aleatoria o dispersa. En el caso del Índice de Moran Global este mide la autocorrelación espacial de un conjunto de datos de forma completa. Entrega una medida única de la tendencia general de los datos, evaluando en nuestro caso si la temperatura de la Región Metropolitana tiene algún tipo de patrón a nivel general.

# Inicializar listas para almacenar resultados

moran_global_results <- list()
moran_mc_results <- list()
moran_plots <- list()

# Meses a analizar

meses <- c(3, 6, 9, 12)  # (Marzo, Junio, Septiembre y Diciembre)

# Realizar el análisis para cada mes

for (mes in meses) {
  # Filtrar los datos de temperatura para el mes correspondiente
  data_temp_filtered <- data_temp %>% filter(month == mes)
  
  # Verificar la existencia de estaciones duplicadas (mismo punto espacial)
  
  if (any(duplicated(data_temp_filtered$geom))) {
    warning(paste("Existen geometrías duplicadas en los datos filtrados para el mes", mes))  # Advertencia si hay duplicados
  }
  
  # Calcular la matriz de distancias entre las estaciones
  
  dist_matrix <- as.matrix(st_distance(data_temp_filtered))  # Generar la matriz de distancias entre estaciones
  w <- 1 / dist_matrix  # Calcular los pesos inversos para las distancias
  diag(w) <- 0  # Poner a cero los elementos diagonales (distancia consigo mismo)
  w <- units::drop_units(w)  # Quitar las unidades de la matriz
  w <- mat2listw(w, style = "W")  # Convertir a formato de lista de vecinos para análisis espacial y especificar el estilo
  
  # Calcular el índice de Moran Global para detectar autocorrelación espacial en la temperatura
  
  moran_global <- moran(data_temp_filtered$temp_promedio, w, n = length(w$neighbours), S0 = Szero(w))
  moran_global_results[[as.character(mes)]] <- moran_global  # Guardar el resultado del índice de Moran Global
  
  # Aplicar la prueba de Monte Carlo para evaluar la significancia de la autocorrelación espacial
  
  moran_mc_result <- moran.mc(data_temp_filtered$temp_promedio, w, nsim = 10000)  # Test Montecarlo con 10,000 simulaciones
  
  # Extraer los resultados relevantes del test Montecarlo
  
  moran_mc_summary <- data.frame(
    statistic = moran_mc_result$statistic,
    p.value = moran_mc_result$p.value,
    month = mes
  )
  
  moran_mc_results[[as.character(mes)]] <- moran_mc_summary  # Guardar los resultados de la prueba
    
  # Graficar el Moran Plot (visualización de la relación espacial)
  
  plot <- moran.plot(data_temp_filtered$temp_promedio, w, main = paste("Moran Plot para el mes de", mes), xlab = "Temperatura Promedio", ylab = "Espacio")
  moran_plots[[as.character(mes)]] <- plot  # Almacenar gráfico
}

Las figuras presentadas anteriormente muestran gráficamente los resultados del Índice de Moran Global donde se puede apreciar gráficamente la correlación espacial positiva existente en la región para todos los meses. En estas representaciones se identifica cómo, en su mayoría, valores altos se relacionan con valores altos y valores bajos se relacionan con valores bajos. En cada uno de los plots se observa como los valores altos se encuentran más concentrados en relación a los valores bajos. Estas representaciones entregan una primera impresión de las métricas que se deberían obtener a nivel numérico para cada mes.

# Convertir los resultados en una tabla y guardarla en un archivo RDS

moran_global_table <- bind_rows(lapply(moran_global_results, as.data.frame), .id = "mes")
print(moran_global_table)
  mes         I        K
1   3 0.3797594 1.574122
2   6 0.4219349 1.344234
3   9 0.4183261 1.305368
4  12 0.4047261 1.428123
write_rds(moran_global_table, "data3/procesada/moran_global_results.rds")  # Guardar en carpeta procesada

En terminos simples para interpretrar el Índice de Moran Global hay que tener en consideración que valores cercanos a 1 indican una autocorrelación positiva, valores cercanos a -1 una autocorrelación negativa y valores cercanos a 0 advierten de ausencia de autocorrelación espacial. Ahora bien, al examinar los resultados del Índice de Moran Global se observa que en todos los meses existe una autocorrelación espacial positiva (como se observó en los Plots mostrados con anterioridad). Lo que sugiere que las estaciones con temperaturas similares tienden a estar cercanas geográficamente. Los valores de este índice para los meses considerados van de 0.37 (mes de marzo) a 0.42 (mes de junio). A su vez, septiembre y diciembre presentan valores de 0.41 y 0.40, respectivamente. Indicando una autocorrelación espacial significativa.

# Convertir los resultados del método de Montecarlo en una tabla y guardarla en un archivo RDS

moran_mc_table <- bind_rows(moran_mc_results)  # No es necesario usar lapply aquí, ya que ya es un data.frame
print(moran_mc_table)
              statistic   p.value month
statistic...1 0.3797594 9.999e-05     3
statistic...2 0.4219349 9.999e-05     6
statistic...3 0.4183261 9.999e-05     9
statistic...4 0.4047261 9.999e-05    12
write_rds(moran_mc_table, "data3/procesada/moran_mc_results.rds")  # Guardar en carpeta procesada

Para conocer la significancia estadística de esta autocorrelación espacial positiva se aplicó la prueba de Montecarlo que evalúa este parámetro mediante simulaciones y probabilidades. En este caso los valores de p.value para todos los meses indican que la autocorrelación obtenida es estadísticamente significativa, por lo tanto, se considera la hipótesis alternativa que indica que efectivamente existe dicha autocorrelación espacial positiva significativa.

En definitiva, el Índice de Moran Global está dando los primeros indicios de que los datos de temperatura disponibles para la Región Metropolitana, en su conjunto, están correlacionados de manera importante a lo largo de todas las estaciones del año. Además, la aplicación de este índice está revelando como los factores geográficos deben ser considerados a la hora de planificar cualquier estudio enfocado en la distribución de la temperatura en la región.

04. Índice de Moran Local para analizar la temperatura de los meses trabajados:

El Índice de Moran Local (Local Indicators of Spatial Association) a diferencia del índice global permite diferenciar zonas geográficas específicas dentro de una misma región, ya que es capaz de detectar focos locales de autocorrelación espacial. En otras palabras, es capaz de identificar grupos de valores similares o distintos en un área definida y determinar cómo se comporta una observación respecto a su entorno.

Para lograr establecer estos grupos de valores se calcula el índice para cada una de las estaciones considerando dos criterios, si el valor de la estación calculada es mayor o menor a la media y si los valores de las estaciones vecinas son mayores o menores a la estación en cuestión. Esto se realiza para cada punto de medición. De estos resultados cada estación se clasifica en uno de los siguientes cuadrantes: high-high, high-low, low-low o low-high. Por ejemplo, para el caso del cuadrante high-high se estaría en presencia de un área de la región donde los valores de las estaciones son altos y están rodeados de otros valores altos, formando clusters de valores altos-altos.También existen los casos de los denominados outliers, o valores anómalos, que ocurren cuando se conjugan valores de Moran Local cercanos a 0 y p.values pequeños. En dichos casos los valores de temperatura serán muy diferentes al de sus vecinos.

#|warning: false
#|message: false

# Inicializar lista para almacenar resultados del índice de Moran local

local_moran_results <- list()

# Realizar el análisis para cada mes
for (mes in meses) {
  # Filtrar los datos de temperatura para el mes correspondiente
  data_temp_filtered <- data_temp %>% filter(month == mes)
  
  # Calcular la matriz de distancias entre las estaciones
  dist_matrix <- as.matrix(st_distance(data_temp_filtered))  # Generar la matriz de distancias entre estaciones
  w <- 1 / dist_matrix  # Calcular los pesos inversos para las distancias
  diag(w) <- 0  # Poner a cero los elementos diagonales (distancia consigo mismo)
  w <- units::drop_units(w)  # Quitar las unidades de la matriz
  w <- mat2listw(w, style = "W")  # Convertir a formato de lista de vecinos para análisis espacial y especificar el estilo
  
  # Calcular el índice de Moran Local
  lm <- localmoran(data_temp_filtered$temp_promedio, w)  # Cálculo del índice de Moran Local
  
  quad <- attr(lm,'quadr') 
  head(quad)
  P_val <- lm[,5]
  
  data_temp_filtered <- data_temp_filtered |> 
    bind_cols(quad = quad$median) |>
    mutate(quad = case_when(
      P_val < 0.05 ~ quad,
      .default = 'No Significativo')) 
  
  # Guardar resultados en la lista
  
  local_moran_results[[as.character(mes)]] <- data_temp_filtered  # Guardar resultados para el mes actual

tmap_mode('plot')
local_moran_map <- tm_shape(region) +  # Añadir capa de fondo (límites o fondo geográfico)
  tm_borders() +  # Dibujar las fronteras de la capa de fondo
  tm_shape(data_temp_filtered) +  # Añadir la capa de puntos
  tm_dots(col = "quad", style = "cat", 
          palette = c('High-High'="red", 'Low-Low'= "blue", 'High-Low'= "green", 'Low-High'="yellow", 'No Significativo'="black"),  
          title = "Cuadrantes Índice de M.L.",
          size = 0.4) +  # Aumentar el tamaño de los puntos
  tm_layout(title = paste("Índice de Moran Local - Temperatura Promedio -", mes),  # Título del mapa
            title.size = 2.5,  # Aumentar el tamaño del título del mapa
            legend.outside = TRUE,  # Colocar la leyenda fuera del mapa
            legend.title.size = 0.9,  # Tamaño del título de la leyenda
            legend.text.size = 0.9)  # Tamaño del texto de la leyenda

# Mostrar el mapa
print(local_moran_map)


}
tmap mode set to plotting

tmap mode set to plotting

tmap mode set to plotting

tmap mode set to plotting

Para el caso del mes de marzo las estaciones totales se clasifican en su gran mayoría en dos grupos de cuadrantes que son del tipo low-low, ubicadas principalmente en las zonas cordilleranas de las comunas de San José de Maipo y Lo Barnechea. Y de clase high-high, las que se ubican en la zona sur-poniente de la región caracterizada por presentar comunas de un marcado carácter rural.

Una única estación ubicada entre las comunas de Colina y Til-Til fue clasificada dentro de un cuadrante high-low, lo que está indicando que es un valor anómalo espacialmente, contrastando con el valor de sus vecinos. También existe el caso de un par de estaciones catalogadas dentro del cuadrante low-high lasque están indicando que se tratan de sectores con valores bajos de temperatura que se encuentran rodeados de sectores con valores altos de temperatura. Los cuadrantes con valores no significativos se ubican en las comunas de San José de Maipo, Pirque, Paine y Buin y uno solitario en la comuna de Colina.

Para el caso del mes de junio solo existen cuadrantes de los tipos high-high (valores altos de temperatura con valores altos de temperatura) y low-low (valores bajos de temperatura con valores bajos de temperatura), ubicados los primeros en la depresión intermedia y en la Cordillera de la Costa, mientras que los segundos se localizan más cerca de la Cordillera de Los Andes. Al listado de estaciones con valores no significativos se suman dos estaciones más.

Para el mes de Septiembre el panorama general de la región, en cuanto a la distribución espacial de los valores de temperatura se mantuvo. Pero existen cuadrantes que están indicando comportamientos distintos en comparación con lo que va corrido del año 2023. El cuadrante más significativo en este sentido corresponde al que se encuentra ubicado en la zona urbana de Santiago que paso de un valor high-high a uno de tipo low-high lo que estaría dando cuenta de una baja en las temperaturas de la zona urbana en contraste a los sectores de carácter más periurbanos del Gran Santiago. El análisis se complica también con la cercanía de valores no significativos. Un segundo valor anómalo se encuentra en la comuna de Colina, que es un sector que muestra variaciones constantes en sus valores de temperatura. Situación que también se pudo apreciar en el mes de diciembre.

Para el último mes del año vuelve a aparecer un cuadrante con un área “outlier” espacial en la mencionada comuna de Colina. Pero gracias al cálculo del Índice de Moran Local fue posible conocer el comportamiento de la temperatura a lo largo de todas las estaciones del año, junto con las zonas específicas donde estos valores, ya sean bajos o altos, se localizaban dentro de la Región Metropolitana.

En definitiva, valores del tipo high-high y low-low sugieren agrupamientos espaciales, donde los valores son similares a los de sus vecinos. Por el contrario, las situaciones cruzadas (high-low y low-high) indican valores atípicos, ósea valores que difieren de manera importante sus vecinos. Siempre es importante tener en consideración que en esta clase de cuadrantes puede haber ocurrido un error en la medición o algún error de otro tipo. Por lo que sería recomendable observar los resultados de esas estaciones con un mayor detenimiento. Aunque ese tipo de consideraciones escapan al objetivo de este tipo de ejercicios prácticos.

# Convertir resultados a un data frame para guardar

local_moran_df <- bind_rows(local_moran_results, .id = "mes")  # Asegurarse de que todos los resultados se guardan

# Guardar los resultados en un archivo RDS

write_rds(local_moran_df, "data3/procesada/local_moran_results.rds")  # Ajusta la ruta según tu estructura de carpetas

# Mostrar los resultados

#| results: hideprint(local_moran_df)

Con el Índice de Moran Global únicamente se puede conocer la existencia de la correlación espacial en un área determinada, lo que permite inferir que entre el total de datos de temperatura de los cuales se dispone existiría una correlación geográfica que se traduce en que valores similares que se ubican en el espacio de manera ordenada y no aleatoria. Pero el índice local va un paso más allá y calcula la relación de cada punto con sus vecinos lo que permite identificar y analizar patrones espaciales escondidos dentro de los datos.

El índice global sirve como una primera aproximación a la definición de la autocorrelación espacial de un espacio geográfico, pero el índice local afina esta correlación puesto que el resultado se obtiene para cada punto, lo que permite diferenciar sectores y comportamientos específicos dentro de una región determinada.

05. Análisis de autocorrelación espacial mediante el uso de variogramas para todos los meses trabajados:

El variograma es una función que describe la dependencia espacial de una variable regionalizada. Se utiliza principalmente para modelar la continuidad espacial de dicha variable. En la práctica es una representación gráfica de la semivarianza. Muestra una dispersión de puntos en el cual se compara la distancia entre pares de puntos (eje x) con la semivarianza (eje y). Estos puntos marcan una tendencia que puede ser modelada con una función, esta función representa la curva que mejor se ajusta a esos puntos. Al primer variograma que se obtiene se le conoce como experimental, mientras que al que se ajusta con la función como modelado.Un variograma típico contiene tres componentes clave: el nugget, el sill, y el range

El nugget (también llamado efecto pepita) es una medida de la componente aleatoria de los datos de una muestra, la componente aleatoria representa lo inesperado. Cuando el efecto nugget es muy alto quiere decir que la componente aleatoria es muy fuerte. Reflejando problemas conceptuales o de metodologías a la hora de abordar los datos.

El sill o meseta representa el punto donde se encuentra la máxima variabilidad que puede observarse en un conjunto de datos. En este punto la línea de tendencia deja de crecer, esto implica que ya no hay relación o influencia de un punto con respecto de otro. También puede ocurrir que el sill sea la posición donde la semivarianza deja de aumentar considerablemente con respecto a la distancia.

El range está muy relacionado con la meseta, puesto que es la distancia a partir de la cual los valores de temperatura ya no están correlacionados geográficamente. En concreto en el punto en el eje x en el que el variograma alcanza la meseta.

5.1 Análisis variográfico para el mes de marzo

Análisis y descripción variográfica isotrópica para el mes de marzo

Un variograma isotrópico es una herramienta estadística usada en geoestadística para analizar la variabilidad espacial de una variable en función de la distancia. En el caso isotrópico, se asume que la variabilidad espacial es la misma en todas las direcciones, es decir, no hay anisotropía o direccionalidad en la dependencia espacial. Esto significa que el variograma solo depende de la distancia entre dos puntos, sin importar en qué dirección se mida esa distancia.

# Filtrar los datos para el mes de marzo
data_temp_filtered_marzo <- data_temp[data_temp$month == 3, ]  # Filtrar datos para marzo

# Crear el variograma experimental de las temperaturas
v_marzo <- variogram(temp_promedio ~ 1, data_temp_filtered_marzo, 
                     cutoff = 160000, width = 6000)  # Variograma experimental

# Graficar el variograma experimental
plot(v_marzo, main = "Variograma Experimental - Mes: Marzo", 
     xlab = "Distancia", ylab = "Semivarianza")  # Visualizar el variograma

# Ajustar el variograma usando autoajuste

auto_fit_marzo <- autofitVariogram(temp_promedio ~ 1, data_temp_filtered_marzo)  # Autoajuste del variograma

#| results: hide print(auto_fit_marzo)  # Mostrar los resultados del ajuste

plot(auto_fit_marzo)  # Graficar el modelo ajustado

# Definir manualmente los parámetros iniciales del modelo Gaussiano
nugget_gau_marzo <- 1.1  # Estimación inicial del nugget
sill_gau_marzo <- 79  # Estimación inicial del sill
psill_gau_marzo <- sill_gau_marzo - nugget_gau_marzo  # Cálculo del psill
range_gau_marzo <- 160000  # Estimación del rango

# Crear el modelo de variograma Gaussiano
model_vgm_gau_marzo <- vgm(psill = psill_gau_marzo, model = "Gau",range = range_gau_marzo, nugget = nugget_gau_marzo)  # Modelo Gaussiano

# Ajustar el modelo al variograma experimental
fit_v_gau_marzo <- fit.variogram(v_marzo, model_vgm_gau_marzo)  
print(fit_v_gau_marzo)  # Ver resultados del ajuste
  model      psill    range
1   Nug  0.9805979      0.0
2   Gau 72.4668290 101298.5
# Graficar el variograma ajustado
plot(v_marzo, fit_v_gau_marzo, main = "Variograma Ajustado - Mes: Marzo", 
     xlab = "Distancia", ylab = "Semivarianza")  # Graficar el variograma ajustado

# Guardar el modelo ajustado en un archivo .rds
saveRDS(fit_v_gau_marzo, "data3/procesada/fit_variograma_gau_marzo.rds")  

Para el caso del mes de marzo al observar el variograma experimental se aprecia como antes de los 50 kilómetros de distancia los valores de variabilidad entre las estaciones de temperatura son relativamente bajos. Esto se explica debido a que, a distancias cortas, los valores de temperatura deberían estar más correlacionados, generando una semivarianza baja. A partir de los 75 kilómetros se observa como la distancia va aumentando de manera significativa disminuyendo la correlación entre los datos y aumentante la semivarianza.

Una vez que se obtiene el variograma experimental es momento de ajustarlo utilizando algunos de los distintos modelos disponibles. En esta oportunidad se utilizó el modelo Gaussiano en el cual a distancias cercanas al origen la variabilidad aumenta de manera lenta, aumentando rápidamente a medida que la distancia comienza a crecer.

Considerando los componentes claves del variograma para el caso de marzo se obtuvo un valor de nugget bastante bajo (cercano a 1), que es un muy buen resultado, puesto que la componente aleatoria del modelo es baja. Para el caso del sill se obtuvo que el valor donde ya no existe influencia entre las estaciones se produce cuando la semivarianza alcanza el valor de 72. La situación del range para este mes en específico mostró que pasado los 100 kilómetros la semivarianza deja de aumentar respecto a la distancia de manera significativa.

Análisis y descripción variográfica anisotrópica para las cuatro direcciones principales del mes de marzo

# Crear el variograma experimental anisotrópico en diferentes direcciones
variograma_anisotropico_marzo <- variogram(temp_promedio ~ 1, data_temp_filtered_marzo, 
                                           alpha = c(0, 45, 90, 135),  # Direcciones a analizar
                                           cutoff = 160000, width = 6000, 
                                           tol.hor = 22.5)  # Tolerancia horizontal en grados

# Graficar el variograma experimental anisotrópico
plot(variograma_anisotropico_marzo, main = "Variograma Anisotrópico - Mes: Marzo", 
     xlab = "Distancia", ylab = "Semivarianza")  # Graficar variograma anisotrópico

# Definir modelo anisotrópico
model_vgm_aniso_marzo <- vgm(psill = psill_gau_marzo, model = "Gau", 
                             range = range_gau_marzo, nugget = nugget_gau_marzo, 
                             anis = c(45, 0.5))  # Ajuste de anisotropía

# Ajustar el modelo anisotrópico al variograma experimental
fit_v_aniso_marzo <- fit.variogram(variograma_anisotropico_marzo, model_vgm_aniso_marzo)
print(fit_v_aniso_marzo)  # Mostrar el modelo ajustado
  model     psill    range ang1 anis1
1   Nug  1.005853      0.0    0   1.0
2   Gau 73.177926 161456.2   45   0.5
# Graficar el variograma ajustado anisotrópico
plot(variograma_anisotropico_marzo, fit_v_aniso_marzo, 
     main = "Variograma Ajustado Anisotrópico - Mes: Marzo", 
     xlab = "Distancia", ylab = "Semivarianza")  # Graficar

# Guardar el modelo ajustado anisotrópico
saveRDS(fit_v_aniso_marzo, "data3/procesada/fit_variograma_aniso_marzo.rds")  # Guardar en un archivo .rds

Para continuar con el análisis variográfico se realizó un modelo que consideró los diferentes azimuts en las 4 direcciones principales. Esto con el fin de conocer el comportamiento de la temperatura tanto en dirección norte-sur (0°), este-oeste (90°) como en sus direcciones intermedias (45° y 135°).

Lo primero a considerar del resultado del variograma experimental direccional es que existe una dirección que presenta un comportamiento totalmente diferente al resto, este es el caso de los 0° que tiene una muy baja variabilidad a medida que se aumenta la distancia. La dirección este-oeste a distancias cortas muestra una semivarianza de magnitud 20 que apartir de los 5 kilometros aumenta de manera considerable estableciendose alrrededor del valor de variabiliada de magnitud 50. Los variogramas de las direciones intermedias son los que presentan una mayor variabiliad espacial para el mes de marzo, a diferencia de los variogramas de diecciones norte-sur y este-oeste. La dirección de 135° presentó los valores más altos de variabilidad.

Lo descrito en el párrafo anterior es un indicativo de la presencia de anisotropía en la conducta de la variable de temperatura al interior de la Región Metropolitana. El concepto de anisotropía se refiere a la presencia de diferentes comportamientos en función de las distintas direcciones. En otras palabras, la direccionalidad de los datos es un factor a tener en cuenta a la hora de construir un modelo de temperatura. La contraparte del modelo anisotrópico es el isotrópico. A una variable se le considera isotrópica si su valor es independiente de la dirección.

El modelo que se utilizó para ajustar los distintos variogarmas experimentales fue el de tipo Gaussiano que se adaptó mejor a la dirección de los 135°, seguida a la de los 90°. Este modelo entrego un sill de 73 y un range de 162 kilometros.

5.2 Análisis variográfico para el mes de junio

Análisis y descripción variográfica isotrópica para el mes de junio

# Filtrar los datos para el mes de junio
data_temp_filtered_junio <- data_temp[data_temp$month == 6, ]  # Filtrar datos para junio

# Crear el variograma experimental de las temperaturas
v_junio <- variogram(temp_promedio ~ 1, data_temp_filtered_junio, 
                     cutoff = 160000, width = 6000)  # Variograma experimental

# Graficar el variograma experimental
plot(v_junio, main = "Variograma Experimental - Mes: Junio", 
     xlab = "Distancia", ylab = "Semivarianza")  # Visualizar el variograma

# Ajustar el variograma usando autoajuste

auto_fit_junio <- autofitVariogram(temp_promedio ~ 1, data_temp_filtered_junio)  # Autoajuste del variograma

#| results: hide print(auto_fit_junio)  # Mostrar los resultados del ajuste

plot(auto_fit_junio)  # Graficar el modelo ajustado

# Definir manualmente los parámetros iniciales del modelo Ste (Stein)
nugget_ste_junio <- 0.4  # Estimación inicial del nugget
sill_ste_junio <- 180  # Estimación inicial del sill
psill_ste_junio <- sill_ste_junio - nugget_ste_junio  # Cálculo del psill
range_ste_junio <- 160000  # Estimación del rango

# Crear el modelo de variograma Stein (Ste)
model_vgm_ste_junio <- vgm(psill = psill_ste_junio, model = "Ste",range = range_ste_junio, nugget = nugget_ste_junio, kappa = 1.3) 

# Ajustar el modelo al variograma experimental
fit_v_ste_junio <- fit.variogram(v_junio, model_vgm_ste_junio)  
print(fit_v_ste_junio)  # Ver resultados del ajuste
  model    psill    range kappa
1   Nug   0.0000      0.0   0.0
2   Ste 704.2124 390139.6   1.3
# Graficar el variograma ajustado
plot(v_junio, fit_v_ste_junio, main = "Variograma Ajustado - Mes: Junio", 
     xlab = "Distancia", ylab = "Semivarianza")  # Graficar el variograma ajustado

# Guardar el modelo ajustado en un archivo .rds
saveRDS(fit_v_ste_junio, "data3/procesada/fit_variograma_ste_junio.rds")  

Para el mes de junio el variograma experimental se aprecia bastante parecido al de marzo hasta más menos los 75 kilómetros, a partir de esa distancia se produce un aumento en la semivarianza, que se mantiene hasta los 150 kilómetros, produciéndose una diferencia significativa entre los valores de las estaciones a esa distancia.

Una vez obtenido el variograma experimental es tiempo de ajustarlo utilizando el modelo que mejor se ajuste a la realidad de la correlación observada. En esta oportunidad se cambió el modelo Gaussiano con el que se venía trabajando por la parametrización de Matern, M. Stein. Este método es una técnica avanzada en geoestadística que permite modelar procesos espaciales complejos de manera más flexible que los enfoques tradicionales, particularmente cuando los datos espaciales muestran anisotropía o patrones de dependencia más complicados. De cierta manera se puede considerar a esta técnica como un método “Gaussiano mejorado”, pero es más correcto referirse a ella como una extensión que mejora el método Gaussiano en contextos donde las suposiciones estándares no son suficientes para capturar las características espaciales del fenómeno que estás modelando.

Considerando los componentes claves del variograma para el caso de junio se obtuvo un valor de nugget de 0. Para el caso del sill el valor quedo establecido en 704 y el range para este mes en específico mostró que pasado los 390 kilómetros la semivarianza deja de aumentar respecto a la distancia.

Análisis y descripción variográfica anisotrópica para las cuatro direcciones principales del mes de junio

# Crear el variograma experimental anisotrópico en diferentes direcciones
variograma_anisotropico_junio <- variogram(temp_promedio ~ 1, data_temp_filtered_junio, 
                                           alpha = c(0, 45, 90, 135),  # Direcciones a analizar
                                           cutoff = 160000, width = 6000, 
                                           tol.hor = 22.5)  # Tolerancia horizontal en grados

# Graficar el variograma experimental anisotrópico
plot(variograma_anisotropico_junio, main = "Variograma Anisotrópico - Mes: Junio", 
     xlab = "Distancia", ylab = "Semivarianza")  # Graficar variograma anisotrópico

# Definir modelo anisotrópico
model_vgm_aniso_junio <- vgm(psill = psill_ste_junio, model = "Ste",range = range_ste_junio, nugget = nugget_ste_junio, kappa = 1.3, anis = c(45, 0.5))  # Ajuste de anisotropía

# Ajustar el modelo anisotrópico al variograma experimental
fit_v_aniso_junio <- fit.variogram(variograma_anisotropico_junio, model_vgm_aniso_junio)
print(fit_v_aniso_junio)  # Mostrar el modelo ajustado
  model    psill    range kappa ang1 anis1
1   Nug   0.0000      0.0   0.0    0   1.0
2   Ste 709.3695 619779.4   1.3   45   0.5
# Graficar el variograma ajustado anisotrópico
plot(variograma_anisotropico_junio, fit_v_aniso_junio, 
     main = "Variograma Ajustado Anisotrópico - Mes: Junio", 
     xlab = "Distancia", ylab = "Semivarianza")  # Graficar

# Guardar el modelo ajustado anisotrópico
saveRDS(fit_v_aniso_junio, "data3/procesada/fit_variograma_aniso_junio.rds")  # Guardar en un archivo .rds

Al igual que en el caso de marzo la dirección norte-sur destaca por presentar valores de semivarianza bastante bajos, evidenciando una vez más que para esta dirección el cambio en los valores de temperatura a medida que aumenta la distancia es poco significativo. Manteniendose valores de temperatura similares a lo largo de toda la región a los 0° de azimut. A diferencia de lo que ocurría en marzo, para junio el resto de las direcciones tienen comportamientos más parecidos con valores de semivarianza que se estabilizan a una distancia cercana a los 100 kilometros, sobretodo las direcciones intermedias de los 45° y los 135°. Para junio el variograma modelado que se definió correspondió a la parametrización de Matern, M. Stein que entregó los mejores resultados, adaptandose de buena manera a los 90° y 135°. Este metodo entregó un valor de sill de 709 y un range de 620 kilómetros.

5.3 Análisis variográfico para el mes de septiembre

Análisis y descripción variográfica isotrópica para el mes de septiembre

# Filtrar los datos para el mes de septiembre
data_temp_filtered_septiembre <- data_temp[data_temp$month == 9, ]  # Filtrar datos para septiembre

# Crear el variograma experimental de las temperaturas
v_septiembre <- variogram(temp_promedio ~ 1, data_temp_filtered_septiembre, 
                          cutoff = 160000, width = 6000)  # Variograma experimental

# Graficar el variograma experimental
plot(v_septiembre, main = "Variograma Experimental - Mes: Septiembre", 
     xlab = "Distancia", ylab = "Semivarianza")  # Visualizar el variograma

# Ajustar el variograma usando autoajuste

auto_fit_septiembre <- autofitVariogram(temp_promedio ~ 1, data_temp_filtered_septiembre)  # Autoajuste del variograma

#| results: hide print(auto_fit_septiembre)  # Mostrar los resultados del ajuste

plot(auto_fit_septiembre)  # Graficar el modelo ajustado

# Definir manualmente los parámetros iniciales del modelo Ste (Stein)
nugget_ste_septiembre <- 0.43  # Estimación inicial del nugget
sill_ste_septiembre <- 225  # Estimación inicial del sill
psill_ste_septiembre <- sill_ste_septiembre - nugget_ste_septiembre  # Cálculo del psill
range_ste_septiembre <- 160000  # Estimación del rango

# Crear el modelo de variograma Stein (Ste)
model_vgm_ste_septiembre <- vgm(psill = psill_ste_septiembre, model = "Ste",range = range_ste_septiembre, nugget = nugget_ste_septiembre, kappa = 1.2) 

# Ajustar el modelo al variograma experimental
fit_v_ste_septiembre <- fit.variogram(v_septiembre, model_vgm_ste_septiembre)  
print(fit_v_ste_septiembre)  # Ver resultados del ajuste
  model    psill    range kappa
1   Nug    0.000      0.0   0.0
2   Ste 1155.812 500522.5   1.2
# Graficar el variograma ajustado
plot(v_septiembre, fit_v_ste_septiembre, main = "Variograma Ajustado - Mes: Septiembre", 
     xlab = "Distancia", ylab = "Semivarianza")  # Graficar el variograma ajustado

# Guardar el modelo ajustado en un archivo .rds
saveRDS(fit_v_ste_septiembre, "data3/procesada/fit_variograma_ste_septiembre.rds")  

El variograma experimental del mes de septiembre al mostrar resultados similares al del mes de junio se optó por modelarlo con la misma técnica de Matern, M. Stein que ajusta de mejor manera la función que el método Gaussiano tradicional. En este caso el valor del nugget fue de 0 (igual que en junio), el sill quedo establecido en 1.155 y por último, el range llego hasta los 500 kilómetros de distancia.

Análisis y descripción variográfica anisotrópica para las cuatro direcciones principales del mes de septiembre

# Crear el variograma experimental anisotrópico en diferentes direcciones
variograma_anisotropico_septiembre <- variogram(temp_promedio ~ 1, data_temp_filtered_septiembre, 
                                                alpha = c(0, 45, 90, 135),  # Direcciones a analizar
                                                cutoff = 160000, width = 6000, 
                                                tol.hor = 22.5)  # Tolerancia horizontal en grados

# Graficar el variograma experimental anisotrópico
plot(variograma_anisotropico_septiembre, main = "Variograma Anisotrópico - Mes: Septiembre", 
     xlab = "Distancia", ylab = "Semivarianza")  # Graficar variograma anisotrópico

# Definir modelo anisotrópico
model_vgm_aniso_septiembre <- vgm(psill = psill_ste_septiembre, model = "Ste",range = range_ste_septiembre, nugget = nugget_ste_septiembre, kappa = 1.2, anis = c(45, 0.5))  # Ajuste de anisotropía

# Ajustar el modelo anisotrópico al variograma experimental
fit_v_aniso_septiembre <- fit.variogram(variograma_anisotropico_septiembre, model_vgm_aniso_septiembre)
print(fit_v_aniso_septiembre)  # Mostrar el modelo ajustado
  model    psill    range kappa ang1 anis1
1   Nug    0.000      0.0   0.0    0   1.0
2   Ste 1168.481 796899.3   1.2   45   0.5
# Graficar el variograma ajustado anisotrópico
plot(variograma_anisotropico_septiembre, fit_v_aniso_septiembre, 
     main = "Variograma Ajustado Anisotrópico - Mes: Septiembre", 
     xlab = "Distancia", ylab = "Semivarianza")  # Graficar

# Guardar el modelo ajustado anisotrópico
saveRDS(fit_v_aniso_septiembre, "data3/procesada/fit_variograma_aniso_septiembre.rds")  # Guardar en un archivo .rds

Para el caso de los variogramas anisotrópicos del mes de septiembre se repite la situación para los 0°. Para los 90° los resultados de septiembre son bastante parecido al del mes anterior con valores de semivarianza muy importantes a los 150 kilometros. El caso de los 45° es interesante puesto que tiene una variabilidad casi nula hasta los 50 kilometros, luego experimenta un aumento y llegando a los 100 kilómetros denota una fuerte caida en sus valores de variabilidad. Para finalmente, alcanzar un aumentro progresivo hasta los 150 kilometros. El comportamiento de los 135° se observa muy similar a los dos meses anteriormente analisados. Para ajustar el variograma anisotropíco en este caso tambien se utilizó la parametrización de Matern, M. Stein. La que dio como resultado una función con un sill de 1.168 y un range de 797 kilometros.

5.4 Análisis variográfico para el mes de diciembre

Análisis y descripción variográfica isotrópica para el mes de diciembre

# Filtrar los datos para el mes de diciembre
data_temp_filtered_diciembre <- data_temp[data_temp$month == 12, ]  # Filtrar datos para diciembre

# Crear el variograma experimental de las temperaturas
v_diciembre <- variogram(temp_promedio ~ 1, data_temp_filtered_diciembre, 
                         cutoff = 160000, width = 6000)  # Variograma experimental

# Graficar el variograma experimental
plot(v_diciembre, main = "Variograma Experimental - Mes: Diciembre", 
     xlab = "Distancia", ylab = "Semivarianza")  # Visualizar el variograma

# Ajustar el variograma usando autoajuste

auto_fit_diciembre <- autofitVariogram(temp_promedio ~ 1, data_temp_filtered_diciembre)  # Autoajuste del variograma

#| results: hide print(auto_fit_diciembre)  # Mostrar los resultados del ajuste

plot(auto_fit_diciembre)  # Graficar el modelo ajustado

# Definir manualmente los parámetros iniciales del modelo Gaussiano para diciembre
nugget_gau_diciembre <- 0.58  # Estimación inicial del nugget
sill_gau_diciembre <- 149  # Estimación inicial del sill
psill_gau_diciembre <- sill_gau_diciembre - nugget_gau_diciembre  # Cálculo del psill
range_gau_diciembre <- 160000  # Estimación del rango

# Crear  el modelo de variograma Gaussiano
model_vgm_gau_diciembre <- vgm(psill = psill_gau_diciembre, model = "Gau", 
                               range = range_gau_diciembre, nugget = nugget_gau_diciembre)  # Modelo Gaussiano

# Ajustar el modelo al variograma experimental
fit_v_gau_diciembre <- fit.variogram(v_diciembre, model_vgm_gau_diciembre)  
print(fit_v_gau_diciembre)  # Ver resultados del ajuste
  model       psill    range
1   Nug   0.2103248      0.0
2   Gau 158.0183803 115931.1
# Graficar el variograma ajustado
plot(v_diciembre, fit_v_gau_diciembre, main = "Variograma Ajustado - Mes: Diciembre", 
     xlab = "Distancia", ylab = "Semivarianza")  # Graficar el variograma ajustado

# Guardar el modelo ajustado en un archivo .rds
saveRDS(fit_v_gau_diciembre, "data3/procesada/fit_variograma_gau_diciembre.rds")  

El mes de diciembre presenta un comportamiento variográfico similar al del mes de marzo, donde se aprecia como los valores a una corta distancia presentan una variabilidad en la temperatura relativamente baja. Esto se explica debido a que, a distancias cortas, los valores de temperatura deberían estar más correlacionados, generando una semivarianza baja. A partir de los 50 kilómetros de distancia entre las estaciones se aprecia como esta variable va aumentando de manera significativa disminuyendo la correlación entre los datos y aumentante la semivarianza.

Una vez obtenido el variograma experimental este se ajustó mediante una serie de modelos existentes para este fin. Como el variograma experimental es similar al de marzo se utilizó el mismo modelo de ajuste para ambos, el modelo Gaussiano.

Tomando en cuenta los componentes claves del variograma para el caso de diciembre se obtuvo un valor de nugget bastante bajo de 0,21, que es un muy buen resultado, puesto que la componente aleatoria del modelo es baja. Para el caso del sill se obtuvo que el valor donde ya no existe influencia entre las estaciones se produce cuando la semivarianza alcanza el valor de 158. La situación del rango para este mes en específico mostró que pasado los 116 kilómetros la semivarianza deja de aumentar respecto a la distancia de manera significativa.

Análisis y descripción variográfica anisotrópica para las cuatro direcciones principales del mes de diciembre

# Crear el variograma experimental anisotrópico en diferentes direcciones
variograma_anisotropico_diciembre <- variogram(temp_promedio ~ 1, data_temp_filtered_diciembre, 
                                               alpha = c(0, 45, 90, 135),  # Direcciones a analizar
                                               cutoff = 160000, width = 6000, 
                                               tol.hor = 22.5)  # Tolerancia horizontal en grados

# Graficar el variograma experimental anisotrópico
plot(variograma_anisotropico_diciembre, main = "Variograma Anisotrópico - Mes: Diciembre", 
     xlab = "Distancia", ylab = "Semivarianza")  # Graficar variograma anisotrópico

# Definir modelo anisotrópico
model_vgm_aniso_diciembre <- vgm(psill = psill_gau_diciembre, model = "Gau", 
                                 range = range_gau_diciembre, nugget = nugget_gau_diciembre, 
                                 anis = c(45, 0.5))  # Ajuste de anisotropía

# Ajustar el modelo anisotrópico al variograma experimental
fit_v_aniso_diciembre <- fit.variogram(variograma_anisotropico_diciembre, model_vgm_aniso_diciembre)
print(fit_v_aniso_diciembre)  # Mostrar el modelo ajustado
  model       psill    range ang1 anis1
1   Nug   0.5300975      0.0    0   1.0
2   Gau 168.8237290 192659.6   45   0.5
# Graficar el variograma ajustado anisotrópico
plot(variograma_anisotropico_diciembre, fit_v_aniso_diciembre, 
     main = "Variograma Ajustado Anisotrópico - Mes: Diciembre", 
     xlab = "Distancia", ylab = "Semivarianza")  # Graficar

# Guardar el modelo ajustado anisotrópico
saveRDS(fit_v_aniso_diciembre, "data3/procesada/fit_variograma_aniso_diciembre.rds")  # Guardar en un archivo .rds

El último mes del año se ajustó con el método de Gaussiano que entregó un buen ajuste para la dirección de los 135° de azimut, que considerando todos los meses trabajados es la dirección que presenta el valor mayor de semivarianza. En este sentido, dicha dirección que en la práctica se puede trazar con una línea desde la Provincia de Chacabuco hacia la parte sur de la comuna de San José de Maipo representa la dirección dentro de los modelos anisotrópicos generados la mayor variabilidad de la temperatura en la región. El modelo Gaussiano para diciembre entregó un sill de 169 y un range de 193 kilómetros.

El comportamiento anisotrópico de la temperatura en la región, por la escala de trabajo, puede ser explicado por la presencia de las tres macroformas principales del relieve que se extienden a lo largo del territorio en sentido norte-sur y que podrían estar actuando como factores homogeneizadores de los valores de temperatura. Cosa que no sucede cuando se analiza el comportamiento de la temperatura dentro del área de estudio en las direcciones donde estas formas no actúan de manera tan fuerte. Direcciones este-oeste (90°) y las direcciones intermedias.

06. Análisis variográfico de los residuos:

6.1 Residuos del mes de marzo

# Filtrar los datos para el mes de marzo
data_temp_filtered_marzo <- data_temp[data_temp$month == 3, ]  # Filtrar datos para MARZO

# Crear el variograma experimental para los residuos del modelo con predictores
vresiduos_marzo <- variogram(temp_promedio ~ MDE + Orientación + Pendiente + DistanciaCosta + NDVI_Diciembre + LST_Diciembre, data_temp_filtered_marzo, cutoff = 80000, width = 6000)

# Graficar el variograma experimental de los residuos     OK
plot(vresiduos_marzo)

# Ajustar el modelo de variograma para los residuos
auto_fit_residuos_marzo <- autofitVariogram(temp_promedio ~ MDE + Orientación + Pendiente + DistanciaCosta + 
NDVI_Diciembre + LST_Diciembre, data_temp_filtered_marzo)

#| results: hide print(auto_fit_residuos_marzo)

plot(auto_fit_residuos_marzo)

# Definir manualmente los parámetros iniciales del modelo Ste (Stein)
nugget_ste_marzo <- 2.18  # Estimación inicial del nugget
sill_ste_marzo <- 6  # Estimación inicial del sill
psill_ste_marzo <- sill_gau_marzo - nugget_gau_marzo  # Cálculo del psill
range_ste_marzo <- 70000  # Estimación del rango

# Crear el modelo de variograma Stein (Ste)
model_vgm_ste_marzo <- vgm(psill = psill_ste_marzo, model = "Ste", range = range_ste_marzo, nugget = nugget_ste_marzo, kappa = 0.3)

# Ajustar el modelo al variograma experimental
fit_v_ste_res_marzo <- fit.variogram(vresiduos_marzo, model_vgm_ste_marzo)  
print(fit_v_ste_res_marzo)  # Ver resultados del ajuste
  model    psill    range kappa
1   Nug 1.708251     0.00   0.0
2   Ste 2.882033 63792.72   0.3
# Graficar el variograma ajustado residuos
plot(vresiduos_marzo, fit_v_ste_res_marzo, main = "Variograma Ajustado Residuos- Mes: Marzo", 
     xlab = "Distancia", ylab = "Semivarianza")  # Graficar el variograma ajustado

# Guardar el modelo ajustado en un archivo .rds
saveRDS(fit_v_ste_res_marzo, "data3/procesada/fit_variograma_ste_res_marzo.rds")  

Para el caso del variograma experimental de los residuos de marzo se aprecia cierta tendencia creciente en los valores, lo que indica que los puntos cercanos tienen residuos más similares entre sí que los puntos distantes. Existiendo una parte de la correlación espacial no explicada. Para modelar el variograma experimental se utilizó el modelo de Matern, M. Stein alcanzando valores de nugget de 1.7 y de sill de 2.8.

6.2 Residuos del mes de junio

# Filtrar los datos para el mes de junio
data_temp_filtered_junio <- data_temp[data_temp$month == 6, ]  # Filtrar datos para JUNIO

# Crear el variograma experimental para los residuos del modelo con predictores
vresiduos_junio <- variogram(temp_promedio ~ MDE + Orientación + Pendiente + DistanciaCosta + NDVI_Junio + LST_Junio, data_temp_filtered_junio, cutoff = 80000, width = 6000)

# Graficar el variograma experimental de los residuos
plot(vresiduos_junio)

# Ajustar el modelo de variograma para los residuos
auto_fit_residuos_junio <- autofitVariogram(temp_promedio ~ MDE + Orientación + Pendiente + DistanciaCosta + NDVI_Junio + LST_Junio, data_temp_filtered_junio)

#| results: hide print(auto_fit_residuos_junio)

plot(auto_fit_residuos_junio)

# Definir manualmente los parámetros iniciales del modelo Ste (Stein)
nugget_ste_junio <- 1.25 # Estimación inicial del nugget
sill_ste_junio <- 4.5  # Estimación inicial del sill
psill_ste_junio <- sill_ste_junio - nugget_ste_junio  # Cálculo del psill
range_ste_junio <- 70000  # Estimación del rango

# Crear el modelo de variograma Stein (Ste)
model_vgm_ste_junio <- vgm(psill = psill_ste_junio, model = "Ste", range = range_ste_junio, nugget = nugget_ste_junio, kappa = 0.3)

# Ajustar el modelo al variograma experimental
fit_v_ste_res_junio <- fit.variogram(vresiduos_junio, model_vgm_ste_junio)  
print(fit_v_ste_res_junio)  # Ver resultados del ajuste
  model     psill    range kappa
1   Nug 0.6380626     0.00   0.0
2   Ste 3.2573781 43455.59   0.3
# Graficar el variograma ajustado residuos
plot(vresiduos_junio, fit_v_ste_res_junio, main = "Variograma Ajustado Residuos- Mes: Junio", 
     xlab = "Distancia", ylab = "Semivarianza")  # Graficar el variograma ajustado

# Guardar el modelo ajustado en un archivo .rds
saveRDS(fit_v_ste_res_junio, "data3/procesada/fit_variograma_ste_res_junio.rds")

Para el caso del variograma experimental de los residuos de junio se observa una tendencia creciente en los valores, lo que indica que los puntos cercanos tienen residuos más similares entre sí que los puntos distantes. Existiendo una parte de la correlación espacial no explicada. Para modelar el variograma experimental se utilizó el modelo de Matern, M. Stein alcanzando valores de nugget de 0.63 y de sill de 3.25.

6.3 Residuos del mes de septiembre

# Filtrar los datos para el mes de septiembre
data_temp_filtered_septiembre <- data_temp[data_temp$month == 9, ]  # Filtrar datos para SEPTIEMBRE

# Crear el variograma experimental para los residuos del modelo con predictores
vresiduos_septiembre <- variogram(temp_promedio ~ MDE + Orientación + Pendiente + DistanciaCosta + NDVI_Septiembre + LST_Septiembre, data_temp_filtered_septiembre, cutoff = 80000, width = 7000)

# Graficar el variograma experimental de los residuos
plot(vresiduos_septiembre)

# Ajustar el modelo de variograma para los residuos
auto_fit_residuos_septiembre <- autofitVariogram(temp_promedio ~ MDE + Orientación + Pendiente + DistanciaCosta + NDVI_Septiembre + LST_Septiembre, data_temp_filtered_septiembre)

#| results: hide  print(auto_fit_residuos_septiembre)

plot(auto_fit_residuos_septiembre)

# Definir manualmente los parámetros iniciales del modelo Ste (Stein)
nugget_ste_septiembre <- 0  # Estimación inicial del nugget
sill_ste_septiembre <- 2.7 # Estimación inicial del sill
psill_ste_septiembre <- sill_ste_septiembre - nugget_ste_septiembre  # Cálculo del psill
range_ste_septiembre <- 70000  # Estimación del rango

# Crear el modelo de variograma Stein (Ste)
model_vgm_ste_septiembre <- vgm(psill = psill_ste_septiembre, model = "Ste", range = range_ste_septiembre, nugget = nugget_ste_septiembre, kappa = 0.3)

# Ajustar el modelo al variograma experimental
fit_v_ste_res_septiembre <- fit.variogram(vresiduos_septiembre, model_vgm_ste_septiembre)  
print(fit_v_ste_res_septiembre)  # Ver resultados del ajuste
  model    psill    range kappa
1   Nug 1.409141      0.0   0.0
2   Ste 6.222621 585875.6   0.3
# Graficar el variograma ajustado residuos
plot(vresiduos_septiembre, fit_v_ste_res_septiembre, main = "Variograma Ajustado Residuos- Mes: Septiembre", 
     xlab = "Distancia", ylab = "Semivarianza")  # Graficar el variograma ajustado

# Guardar el modelo ajustado en un archivo .rds
saveRDS(fit_v_ste_res_septiembre, "data3/procesada/fit_variograma_ste_res_septiembre.rds")

Para el caso del variograma experimental de los residuos de septiembre se observa una tendencia creciente en los valores, lo que indica que los puntos cercanos tienen residuos más similares entre sí que los puntos distantes. Existiendo una parte de la correlación espacial no explicada. Para modelar el variograma experimental se utilizó el modelo de Matern, M. Stein alcanzando valores de nugget de 1.40 y de sill de 6.22.

6.4 Residuos del mes de diciembre

# Filtrar los datos para el mes de diciembre
data_temp_filtered_diciembre <- data_temp[data_temp$month == 12, ]  # Filtrar datos para DICIEMBRE

# Crear el variograma experimental para los residuos del modelo con predictores
vresiduos_diciembre <- variogram(temp_promedio ~ MDE + Orientación + Pendiente + DistanciaCosta + NDVI_Diciembre + LST_Diciembre, data_temp_filtered_diciembre, cutoff = 80000, width = 6000)

# Graficar el variograma experimental de los residuos
plot(vresiduos_diciembre)

# Ajustar el modelo de variograma para los residuos
auto_fit_residuos_diciembre <- autofitVariogram(temp_promedio ~ MDE + Orientación + Pendiente + DistanciaCosta + NDVI_Diciembre + LST_Diciembre, data_temp_filtered_diciembre)

#| results: hide print(auto_fit_residuos_diciembre)

plot(auto_fit_residuos_diciembre)

# Definir manualmente los parámetros iniciales del modelo Ste (Stein)
nugget_ste_diciembre <- 2.09  # Estimación inicial del nugget
sill_ste_diciembre <- 6  # Estimación inicial del sill
psill_ste_diciembre <- sill_ste_diciembre - nugget_ste_diciembre  # Cálculo del psill
range_ste_diciembre <- 70000  # Estimación del rango

# Crear el modelo de variograma Stein (Ste)
model_vgm_ste_diciembre <- vgm(psill = psill_ste_diciembre, model = "Ste", range = range_ste_diciembre, nugget = nugget_ste_diciembre, kappa = 0.3)

# Ajustar el modelo al variograma experimental
fit_v_ste_res_diciembre <- fit.variogram(vresiduos_diciembre, model_vgm_ste_diciembre)  
print(fit_v_ste_res_diciembre)  # Ver resultados del ajuste
  model    psill    range kappa
1   Nug 1.499005      0.0   0.0
2   Ste 5.809386 138736.8   0.3
# Graficar el variograma ajustado residuos
plot(vresiduos_diciembre, fit_v_ste_res_diciembre, main = "Variograma Ajustado Residuos- Mes: Diciembre", 
     xlab = "Distancia", ylab = "Semivarianza")  # Graficar el variograma ajustado

# Guardar el modelo ajustado en un archivo .rds
saveRDS(fit_v_ste_res_diciembre, "data3/procesada/fit_variograma_ste_res_diciembre.rds")

Para el caso del variograma experimental de los residuos de diciembre se observa una tendencia bastante creciente en los valores, lo que indica que los puntos cercanos tienen residuos más similares entre sí que los puntos distantes. Existiendo una parte de la correlación espacial no explicada. Para modelar el variograma experimental se utilizó el modelo de Matern, M. Stein alcanzando valores de nugget de 1.49 y de sill de 5.80.

l’esercizio è finito