Taller 4 Región Metropolitana

Magister en Teledetección. Análisis de Datos Espaciales con R. Interpolación espacial con Kriging

Introducción:

La interpolación espacial es una técnica fundamental en análisis geoespacial, que permite estimar valores de una variable en ubicaciones desconocidas a partir de puntos de datos conocidos. En el taller n°3, de análisis de autocorrelación espacial, se calcularon y modelaron una serie de variogramas necesarios para la continuidad del curso y alcanzar su objetivo final. En esta ocasión estos variogramas se utilizaron como insumos para realizar una serie de interpolaciones espaciales con distintos métodos de Kriging.

Los Kriging escogidos para ser trabajados correspondieron al Kriging Ordinario (OK) que es una técnica que asume los valores más cercanos tienden a ser más similares que los que están más alejados. Para capturar esta relación espacial, el Kriging utiliza el variograma. Por esta razón esta función es clave en el cálculo del Kriging, puesto que es esta la que asigna los pesos al conjunto de datos a diferencia de métodos no geoestadísticos que solo consideran la distancia. Para su obtención se utilizaron los variogramas modelados sin considerar ninguno de los predictores generados en el primer taller.

La segunda técnica utilizada fue la del Regresión-Kriging (RK) que es un método que junto con el Kriging mismo, en primer lugar, ajusta un modelo de regresión mediante predictores que aportan nueva información en la creación del producto final. Luego de terminado este proceso se calculan los residuos del análisis de esta regresión donde se considera la diferencia entre los valores observados y los predichos que son interpolados utilizando un Kriging Ordinario para considerar la variabilidad que la regresión no es capaz de determinar. Para este caso se utilizaron los variogramas modelados de los residuos, incorporando en el cálculo todos los predictores espaciales para que, finalmente, la función krige de R detecte automáticamente que el modelo que se busca obtener es un Regresión-Kriging.

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 que instala otros paquete 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 autentificarse 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 en esta ocasión se trató de cargar los datos de interés como los son el límite de la Región Metropolitana, los datos de temperatura con toda la información de los predictores y el ráster de los predictores mismos. Pero lo más importante fue incorporar al proceso los distintos variogramas para cada uno de los meses calculados y modelados en en taller anterior.

# Cargar el archivo geoespacial de la Región Metropolitana en formato GPKG
region <- st_read('data4/region_RM_UTM.gpkg')  

# Cargar los datos de temperatura con predictores
data_temp <- read_rds('data4/data_estas_temp_con_predictores.rds')

# Cargar los datos raster de predictores (MDE, orientación, pendiente, etc.)
preds <- rast('data4/predictores.tif')  

# Cargar el modelo de variograma ajustado para los meses de interes .rds
fit_v_gau_mar <- read_rds('data4/fit_variograma_gau_marzo.rds')  
fit_v_ste_jun <- read_rds('data4/fit_variograma_ste_junio.rds')
fit_v_ste_sep <- read_rds('data4/fit_variograma_ste_septiembre.rds') 
fit_v_gau_dic<- read_rds('data4/fit_variograma_gau_diciembre.rds') 

# Cargar el modelo de variograma de residuos ajustado para los meses de interes .rds
fit_v_ste_res_mar <- read_rds('data4/fit_variograma_ste_res_marzo.rds')
fit_v_ste_res_jun <- read_rds('data4/fit_variograma_ste_res_junio.rds')  
fit_v_ste_res_sep <- read_rds('data4/fit_variograma_ste_res_septiembre.rds')  
fit_v_ste_res_dic <- read_rds('data4/fit_variograma_ste_res_diciembre.rds')  

Debido al escaso conjunto de datos con los que se contó en la Región Metropolitana fue necesario incorporar una serie de “estaciones ficticias” con el objeto de complementar la información proporcionada por Agromet. Las estaciones ficticias fueron incorporaron, precisamente, en sectores con escasa información para mejorar los resultados de las prácticas en general. La ubicación de las estaciones de Agromet como de las nuevas estaciones incorporadas se presenta acontinuación:

# 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

03. Interpolación espacial mediante el método de Kriging Ordinario (OK):

En este punto, el objetivo fue generar un modelo predictivo de temperatura para los meses trabajados utilizando el método de Kriging Ordinario. Este método se basa únicamente en la autocorrelación espacial de los datos y requiere el ajuste de un modelo de variograma para describir la estructura espacial de la variable de interés. Los modelos de variogramas ajustados en el Taller 3 fueron utilizados para capturar adecuadamente la dependencia espacial de los datos de temperatura.

Para evaluar la precisión del modelo generado mediante Kriging Ordinario, se aplicó el método de validación cruzada leave-one-out (LOOCV). Este método permite medir el rendimiento del modelo comparando las predicciones en cada punto con el valor observado correspondiente, eliminando cada punto de manera iterativa y utilizando los demás para la predicción. El resultado de este análisis brindará métricas tales como el RMSE y R² que indicarán la calidad del ajuste.

tmap_mode('plot')

# Definir los meses de interés
meses <- c(3, 6, 9, 12)  # Marzo, Junio, Septiembre, Diciembre
nombre_meses <- c("Marzo", "Junio", "Septiembre", "Diciembre")

# Definir la función de RMSE
rmse <- function(obs, pred) sqrt(sum((obs - pred)^2) / length(obs))

# Definir las variables predictoras por mes
pred_vars_list <- list(
  c('MDE', 'Orientación', 'Pendiente', 'DistanciaCosta', 'NDVI_Marzo', 'LST_Marzo'),    # Marzo
  c('MDE', 'Orientación', 'Pendiente', 'DistanciaCosta', 'NDVI_Junio', 'LST_Junio'),    # Junio
  c('MDE', 'Orientación', 'Pendiente', 'DistanciaCosta', 'NDVI_Septiembre', 'LST_Septiembre'),  # Septiembre
  c('MDE', 'Orientación', 'Pendiente', 'DistanciaCosta', 'NDVI_Diciembre', 'LST_Diciembre')  # Diciembre
)

# Cargar los modelos de variograma por mes
fit_variogramas <- list(fit_v_gau_mar, fit_v_ste_jun, fit_v_ste_sep, fit_v_gau_dic)

# Crear un data.frame para almacenar todas las métricas de rendimiento
metricas_totales_OK <- data.frame(
  Mes = character(),
  RMSE = numeric(),
  R2 = numeric(),
  stringsAsFactors = FALSE
)

# Bucle para realizar interpolación por cada mes
for (i in 1:4) {
  cat(paste0("[using ordinary kriging] ", nombre_meses[i], "\n"))
  
  # Filtrar las capas predictoras por mes
  preds_mes <- subset(preds, pred_vars_list[[i]])
  preds_mes_st <- st_as_stars(preds_mes) |> split('band')
  
  # Filtrar los datos de temperatura por cada mes correspondiente
  data_temp_filtered <- data_temp |> filter(month == meses[i])
  
  # Obtener las coordenadas para las estaciones
  coords <- st_coordinates(data_temp_filtered)
  
  # Interpolación con OK (Kriging Ordinario)
  OK_model <- krige(temp_promedio ~1, data_temp_filtered, preds_mes_st, model = fit_variogramas[[i]])
  
  # Convertir la predicción a formato ráster
  OK_raster <- rast(OK_model)
  
  # Validación cruzada Leave-one-out para OK
  OK_loocv <- krige.cv(temp_promedio ~1, data_temp_filtered, model = fit_variogramas[[i]])
  
  # Calcular métricas de rendimiento (RMSE y R²)
  rmse_OK_loocv <- rmse(OK_loocv$observed, OK_loocv$var1.pred)
  r2_OK_loocv <- cor(OK_loocv$observed, OK_loocv$var1.pred)^2
  
  # Agregar las métricas al data.frame
  metricas_totales_OK <- rbind(
    metricas_totales_OK,
    data.frame(Mes = nombre_meses[i], RMSE = rmse_OK_loocv, R2 = r2_OK_loocv)
  )
  
  # Imprimir métricas de rendimiento
  cat(paste("RMSE para OK -", nombre_meses[i], ":", rmse_OK_loocv, "\n"))
  cat(paste("R² para OK -", nombre_meses[i], ":", r2_OK_loocv, "\n"))
  
  # Guardar el raster de predicciones
  writeRaster(OK_raster, filename = paste0("data4/procesada/OK_", tolower(nombre_meses[i]), "_predicciones_temp_promedio.tif"), 
              filetype = "GTiff", overwrite = TRUE)
  
  # Guardar resultados de validación cruzada en un RDS
  results_loocv_OK <- data.frame(
    Observed = OK_loocv$observed,
    Predicted = OK_loocv$var1.pred,
    Residuals = OK_loocv$observed - OK_loocv$var1.pred,
    ZScore = (OK_loocv$observed - OK_loocv$var1.pred) / sqrt(OK_loocv$var1.var),
    Mes = nombre_meses[i]  # Añadir mes para el análisis posterior
  )
  
  saveRDS(results_loocv_OK, paste0("data4/procesada/OK_", tolower(nombre_meses[i]), "_loocv_results.rds"))
}

[using ordinary kriging] Marzo
[using ordinary kriging]
RMSE para OK - Marzo : 1.50359864334868 
R² para OK - Marzo : 0.917209562803111 
[using ordinary kriging] Junio
[using ordinary kriging]
RMSE para OK - Junio : 1.78945493285607 
R² para OK - Junio : 0.948803723188103 
[using ordinary kriging] Septiembre
[using ordinary kriging]
RMSE para OK - Septiembre : 1.74259535615198 
R² para OK - Septiembre : 0.958797438959378 
[using ordinary kriging] Diciembre
[using ordinary kriging]
RMSE para OK - Diciembre : 1.39721752018912 
R² para OK - Diciembre : 0.959606546974774 

# Guardar todas las métricas de rendimiento en un archivo RDS
saveRDS(metricas_totales_OK, "data4/procesada/metricas_totales_OK.rds")

Para continuar con el taller se presentarán una serie de parámetros que tienen que ver con los resultados de las métricas del RMSE y R², la distribución de los residuos y del Z-Scores por mes y un gráfico con los valores observados vs los valores predichos, derivados del cálculo del Kriging Ordinario.

# Gráficos para comparar RMSE entre meses
ggplot(metricas_totales_OK, aes(x = Mes, y = RMSE, fill = Mes)) +
  geom_bar(stat = "identity", position = "dodge", na.rm = TRUE) +
  geom_text(aes(label = round(RMSE, 2)), position = position_dodge(width = 0.9), vjust = -0.3, na.rm = TRUE) +
  ggtitle("Comparación de RMSE entre Modelos por Mes") +
  xlab("MESES") +
  ylab("RMSE") +
  scale_fill_manual(name = "Mes", values = c("Marzo" = "steelblue", "Junio" = "orange", "Septiembre" = "green", "Diciembre" = "red")) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Gráficos para comparar R² entre meses
ggplot(metricas_totales_OK, aes(x = Mes, y = R2, fill = Mes)) +
  geom_bar(stat = "identity", position = "dodge", na.rm = TRUE) +
  geom_text(aes(label = round(R2, 2)), position = position_dodge(width = 0.9), vjust = -0.3, na.rm = TRUE) +
  ggtitle("Comparación de R² entre Modelos por Mes") +
  xlab("MESES") +
  ylab("R²") +
  scale_fill_manual(name = "Mes", values = c("Marzo" = "steelblue", "Junio" = "orange", "Septiembre" = "green", "Diciembre" = "red")) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Graficar manteniendo el orden de los modelos por mes
data_cv_OK_long <- metricas_totales_OK |> 
  pivot_longer(cols = c(RMSE, R2), names_to = "Métrica", values_to = "Valor")

# Kriging Ordinario (RMSE y R² entre meses)
(plot_OK <- ggplot(data_cv_OK_long, aes(x = Mes, y = Valor, fill = Métrica)) +
    geom_col(alpha = 0.7, na.rm = TRUE, position = "dodge") +  # Usar position = "dodge" para separar las barras
    ggtitle("Ordinary Kriging (RMSE y R² entre meses)") +
    xlab("MESES") +
    ylab("Valor de la Métrica") +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
    scale_fill_manual(name = "Métrica", values = c("RMSE" = "steelblue", "R2" = "orange"))  # Colores por métrica
)

# Convertir la lista a un dataframe para cada mes
loocv_data_list <- lapply(meses, function(mes) {
  readRDS(paste0("data4/procesada/OK_", tolower(nombre_meses[mes/3]), "_loocv_results.rds"))
})

# Unir los datos en un solo dataframe
loocv_data <- bind_rows(loocv_data_list, .id = "mes")

# Histograma de residuales
ggplot(loocv_data, aes(x = Residuals)) +
  geom_histogram(bins = 30, fill = "blue", alpha = 0.7) +
  facet_wrap(~ Mes) +
  labs(title = "Distribución de Residuos por Mes",
       x = "Residuales",
       y = "Frecuencia") +
  theme_minimal()

# Visualización de Z-Scores por mes
ggplot(loocv_data, aes(x = ZScore)) +
  geom_histogram(bins = 30, fill = "blue", alpha = 0.7) +
  facet_wrap(~ Mes) +
  labs(title = "Distribución de Z-Scores por Mes",
       x = "Z-Score",
       y = "Frecuencia") +
  theme_minimal()

# Gráfico de dispersión para valores observados vs predichos por mes
ggplot(loocv_data, aes(x = Observed, y = Predicted, color = Mes)) +
  geom_point(alpha = 0.7) +
  geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +  # Línea de referencia
  labs(title = "Predicciones vs Observaciones por Mes",
       x = "Valores Observados",
       y = "Valores Predichos") +
  theme_minimal() +
  facet_wrap(~ Mes)

Creación y representación de los modelos predictivos para los meses trabajados junto con sus valores de varianza

# Mostrar Mapas Ordinary Kriging de Predición y Variancia de temperatura en la RM.

# Crear una lista para almacenar los mapas
mapas <- list()

# Cargar los mapas desde los archivos guardados
for (i in 1:4) {
  mes_nombre <- nombre_meses[i]
  
  # Cargar el ráster de predicción
  pred_raster <- rast(paste0("data4/procesada/OK_", tolower(mes_nombre), "_predicciones_temp_promedio.tif"))
  
  # La varianza está en la segunda banda del mismo ráster
  var_raster <- pred_raster[[2]]  # Accede a la segunda banda del raster
  
  # Crear el mapa para la predicción de temperatura promedio
  mapa_pred <- tm_shape(pred_raster[[1]]) +
    tm_raster(style = 'quantile', 
              palette = viridis(20), 
              title = "Temperatura Promedio (°C)", 
              colorNA = "grey90", 
              alpha = 0.9) +
    tm_shape(region) +
    tm_borders(lwd = 3) +
    tm_shape(data_temp_filtered) +
    tm_dots(col = 'temp_promedio', 
            title = "Temperatura Media (°C)", 
            style = 'jenks', 
            palette = '-RdYlBu', 
            size = 0.1, 
            alpha = 1) +
    tm_compass(position = c("left", "top"), type = "4star", size = 0.9) +
    tm_grid() +
    tm_layout(
      main.title = paste("OK Predicción Temperatura -", mes_nombre),  
      main.title.size = 1,  
      main.title.position = "center",  
      main.title.fontface = "bold",  
      legend.outside = TRUE,  
      legend.title.size = 0.9,  
      legend.text.size = 0.7,  
      legend.position = c("left", "bottom"),  
      frame = FALSE  
    )
  
  # Crear el mapa para la varianza de la predicción
  mapa_var <- tm_shape(var_raster) +
    tm_raster(style = 'quantile', 
              palette = viridis(20), 
              title = "Varianza de Predicción", 
              colorNA = "grey90",  
              alpha = 0.9) +
    tm_shape(region) +
    tm_borders(lwd = 3) +
    tm_shape(data_temp_filtered) +
    tm_dots(col = 'temp_promedio', 
            title = "Temperatura Media (°C)", 
            style = 'jenks', 
            palette = '-RdYlBu', 
            size = 0.1, 
            alpha = 1) +
    tm_compass(position = c("left", "top"), type = "4star", size = 0.9) +
    tm_grid() +
    tm_layout(
      main.title = paste("OK Varianza de la Predicción -", mes_nombre),  
      main.title.size = 1,  
      main.title.position = "center",  
      main.title.fontface = "bold",  
      legend.outside = TRUE,  
      legend.title.size = 0.9,  
      legend.text.size = 0.7,  
      legend.position = c("left", "bottom"),  
      frame = FALSE  
    )
  
  # Almacenar mapas en la lista
  mapas[[paste("Predicción -", mes_nombre)]] <- mapa_pred
  mapas[[paste("Varianza -", mes_nombre)]] <- mapa_var
}

# Mostrar en dos etapas, para mejorar la visualizaciones por meses evaluados.

# Primera mitad
cat("Mostrando mapas para los meses: Marzo y Junio\n")

Mostrando mapas para los meses: Marzo y Junio

tmap_arrange(mapas[1:4], ncol = 2, nrow = 2)

# Segunda mitad
cat("Mostrando mapas para los meses: Septiembre y Diciembre\n")

Mostrando mapas para los meses: Septiembre y Diciembre

tmap_arrange(mapas[5:8], ncol = 2, nrow = 2)

04. Interpolación espacial mediante el método de Regresión-Kriging (RK):

El siguiente paso consistió en realizar un segundo medelo predictivo para cada mes utilizando el método de regresión-kriging, un método que combina dos técnicas la regresión lineal múltiple y el Kriging. Para este caso, se utilizó el modelo de regresión que obtuvo el mayor R² y el menor RMSE en el Taller n°2, lo que ayudó a incorporar predictores externos como lo fueron los de tipo topográfico, NDVI y la LST que mejoraron el resultado de las predicciones. Además, se emplearon los variogramas modelados que fueron realizados en el Taller n°3 para modelar los residuos espaciales. Finalmente, se generarán gráficos interactivos utilizando {tmap} para visualizar tanto los valores estimados como la varianza del error de predicción.

Al igual que con Kriging Ordinario, se aplicó el método de validación cruzada leave-one-out (LOOCV) para evaluar el desempeño de la interpolación basada en Regresión-Kriging. Esta validación permitió analizar la precisión del modelo combinando los predictores externos con la autocorrelación espacial de los residuos, proporcionando métricas clave como el RMSE y el R² para medir su rendimiento.

tmap_mode('plot')

# Definir los meses de interés
meses <- c(3, 6, 9, 12)  # Marzo, Junio, Septiembre, Diciembre
nombre_meses <- c("Marzo", "Junio", "Septiembre", "Diciembre")

# Definir la función de RMSE
rmse <- function(obs, pred) sqrt(sum((obs - pred)^2) / length(obs))

# Definir las variables predictoras por mes
pred_vars_list <- list(
  c('MDE', 'Orientación', 'Pendiente', 'DistanciaCosta', 'NDVI_Marzo', 'LST_Marzo'),    # Marzo
  c('MDE', 'Orientación', 'Pendiente', 'DistanciaCosta', 'NDVI_Junio', 'LST_Junio'),    # Junio
  c('MDE', 'Orientación', 'Pendiente', 'DistanciaCosta', 'NDVI_Septiembre', 'LST_Septiembre'),  # Septiembre
  c('MDE', 'Orientación', 'Pendiente', 'DistanciaCosta', 'NDVI_Diciembre', 'LST_Diciembre')  # Diciembre
)
# Cargar los modelos de variogramas ajustados de residuos por mes
fit_variogramas_res <- list(fit_v_ste_res_mar, fit_v_ste_res_jun, fit_v_ste_res_sep, fit_v_ste_res_dic)

# Crear un data.frame para almacenar todas las métricas de rendimiento
metricas_totales_RK <- data.frame(
  Mes = character(),
  RMSE = numeric(),
  R2 = numeric(),
  stringsAsFactors = FALSE
)

# Bucle para realizar interpolación por cada mes
for (i in 1:4) {
  cat(paste0("[using Regresión kriging] ", nombre_meses[i], "\n"))
  
  # Filtrar las capas predictoras por mes
  preds_mes <- subset(preds, pred_vars_list[[i]])
  preds_mes_st <- st_as_stars(preds_mes) |> split('band')
  
  # Filtrar los datos de temperatura por cada mes correspondiente
  data_temp_filtered <- data_temp |> filter(month == meses[i])
  
  # Obtener las coordenadas para las estaciones
  coords <- st_coordinates(data_temp_filtered)
  
  # Interpolación con RK (Regresión-kriging)
  RK_model <- krige(temp_promedio ~ MDE + Orientación + Pendiente + DistanciaCosta + 
                      get(paste0('NDVI_', nombre_meses[i])) + get(paste0('LST_', nombre_meses[i])), 
                    data_temp_filtered, preds_mes_st, model = fit_variogramas_res[[i]])
  
  # Convertir la predicción a formato ráster
  RK_raster <- rast(RK_model)
  
  # Validación cruzada Leave-one-out para RK
  RK_loocv <- krige.cv(temp_promedio ~ MDE + Orientación + Pendiente + DistanciaCosta + 
                         get(paste0('NDVI_', nombre_meses[i])) + get(paste0('LST_', nombre_meses[i])), 
                       data_temp_filtered, model = fit_variogramas_res[[i]])
  
  # Calcular métricas de rendimiento (RMSE y R²)
  rmse_RK_loocv <- rmse(RK_loocv$observed, RK_loocv$var1.pred)
  r2_RK_loocv <- cor(RK_loocv$observed, RK_loocv$var1.pred)^2
  
  # Agregar las métricas al data.frame
  metricas_totales_RK <- rbind(
    metricas_totales_RK,
    data.frame(Mes = nombre_meses[i], RMSE = rmse_RK_loocv, R2 = r2_RK_loocv)
  )
  
  # Imprimir métricas de rendimiento
  cat(paste("RMSE para RK -", nombre_meses[i], ":", rmse_RK_loocv, "\n"))
  cat(paste("R² para RK -", nombre_meses[i], ":", r2_RK_loocv, "\n"))
  
  # Guardar el raster de predicciones
  writeRaster(RK_raster, filename = paste0("data4/procesada/RK_", tolower(nombre_meses[i]), "_predicciones_temp_promedio.tif"), 
              filetype = "GTiff", overwrite = TRUE)
  
  # Guardar resultados de validación cruzada en un RDS
  results_loocv_RK <- data.frame(
    Observed = RK_loocv$observed,
    Predicted = RK_loocv$var1.pred,
    Residuals = RK_loocv$observed - RK_loocv$var1.pred,
    ZScore = (RK_loocv$observed - RK_loocv$var1.pred) / sqrt(RK_loocv$var1.var),
    Mes = nombre_meses[i]  # Añadir mes para el análisis posterior
  )
  
  saveRDS(results_loocv_RK, paste0("data4/procesada/RK_", tolower(nombre_meses[i]), "_loocv_results.rds"))
}

[using Regresión kriging] Marzo
[using universal kriging]
RMSE para RK - Marzo : 1.66765631306714 
R² para RK - Marzo : 0.898319651922219 
[using Regresión kriging] Junio
[using universal kriging]
RMSE para RK - Junio : 1.36188255039261 
R² para RK - Junio : 0.969309584226408 
[using Regresión kriging] Septiembre
[using universal kriging]
RMSE para RK - Septiembre : 1.57647206113482 
R² para RK - Septiembre : 0.965650793551374 
[using Regresión kriging] Diciembre
[using universal kriging]
RMSE para RK - Diciembre : 1.63277319791821 
R² para RK - Diciembre : 0.944868085810158 

# Guardar todas las métricas de rendimiento en un archivo RDS
saveRDS(metricas_totales_RK, "data4/procesada/metricas_totales_RK.rds")

A continuación se presentarán una serie de métricas y representaciones gráficas obtenidas mediante el método de Regresión-Kriging que tienen que ver con los resultados del RMSE y R², la distribución de los residuos y del Z-Scores por mes y los valores observados vs los valores predichos.

# Gráficos para comparar RMSE entre meses
ggplot(metricas_totales_RK, aes(x = Mes, y = RMSE, fill = Mes)) +
  geom_bar(stat = "identity", position = "dodge", na.rm = TRUE) +
  geom_text(aes(label = round(RMSE, 2)), position = position_dodge(width = 0.9), vjust = -0.3, na.rm = TRUE) +
  ggtitle("Comparación de RMSE entre Modelos por Mes") +
  xlab("MESES") +
  ylab("RMSE") +
  scale_fill_manual(name = "Mes", values = c("Marzo" = "steelblue", "Junio" = "orange", "Septiembre" = "green", "Diciembre" = "red")) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Gráficos para comparar R² entre meses
ggplot(metricas_totales_RK, aes(x = Mes, y = R2, fill = Mes)) +
  geom_bar(stat = "identity", position = "dodge", na.rm = TRUE) +
  geom_text(aes(label = round(R2, 2)), position = position_dodge(width = 0.9), vjust = -0.3, na.rm = TRUE) +
  ggtitle("Comparación de R² entre Modelos por Mes") +
  xlab("MESES") +
  ylab("R²") +
  scale_fill_manual(name = "Mes", values = c("Marzo" = "steelblue", "Junio" = "orange", "Septiembre" = "green", "Diciembre" = "red")) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Graficar manteniendo el orden de los modelos por mes
data_cv_RK_long <- metricas_totales_RK |> 
  pivot_longer(cols = c(RMSE, R2), names_to = "Métrica", values_to = "Valor")

# Regresión-Kriging (RMSE y R² entre meses)
(plot_RK <- ggplot(data_cv_RK_long, aes(x = Mes, y = Valor, fill = Métrica)) +
    geom_col(alpha = 0.7, na.rm = TRUE, position = "dodge") +  # Usar position = "dodge" para separar las barras
    ggtitle("Regresión kriging (RMSE y R² entre meses)") +
    xlab("MESES") +
    ylab("Valor de la Métrica") +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
    scale_fill_manual(name = "Métrica", values = c("RMSE" = "steelblue", "R2" = "orange"))  # Colores por métrica
)

# Convertir la lista a un dataframe para cada mes
loocv_data_list <- lapply(meses, function(mes) {
  readRDS(paste0("data4/procesada/RK_", tolower(nombre_meses[mes/3]), "_loocv_results.rds"))
})

# Unir los datos en un solo dataframe
loocv_data <- bind_rows(loocv_data_list, .id = "mes")

# Histograma de residuales
ggplot(loocv_data, aes(x = Residuals)) +
  geom_histogram(bins = 30, fill = "blue", alpha = 0.7) +
  facet_wrap(~ Mes) +
  labs(title = "Distribución de Residuos por Mes",
       x = "Residuales",
       y = "Frecuencia") +
  theme_minimal()

# Visualización de Z-Scores por mes
ggplot(loocv_data, aes(x = ZScore)) +
  geom_histogram(bins = 30, fill = "blue", alpha = 0.7) +
  facet_wrap(~ Mes) +
  labs(title = "Distribución de Z-Scores por Mes",
       x = "Z-Score",
       y = "Frecuencia") +
  theme_minimal()

# Gráfico de dispersión para valores observados vs predichos por mes
ggplot(loocv_data, aes(x = Observed, y = Predicted, color = Mes)) +
  geom_point(alpha = 0.7) +
  geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +  # Línea de referencia
  labs(title = "Predicciones vs Observaciones por Mes",
       x = "Valores Observados",
       y = "Valores Predichos") +
  theme_minimal() +
  facet_wrap(~ Mes)

Creación y representación de los modelos predictivos para los meses trabajados junto con sus valores de varianza

# Mostrar Mapas Regresión kriging de Predición y Variancia de temperatura en la RM.


# Crear una lista para almacenar los mapas
mapas <- list()

# Cargar los mapas desde los archivos guardados
for (i in 1:4) {
  mes_nombre <- nombre_meses[i]
  
  # Cargar el ráster de predicción
  pred_raster <- rast(paste0("data4/procesada/RK_", tolower(mes_nombre), "_predicciones_temp_promedio.tif"))
  
  # La varianza está en la segunda banda del mismo ráster
  var_raster <- pred_raster[[2]]  # Accede a la segunda banda del ráster
  
  # Crear el mapa para la predicción de temperatura promedio
  mapa_pred <- tm_shape(pred_raster[[1]]) +
    tm_raster(style = 'quantile', 
              palette = viridis(20), 
              title = "Temperatura Promedio (°C)", 
              colorNA = "grey90", 
              alpha = 0.9) +
    tm_shape(region) +
    tm_borders(lwd = 3) +
    tm_shape(data_temp_filtered) +
    tm_dots(col = 'temp_promedio', 
            title = "Temperatura Media (°C)", 
            style = 'jenks', 
            palette = '-RdYlBu', 
            size = 0.1, 
            alpha = 1) +
    tm_compass(position = c("left", "top"), type = "4star", size = 0.9) +
    tm_grid() +
    tm_layout(
      main.title = paste("RK Predicción Temperatura -", mes_nombre),  
      main.title.size = 1,  
      main.title.position = "center",  
      main.title.fontface = "bold",  
      legend.outside = TRUE,  
      legend.title.size = 0.9,  
      legend.text.size = 0.7,  
      legend.position = c("left", "bottom"),  
      frame = FALSE  
    )
  
  # Crear el mapa para la varianza de la predicción
  mapa_var <- tm_shape(var_raster) +
    tm_raster(style = 'quantile', 
              palette = viridis(20), 
              title = "Varianza de Predicción", 
              colorNA = "grey90",  
              alpha = 0.9) +
    tm_shape(region) +
    tm_borders(lwd = 3) +
    tm_shape(data_temp_filtered) +
    tm_dots(col = 'temp_promedio', 
            title = "Temperatura Media (°C)", 
            style = 'jenks', 
            palette = '-RdYlBu', 
            size = 0.1, 
            alpha = 1) +
    tm_compass(position = c("left", "top"), type = "4star", size = 0.9) +
    tm_grid() +
    tm_layout(
      main.title = paste("RK Varianza de la Predicción -", mes_nombre),  
      main.title.size = 1,  
      main.title.position = "center",  
      main.title.fontface = "bold",  
      legend.outside = TRUE,  
      legend.title.size = 0.9,  
      legend.text.size = 0.7,  
      legend.position = c("left", "bottom"),  
      frame = FALSE  
    )
  
  # Almacenar mapas en la lista
  mapas[[paste("Predicción -", mes_nombre)]] <- mapa_pred
  mapas[[paste("Varianza -", mes_nombre)]] <- mapa_var
}

# Mostrar en dos etapas, para mejorar la visualizaciones por meses evaluados.

# Primera mitad
cat("Mostrando mapas para los meses: Marzo y Junio\n")

Mostrando mapas para los meses: Marzo y Junio

tmap_arrange(mapas[1:4], ncol = 2, nrow = 2)

# Segunda mitad
cat("Mostrando mapas para los meses: Septiembre y Diciembre\n")

Mostrando mapas para los meses: Septiembre y Diciembre

tmap_arrange(mapas[5:8], ncol = 2, nrow = 2)

05. Comparación de los resultados de validación cruzada entre los métodos de Kriging Ordinario y Regresión-Kriging:

En este último paso, se realizó una comparación detallada entre los resultados obtenidos de las validaciones cruzadas de Kriging Ordinario y Regresión-Kriging. Se analizaron las diferencias en términos de precisión (RMSE y R²) para determinar qué método se adaptó mejor a la predicción de la temperatura en los distintos meses estudiados. Esta comparación ofrecerá información valiosa sobre las ventajas y limitaciones de cada enfoque, así como sobre la influencia de los predictores externos en las estimaciones finales.

# Cargar las métricas guardadas en archivos RDS
metricas_OK <- read_rds("data4/procesada/metricas_totales_OK.rds")
metricas_RK <- read_rds("data4/procesada/metricas_totales_RK.rds")

# Añadir una columna de método para distinguir los resultados
metricas_OK <- metricas_OK  |> 
  mutate(Metodo = "OK")

metricas_RK <- metricas_RK |> 
  mutate(Metodo = "RK")

# Unir ambas tablas de métricas
metricas_comparadas <- bind_rows(metricas_OK, metricas_RK)

# Ver la tabla comparativa
print(metricas_comparadas)

         Mes     RMSE        R2 Metodo
1      Marzo 1.503599 0.9172096     OK
2      Junio 1.789455 0.9488037     OK
3 Septiembre 1.742595 0.9587974     OK
4  Diciembre 1.397218 0.9596065     OK
5      Marzo 1.667656 0.8983197     RK
6      Junio 1.361883 0.9693096     RK
7 Septiembre 1.576472 0.9656508     RK
8  Diciembre 1.632773 0.9448681     RK

# Crear gráficos comparativos de RMSE y R2 por mes y método

# Gráfico RMSE
rmse_plot <- ggplot(metricas_comparadas, aes(x = Mes, y = RMSE, fill = Metodo)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(title = "Comparación de RMSE entre OK y RK",
       x = "Mes", y = "RMSE") +
  theme_minimal()

# Gráfico R²
r2_plot <- ggplot(metricas_comparadas, aes(x = Mes, y = R2, fill = Metodo)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(title = "Comparación de R² entre OK y RK",
       x = "Mes", y = "R²") +
  theme_minimal()

# Mostrar los gráficos
print(rmse_plot)

print(r2_plot)

Ambos métodos, Kriging Ordinario (OK) y Regresión-Kriging (RK), fueron evaluados utilizando la validación cruzada leave-one-out (LOOCV), generando las métricas de RMSE (error cuadrático medio) y R² (coeficiente de determinación) para los meses de marzo, junio, septiembre y diciembre. Al observar sus resultados es posible apreciar que para marzo el modelo construido a través del método de OK tiene un mejor desempeño, ya que presenta un menor RMSE (1.50 vs 1.67) y un mayor R² (0.92 vs 0.90). Para el caso de junio el panorama es diferente y RK impone su modelo con un RMSE menor (1.36 vs 1.79) y un R² mayor (0.97 vs 0.95). En el mes de septiembre el modelo de RK vuelve a tener cifras superiores a su par de OK con un menor RMSE (1.58 vs 1.74) y un mayor R² (0.97 vs 0.96). Finalmente, para diciembre OK vuelve a mostrar un mejor rendimiento como lo hizo en el mes de marzo con valores de RMSE de 1.40 vs 1.63 y R² de 0.96 vs 0.94.

Discusión:

En general el enfoque RK puede ofrecer estimaciones más precisas al incluir predictores adicionales como lo fueron la pendiente, orientación, NDVI, entre otros. Esto permitió que se capturara mejor la variabilidad espacial de la temperatura en áreas donde estos predictores tienen una relación fuerte con el fenómeno de interés, como en junio y septiembre, donde RK muestra mejor ajuste (menor RMSE y mayor R²).

En términos más concretos el Kriging Ordinario como se basa únicamente en la autocorrelación espacial de los datos, sin considerar predictores externos. En situaciones donde no se disponen de buenos predictores, o donde la relación entre estos y la variable de interés no es fuerte, el OK es una opción robusta. En los meses de marzo y diciembre, OK obtiene mejores resultados, lo que sugiere que, en esos periodos, la autocorrelación espacial por sí sola es suficiente para describir el comportamiento de la temperatura.

Como se mencionó anteriormente el Regresión-Kriging combina un modelo de regresión sobre predictores auxiliares con Kriging para modelar los residuos espaciales. RK puede producir estimaciones más precisas cuando los predictores tienen una relación clara con la variable de interés. En junio y septiembre, RK mostró un mejor desempeño, sugiriendo que, en esos meses, los predictores utilizados capturaron bien la variabilidad espacial de la temperatura, mejorando los resultados de los modelos en comparación con los de OK.

Si los resultados de la validación cruzada indicaron que RK tiene un menor RMSE y mayor R² en los meses de junio y septiembre, se puede concluir que este método es más efectivo en esos periodos al combinar los predictores adicionales y la autocorrelación espacial. Por otro lado, en marzo y diciembre, el OK muestra un mejor ajuste, lo que podría ser debido a una menor dependencia de los predictores en dichos meses.

En este punto no solo es importante referirse a los resultados de los métodos de interpolación en si mismos, sino que también a los variogramas que están detrás de ellos y permiten generar los modelos con una alta precisión.

En el caso del Kriging Ordinario los variogramas ajustados jugaron un papel crucial en la modelización de la autocorrelación espacial. Para los meses de marzo y diciembre, se utilizó un modelo Gaussiano (“Gau”), que resultó ser el que mejor representó el comportamiento de la temperatura en esos periodos, caracterizados por una suavización gradual de la variabilidad espacial.

Sin embargo, en los meses de junio y septiembre, el modelo Gaussiano no logró capturar adecuadamente la variabilidad de la temperatura. Por lo tanto, se optó cambiar este modelo por la parametrización de Matern, M. Stein (“Ste”), la cual presentó una mayor flexibilidad permitiendo un mejor ajuste a las características espaciales de la temperatura en estos meses. Esto es especialmente relevante porque el modelo de Stein pudo manejar mejor los cambios abruptos en la variabilidad espacial, que aparentemente caracterizaron a estos meses.

Para el método de Regresión-Kriging, que incorpora predictores externos y modela los residuos espaciales, se optó por el mismo modelo de Matern, M. Stein para todos los meses. Esto aseguró que los residuos, tras ser ajustados por los predictores, estuvieran modelados con un enfoque más flexible que pueda capturar mejor la complejidad espacial remanente.

Entonces, es muy importante tener en cuenta que al comparar los resultados entre OK y RK, no solo estamos comparando sus resultados finales, sino también el ajuste de los modelos de variogramas utilizados en cada caso. La calidad de los modelos predictivos observados en junio y septiembre utilizando el método RK puede ser atribuida en buena parte al uso del modelo de Stein para los residuos, que permitió capturar mejor la variabilidad espacial una vez fue ajustada por los predictores.

De esta manera, las diferencias en los modelos de variogramas aplicados fueron un factor clave al interpretar los resultados y las métricas de validación cruzada. Los buenos resultados de OK en marzo y diciembre reflejaron un ajuste más adecuado con el modelo Gaussiano, mientras que RK mostró un mejor rendimiento en los meses de junio y septiembre debido tanto al uso de predictores auxiliares como a la modelización de residuos mediante el modelo de Stein.

l’esercizio è finito