Taller 2 Región Metropolitana

Magister en Teledetección. Análisis de Datos Espaciales con R. Métodos de Predicción Espacial Determinísticos y de Regresión

Introducción:

El objetivo de este segundo taller es generar una serie de modelos espaciales predictivos para la Región Metropolitana de Santiago considerando la variable temperatura para los meses de marzo, junio, septiembre y diciembre del año 2023. A partir de información descargada, procesada y presentada en el taller anterior. Estos modelos fueron construidos a través de métodos de interpolación clásicos considerados no geoestadísticos, puesto que para realizar sus cálculos utilizan algoritmos matemáticos fijos que no toman en cuenta la variabilidad espacial del conjunto de datos analizados.

Entre los métodos de interpolación utilizados se encuentran aquellos conocidos como determinísticos los cuales no consideran la estructura de autocorrelación espacial de los datos y no asumen que exista una correlación entre ellos, basando sus estimaciones únicamente en distancias o superficies de ajuste. Para el caso concreto del taller se trabajó con los métodos del vecino más próximo e inverso a la distancia (IDW). También se utilizaron métodos de regresión lineal aprovechando todos los predictores construidos con anterioridad. En este caso este tipo de métodos se basan en la definición de una variable dependiente y una o más independientes.

Para cada uno de los métodos mencionados se realizó un proceso de validación cruzada denominado “Leave one out” (LOOCV) para obtener las métricas del error medio cuadrático (RMSE) y del coeficiente de determinación (R²).

00. Instalación de los paquetes requeridos:

#install.packages(c('tidyverse','sf','terra','tmap','remotes','geodata','rstac','earthdatalogin','viridis','gt','quarto','tinytex', 'knitr', 'rmarkdown','gridExtra'))
#remotes::install_github("ODES-CHILE/agrometR")

01. Carga de las librerías requeridas:

library(remotes)          # Paquete que trae función para instalar otro paquete desde github
library(tidyverse)        # Colección de paquetes para análisis de datos
library(sf)               # Manejo de datos vectoriales
library(terra)            # Manejo de datos raster y vectoriales
library(tmap)             # Creación de mapas temáticos
library(geodata)          # Acceso a descarga de diferentes datos raster y vectoriales
library(rstac)            # Permite acceder a información de spatial temporal assets catalogs
library(agrometR)         # Permite acceder a datos de estaciones red Agromet
library(dplyr)            # Librería de manipulación de datos
library(earthdatalogin)   # Sirve para acceder a productos NASA 'EarthData' 
library(viridis)          # Proporciona paletas de colores perceptiblemente uniformes
library(gt)               # Paquete que permite la creación de tablas
library(quarto)           # Permite la creación de documentos en varios formatos
library(tinytex)          # Proporciona una versión ligera de LaTeX para compilar a PDF
library(knitr)            # Para la creación de documentos dinámicos
library(rmarkdown)        # Para generar documentos reproducibles en diferentes formatos
library(gstat)            # Se utiliza principalmente para análisis geoestadístico
library(stars)            # Esta diseñada para el manejo y análisis de datos raster y vectoriales
library(gridExtra)        #utiliza para facilitar la disposición y organización de gráficos

02. Preparación de los datos para las interpolaciones:

Para realizar una interpolación, bajo cualquier método, es necesario contar con un conjunto de datos de referencia (datos de temperatura en este caso) y un área donde se desee conocer la información de dichos datos de manera continua para toda su extensión. En el caso de este punto se cargó la información generada en el taller anterior, y lo más importante se construyó una grilla que sirvió de sustento para todos los resultados que serán obtenidos.

# Carga de información anteriormente generada

region <- st_read('data2/region_RM_UTM.gpkg')

data_temp <- read_rds('data2/data_estas_temp_con_predictores.rds')

preds <- rast('data2/predictores.tif')

# Creación de la grilla para las interpolaciones

bb <- st_bbox(region)
grilla <- rast(xmin = bb$xmin, xmax = bb$xmax, ymin = bb$ymin, ymax = bb$ymax, 
               res = 500, crs = "EPSG:32719") 

grilla_st <- grilla |> st_as_stars()  # Convertir a clase stars

Representación cartográfica de las estaciones Agromet y las nuevas estaciones incorporadas como complemento para el ejercicio práctico de los talleres

El objetivo último del curso se centró en la necesidad de analizar datos de temperatura de distintas estaciones meteorológicas mediante diferentes métodos. Debido al escaso conjunto de datos con los que se cuentó 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.

# Agregar una columna 'estacion_ficticia' que indica si la estación es ficticia o real
data_temp <- data_temp |> 
  mutate(estacion_ficticia = ifelse(station_id >= 1000, "Ficticia", "Real"))

# Agregar las estaciones de temperatura, diferenciando entre reales y ficticias
tmap_mode('plot')
tm_shape(region) + 
  tm_borders() +
  tm_shape(data_temp) +
  tm_dots(col = "estacion_ficticia", size = 0.2, 
          palette = c("blue", "red"),  # Color azul para reales, rojo para ficticias
          title = "Tipo de Estación") +
  tm_layout(legend.outside = TRUE,       # Colocar la leyenda afuera
            legend.title.size = 3,     # Tamaño del título de la leyenda
            legend.text.size = 1)      # Tamaño del texto de la leyenda

03. Métodos de predicción espacial determinísticos:

Los métodos de predicción espacial determinísticos son aquellos que no incluyen una componente de incertidumbre en sus cálculos, basandose en relaciones matematicas fijas. Para el caso del taller se trabajó con el método del vecino más próximo y el método inverso a la distancia.

3.1 Cálculo método del vecino más próximo:

El método del vecino más próximo es un metodo clásico de interpolación de caracter local.En esta técnica el valor de una zona desconocida se estima considerando su relación con el valor del punto más cercano que tenga un valor conocido. No se realiza ningún cálculo o promedio adicional, simplemente el resultado se ajusta en base al valor del vecino más cercano. Para el caso del taller en cuestión se realizaron 3 tipos de iteraciones para todos los meses trabajados. Una considerando los 3 punto más cercano (k=3), otra los 5 puntos más próximos (k=5) y finalmente, una última tomando como referencia los 10 puntos más cercanos (k=10). En definitiva, este método de estimación atribuye toda la ponderación en su cálculo al dato más cercano del área donde se desea interpolar, considerando en dicho cálculo las n observaciones disponibles.

# Definir la función para el proceso de interpolación y validación cruzada

interpolar_por_mes <- function(mes) {
  
  # Filtrar los datos para el mes específico
  
  data_temp_filtered <- data_temp |> filter(month == mes)
  
  # Obtener las coordenadas para las estaciones
  
  coords <- st_coordinates(data_temp_filtered)
  
  # Definir los valores de k para los vecinos
  
  kas <- c(3, 5, 10)
  
  # Aplicar validación cruzada y obtener las predicciones
  
  out_kas <- purrr::map_df(1:nrow(data_temp_filtered), \(i) { 
    purrr::map_df(kas, \(k) {  # Iterar en los valores de k
      # Interpolación con IDW excluyendo el i-ésimo punto usando la misma grilla
      temp_nn <- idw(temp_promedio ~ 1, data_temp_filtered[-i,], grilla_st, nmax = k)
      # Aplicar máscara para la región
      temp_nn <- rast(temp_nn) |> mask(vect(region)) |> trim()  
      # Extraer el valor interpolado en la posición i
      df <- extract(temp_nn, data_temp_filtered[i,])
      
      # Guardar el ráster para cada k y mes
      filename <- paste0("RasterVMP_mes_", mes, "_k", k, "_idw_", i, ".tif")
      writeRaster(temp_nn, file.path('data2/procesada/VMP', filename), overwrite = TRUE) 
      
      df |> mutate(k = k)
    })
  })
  
  # Repetir los valores observados para todos los k
  
  observados <- rep(data_temp_filtered$temp_promedio, each = length(kas))
  
  # Combinar datos de predicciones con los observados y las coordenadas
  
  data_cv_kas <- cbind(coords, out_kas, observados)
  
  # Definir la función de RMSE
  rmse <- function(obs, pred) sqrt(sum((obs - pred)^2) / length(obs))
  
  # Calcular los residuos y el Z-Score
  
  data_cv_kas <- data_cv_kas |> 
    mutate(residuos = observados - var1.pred_var1.pred,  # Ajustar el nombre de la columna
           zscore = scale(residuos))
  
  # Calcular métricas de evaluación
  
  data_cv_kas_sum <- data_cv_kas |> 
    group_by(k) |> 
    summarize(rmse = rmse(observados, var1.pred_var1.pred),  # Usar el nombre correcto de la columna
              r2 = cor(observados, var1.pred_var1.pred)^2) |> 
    mutate(mes = mes)  # Agregar el mes a las métricas
  
  return(list(data_cv_kas_sum = data_cv_kas_sum, data_cv_kas = data_cv_kas))
}

# Lista de meses para iterar: marzo (3), junio (6), septiembre (9), diciembre (12)

meses <- c(3, 6, 9, 12)

# Aplicar la función para cada mes y guardar los resultados

resultados_meses <- purrr::map(meses, \(mes) {
  cat("Procesando el mes:", mes, "\n")  # Mostrar en qué mes está trabajando
  interpolar_por_mes(mes)
})

# Extraer resultados

data_cv_kas_sum_lista <- purrr::map(resultados_meses, "data_cv_kas_sum")
data_cv_kas_lista <- purrr::map(resultados_meses, "data_cv_kas")

# Combinar todos los resultados en un solo data frame

data_cv_kas_sum_total <- bind_rows(data_cv_kas_sum_lista, .id = "mes")
data_cv_kas_total <- bind_rows(data_cv_kas_lista, .id = "mes")

# Asegúrate de que la columna 'mes' sea un factor con etiquetas correctas

data_cv_kas_sum_total$mes <- factor(data_cv_kas_sum_total$mes,
                                    levels = c(1, 2, 3, 4),
                                    labels = c("Marzo", "Junio", "Septiembre", "Diciembre"))

# Convertir mes a factor con etiquetas para data_cv_kas_total

data_cv_kas_total$mes <- factor(data_cv_kas_total$mes,
                                levels = c(1, 2, 3, 4),
                                labels = c("Marzo", "Junio", "Septiembre", "Diciembre"))

# Mostrar los resultados de métricas con nombres de meses

Resultados de las Métricas para el cálculo del vecino más próximo

Las métricas presentadas a continuación para cada mes correspondientes al RMSE que es una medida de la magnitud del error en las predicciones de un modelo y el R² que mide la proporción de la variabilidad que es explicada por el modelo son de gran utilidad para medir el rendimiento y presición de los productos construidos.

print(data_cv_kas_sum_total)

# A tibble: 12 × 4
       k  rmse    r2 mes       
   <dbl> <dbl> <dbl> <fct>     
 1     3  1.55 0.912 Marzo     
 2     5  1.54 0.914 Marzo     
 3    10  1.56 0.914 Marzo     
 4     3  1.44 0.966 Junio     
 5     5  1.46 0.965 Junio     
 6    10  1.56 0.962 Junio     
 7     3  1.54 0.968 Septiembre
 8     5  1.59 0.966 Septiembre
 9    10  1.68 0.962 Septiembre
10     3  1.50 0.954 Diciembre 
11     5  1.54 0.952 Diciembre 
12    10  1.63 0.947 Diciembre 

Comparación de error medio cuadrático entre modelos por mes

ggplot(data_cv_kas_sum_total, aes(x = as.factor(k), 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("k") +
  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))

Comparación del coeficiente de determinación entre modelos por mes

ggplot(data_cv_kas_sum_total, aes(x = as.factor(k), 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("k") +
  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))

Visualización de residuos por mes y k

ggplot(data_cv_kas_total, aes(x = residuos, fill = as.factor(k))) +
  geom_histogram(bins = 30, alpha = 0.7, position = "identity") +
  labs(title = "Distribución de Residuos por Mes y k",
       x = "Residuos",
       y = "Frecuencia") +
  theme_minimal() +
  facet_grid(mes ~ k) +  # Facetas por mes y k
  scale_fill_viridis_d(name = "k")  # Colores diferenciados para k

Visualización de Z-Scores por mes y k

ggplot(data_cv_kas_total, aes(x = zscore, fill = as.factor(k))) +
  geom_histogram(bins = 30, alpha = 0.7, position = "identity") +
  labs(title = "Distribución de Z-Scores por Mes y k",
       x = "Z-Score",
       y = "Frecuencia") +
  theme_minimal() +
  facet_grid(mes ~ k) +  # Facetas por mes y k
  scale_fill_viridis_d(name = "k")  # Colores diferenciados para k

Gráfico de dispersión para valores observados vs valores predichos por mes

ggplot(data_cv_kas_total, aes(x = observados, y = var1.pred_var1.pred, color = as.factor(k))) +
  geom_point(alpha = 0.7) +
  geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +  # Línea de referencia
  labs(title = "Predicciones vs Observados por Mes y k",
       x = "Valores Observados",
       y = "Valores Predichos") +
  theme_minimal() +
  facet_grid(mes ~ k) +  # Facetas por mes y k
  scale_color_viridis_d(name = "k")  # Colores diferenciados para k

Creación y promedios para cada uno de los meses trabajados y los valores k considerados

#Lista de meses y valores de k
meses <- c(3, 6, 9, 12)
kas <- c(3, 5, 10)

# Iterar sobre meses y valores de k
for (mes in meses) {
  for (k in kas) {
    # Cargar los rásteres de las diferentes iteraciones para un mes y un valor de k
    raster_files <- list.files('data2/procesada/VMP', 
                               pattern = paste0("RasterVMP_mes_", mes, "_k", k, "_idw_"), 
                               full.names = TRUE)
    
    # Cargar todos los rásteres en una lista, asegurándose de que son objetos raster válidos
    rasters_list <- lapply(raster_files, function(file) {
      raster <- tryCatch({
        rast(file)
      }, error = function(e) {
        message("Error al cargar el archivo: ", file)
        return(NULL) # Devolver NULL si hay un error
      })
      return(raster)
    })
    
    # Filtrar los objetos NULL (los que no se pudieron cargar)
    rasters_list <- Filter(Negate(is.null), rasters_list)
    
    # Verificar si la lista de rásteres no está vacía
    if (length(rasters_list) > 0) {
      # Convertir la lista a un objeto SpatRaster
      raster_stack <- do.call(c, rasters_list) # Combina las listas de rásteres
      
      # Calcular el promedio de los rásteres usando terra::app
      raster_promedio <- app(raster_stack, fun = mean, na.rm = TRUE)
      
      # Guardar el ráster promedio
      filename_promedio <- paste0("RasterVMP_mes_", mes, "_k", k, "_promedio.tif")
      writeRaster(raster_promedio, file.path('data2/procesada/VMP', filename_promedio), overwrite = TRUE)
      
      cat("Promedio guardado para mes:", mes, "y k:", k, "\n")
    } else {
      cat("No se encontraron rásteres válidos para mes:", mes, "y k:", k, "\n")
    }
  }
}

# Cargar los rasters desde los archivos por mes y k seleccionados
raster_marzo <- rast("data2/procesada/VMP/RasterVMP_mes_3_k3_promedio.tif")    # Raster de marzo
raster_junio <- rast("data2/procesada/VMP/RasterVMP_mes_6_k3_promedio.tif")    # Raster de junio
raster_septiembre <- rast("data2/procesada/VMP/RasterVMP_mes_9_k3_promedio.tif") # Raster de septiembre
raster_diciembre <- rast("data2/procesada/VMP/RasterVMP_mes_12_k3_promedio.tif") # Raster de diciembre

# Crear un stack de rasters
RasterVMP_stack <- c(raster_marzo, raster_junio, raster_septiembre, raster_diciembre)

# Nombres para los rasters en el stack
names(RasterVMP_stack) <- c("Marzo", "Junio", "Septiembre", "Diciembre")

Cartografías utilizando el método del vecino más próximo (k=3)

tm_shape(RasterVMP_stack) +  # Añadir el ráster de temperatura
  tm_raster(style = 'quantile', 
            palette = viridis(20), 
            title = "Rango de Temperatura (°C)", 
            colorNA = "grey90",  # Color para valores NA
            alpha = 0.9) +  # Ajustar la transparencia del ráster
  tm_shape(region) +  # Añadir la capa de la región
  tm_borders(lwd = 3) +  # Dibujar los bordes de la región con mayor grosor
  tm_shape(data_temp) +  # Añadir los puntos de datos de temperatura
  tm_dots(col = 'temp_promedio', 
          title = "Temperatura Media (°C)",  # Título descriptivo de la leyenda
          style = 'jenks', 
          palette = '-RdYlBu', 
          size = 0.3,  # Aumentar el tamaño de los puntos
          alpha = 1) +  # Ajustar la transparencia de los puntos
  tm_compass(position = c("left", "top"), type = "4star", size = 0.7) +  # Añadir flecha del norte
  tm_grid() +  # Añadir una cuadrícula
  tm_layout(
    main.title = "Vecino más Cercano (k = 3)",  # Título principal
    main.title.size = 1,  # Tamaño del título
    main.title.position = "center",  # Centrar el título
    main.title.fontface = "bold",  # Negrita para el título
    legend.outside = TRUE,  # Colocar la leyenda afuera
    legend.title.size = 0.9,  # Tamaño del título de la leyenda
    legend.text.size = 0.7,  # Tamaño del texto de la leyenda
    legend.position = c("left", "bottom"),  # Posición de la leyenda
    frame = FALSE  # Eliminar el marco del mapa
  )

tmap_mode('plot')

3.2 Cálculo método del inverso a la distancia (IDW):

El método de interpolación del inverso a la distancia es un segundo método clásico de caracter local. La principal diferencia con el método anterior es que es este caso el resultado para cada celda se asigna en función de una media ponderada de varios puntos conocidos. La influencia de cada punto depende de su distancia al punto de estimación, con los puntos más cercanos teniendo más peso. Una desventaja importante en este caso es que la técnica al trabajar com promedios no puede obtener valores mayores ni menores que los valores máximos y minimos de las observaciones, limitando sus resultados.

# Definir la función para el proceso de interpolación y validación cruzada

interpolar_por_mes <- function(mes) {
  
  # Filtrar los datos para el mes específico
  
  data_temp_filtered <- data_temp |> 
    filter(month == mes)
  
  # Obtener las coordenadas para las estaciones
  
  coords <- st_coordinates(data_temp_filtered)
  
  # Aplicar validación cruzada y obtener las predicciones
  
  out_IDW <- purrr::map_df(1:nrow(data_temp_filtered), \(i) { 
    # Interpolación con IDW excluyendo el i-ésimo punto
    temp_nn <- idw(temp_promedio ~ 1, data_temp_filtered[-i,], grilla_st)
    # Aplicar máscara para la región y ajustar el raster
    temp_nn <- rast(temp_nn) |> mask(vect(region)) |> trim()  
    # Extraer el valor interpolado en la posición i
    df <- extract(temp_nn, data_temp_filtered[i,])
    
    # Guardar el ráster para cada mes
    filename <- paste0("RasterIDW_mes_", mes,"_idw_", i, ".tif")
    writeRaster(temp_nn, file.path('data2/procesada/IDW', filename), overwrite = TRUE)
    
    df
    
  })
  
  # Repetir los valores observados para todos los k
  
  observados_IDW <- data_temp_filtered$temp_promedio
  
  # Combinar datos de predicciones con los observados
  
  data_cv_IDW <- cbind(out_IDW, observados_IDW)
  
  # Definir la función de RMSE
  
  rmse <- function(obs, pred) sqrt(sum((obs - pred)^2) / length(obs))
  
  # Calcular los residuos y el Z-Score
  
  data_cv_IDW <- data_cv_IDW |> 
    mutate(residuos = observados_IDW - var1.pred_var1.pred,
           zscore = scale(residuos))
  
  # Calcular métricas de evaluación
  
  data_cv_IDW_sum <- data_cv_IDW |> 
    summarize(rmse = rmse(observados_IDW, var1.pred_var1.pred),
              r2 = cor(observados_IDW, var1.pred_var1.pred)^2) |> 
    mutate(mes = mes)  # Agregar el mes a las métricas
  
  return(list(data_cv_IDW_sum = data_cv_IDW_sum, data_cv_IDW = data_cv_IDW))
}

# Lista de meses para iterar: marzo (3), junio (6), septiembre (9), diciembre (12)

meses <- c(3, 6, 9, 12)

# Aplicar la función para cada mes y guardar los resultados

resultados_meses <- purrr::map(meses, \(mes) {
  cat("Procesando el mes:", mes, "\n")  # Mostrar en qué mes está trabajando
  interpolar_por_mes(mes)
})

# Extraer resultados

data_cv_IDW_sum_lista <- purrr::map(resultados_meses, "data_cv_IDW_sum")
data_cv_IDW_lista <- purrr::map(resultados_meses, "data_cv_IDW")

# Combinar todos los resultados en un solo data frame

data_cv_IDW_sum_total <- bind_rows(data_cv_IDW_sum_lista, .id = "mes")
data_cv_IDW_total <- bind_rows(data_cv_IDW_lista, .id = "mes")

# Asegúrate de que la columna 'mes' sea un factor con etiquetas correctas

data_cv_IDW_sum_total$mes <- factor(data_cv_IDW_sum_total$mes,
                                    levels = c(1, 2, 3, 4),
                                    labels = c("Marzo", "Junio", "Septiembre", "Diciembre"))

# Convertir mes a factor con etiquetas para data_cv_IDW_total

data_cv_IDW_total$mes <- factor(data_cv_IDW_total$mes,
                                levels = c(1, 2, 3, 4),
                                labels = c("Marzo", "Junio", "Septiembre", "Diciembre"))

Resultados de las métricas para el cálculo del inverso a la distancia (IDW)

print(data_cv_IDW_sum_total)

      rmse        r2        mes
1 1.944778 0.8956168      Marzo
2 2.212957 0.9503709      Junio
3 2.393987 0.9541756 Septiembre
4 2.156680 0.9353780  Diciembre

Comparación de error medio cuadrático entre modelos por mes

ggplot(data_cv_IDW_sum_total, aes(x = 'IDW', 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))

Comparación del coeficiente de determinación entre modelos por mes

ggplot(data_cv_IDW_sum_total, aes(x = 'IDW', 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))

Visualización de los residuos por mes

ggplot(data_cv_IDW_total, aes(x = residuos, fill = mes)) +
  geom_histogram(bins = 30, alpha = 0.7) +
  labs(title = "Distribución de Residuos por Mes",
       x = "Residuos",
       y = "Frecuencia") +
  theme_minimal() +
  facet_wrap(~mes)

Visualización de Z-Scores por mes

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

Gráfico de dispersión para valores observados vs valores predichos por mes

ggplot(data_cv_IDW_total, aes(x = observados_IDW, y = var1.pred_var1.pred, 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 Observados por Mes",
       x = "Valores Observados",
       y = "Valores Predichos") +
  theme_minimal() +
  facet_wrap(~mes)

Iteraciones para cada mes para obtener los promedios de temperatura

# Lista de meses
meses <- c(3, 6, 9, 12)

# Iterar sobre meses
for (mes in meses) {
  # Cargar los rásteres de las diferentes iteraciones para un mes
  raster_files <- list.files('data2/procesada/IDW', 
                             pattern = paste0("RasterIDW_mes_", mes, "_idw_"), 
                             full.names = TRUE)
  
  # Cargar todos los rásteres en una lista, asegurándose de que son objetos raster válidos
  rasters_list <- lapply(raster_files, function(file) {
    raster <- tryCatch({
      rast(file)
    }, error = function(e) {
      message("Error al cargar el archivo: ", file)
      return(NULL) # Devolver NULL si hay un error
    })
    return(raster)
  })
  
  # Filtrar los objetos NULL (los que no se pudieron cargar)
  rasters_list <- Filter(Negate(is.null), rasters_list)
  
  # Verificar si la lista de rásteres no está vacía
  if (length(rasters_list) > 0) {
    # Convertir la lista a un objeto SpatRaster
    raster_stack <- do.call(c, rasters_list) # Combina las listas de rásteres
    
    # Calcular el promedio de los rásteres usando terra::app
    raster_promedio <- app(raster_stack, fun = mean, na.rm = TRUE)
    
    # Guardar el ráster promedio
    filename_promedio <- paste0("RasterIDW_mes_", mes, "_promedio.tif")
    writeRaster(raster_promedio, file.path('data2/procesada/IDW', filename_promedio), overwrite = TRUE)
    
    cat("Promedio guardado para mes:", mes, "\n")
  } else {
    cat("No se encontraron rásteres válidos para mes:", mes, "\n")
  }
}

Carga de los rásters promedio y creación del stack

# Cargar los rasters promedios
raster_marzo <- rast("data2/procesada/IDW/RasterIDW_mes_3_promedio.tif")
raster_junio <- rast("data2/procesada/IDW/RasterIDW_mes_6_promedio.tif")
raster_septiembre <- rast("data2/procesada/IDW/RasterIDW_mes_9_promedio.tif")
raster_diciembre <- rast("data2/procesada/IDW/RasterIDW_mes_12_promedio.tif")

# Crear un stack de rasters
RasterIDW_stack <- c(raster_marzo, raster_junio, raster_septiembre, raster_diciembre)

# Nombres para los rasters en el stack
names(RasterIDW_stack) <- c("Marzo", "Junio", "Septiembre", "Diciembre")

Cartografías utilizando el método del inverso a la distancia (IDW)

tm_shape(RasterIDW_stack) +  # Añadir el ráster de temperatura
  tm_raster(style = 'quantile', 
            palette = viridis(20), 
            title = "Rango de Temperatura (°C)", 
            colorNA = "grey90",  # Color para valores NA
            alpha = 0.9) +  # Ajustar la transparencia del ráster
  tm_shape(region) +  # Añadir la capa de la región
  tm_borders(lwd = 3) +  # Dibujar los bordes de la región con mayor grosor
  tm_shape(data_temp) +  # Añadir los puntos de datos de temperatura
  tm_dots(col = 'temp_promedio', 
          title = "Temperatura Media (°C)",  # Título descriptivo de la leyenda
          style = 'jenks', 
          palette = '-RdYlBu', 
          size = 0.3,  # Aumentar el tamaño de los puntos
          alpha = 1) +  # Ajustar la transparencia de los puntos
  tm_compass(position = c("left", "top"), type = "4star", size = 0.7) +  # Añadir flecha del norte
  tm_grid() +  # Añadir una cuadrícula
  tm_layout(
    main.title = "Inverso de la Distancia IDW",  # Título principal
    main.title.size = 1,  # Tamaño del título
    main.title.position = "center",  # Centrar el título
    main.title.fontface = "bold",  # Negrita para el título
    legend.outside = TRUE,  # Colocar la leyenda afuera
    legend.title.size = 0.9,  # Tamaño del título de la leyenda
    legend.text.size = 0.7,  # Tamaño del texto de la leyenda
    legend.position = c("left", "bottom"),  # Posición de la leyenda
    frame = FALSE  # Eliminar el marco del mapa
  )

04. Regresión líneal considerando todos los predictores:

La regresión líneal centra sus esfuerzos en ajustar una línea o un plano que minimice las diferencias entre los valores observados y los valores predichos. Para el caso del presente taller se utilizó la función “krige” que se usa para realizar interpolaciones mediante el método de kriging. Reemplazando la utilización de modelos de regresión simples para las estimaciones. La diferencia en este caso es que solo se utilizó la parte deterministica del método. Por lo tanto el enfoque que se le dió a la función “krige” se centró únicamente en la tendencia general de los datos. En otras palabras, se dejo de lado la componenete estocástica de la técnica, asumiendo que la tendencia general es sufucientes para entender y modelar el fenomeno en cuestión.

# Utilizando la función krige de {gstat}, utiliza modelos de regresión lineal generalizada en vez de modelos de regresión simple.

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

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

# Crear una lista para almacenar los resultados
resultados <- list()
metricas <- data.frame(Mes = character(), RMSE = numeric(), R2 = numeric(), stringsAsFactors = FALSE)

for (i in seq_along(meses)) {
  mes <- meses[i]
  mes_nombre <- nombre_meses[i]
  
  # Filtrar las capas predictoras por mes
  pred_vars <- c('MDE', 'Orientación', 'Pendiente', 'DistanciaCosta',
                 paste0('NDVI_', mes_nombre), paste0('LST_', mes_nombre))
  
  preds_mes <- subset(preds, pred_vars)
  preds_mes_st <- preds_mes |> st_as_stars() |> split('band')
  
  # Filtrar los datos de temperatura por mes
  data_temp_filtered <- data_temp |> filter(month == mes)
  
  # Regresión lineal (Krige)
  temp_lm6 <- krige(temp_promedio ~ MDE + Orientación + Pendiente + DistanciaCosta +
                      get(paste0('NDVI_', mes_nombre)) + get(paste0('LST_', mes_nombre)),
                    data_temp_filtered, preds_mes_st)
  
  # Convertir la predicción a formato raster
  temp_lm6_raster <- rast(temp_lm6)
  
  # Validación cruzada Leave-one-out
  temp_lm6_loocv <- krige.cv(temp_promedio ~ MDE + Orientación + Pendiente + DistanciaCosta +
                               get(paste0('NDVI_', mes_nombre)) + get(paste0('LST_', mes_nombre)),
                             data_temp_filtered)
  
  # Calcular métricas de rendimiento (RMSE y R^2)
  rmse_loocv <- rmse(temp_lm6_loocv$observed, temp_lm6_loocv$var1.pred)
  r2_loocv <- cor(temp_lm6_loocv$observed, temp_lm6_loocv$var1.pred)^2
  
  # Almacenar las métricas en el data frame
  metricas <- rbind(metricas, data.frame(Mes = mes_nombre, RMSE = rmse_loocv, R2 = r2_loocv))
  
  # Almacenar los resultados en la lista
  resultados[[mes_nombre]] <- list(
    prediccion_raster = temp_lm6_raster,
    loocv = list(modelo = temp_lm6_loocv, metricas = c(rmse_loocv, r2_loocv))
  )
}

Cálculo de Métricas del modelo de regresión considerando todos los predictores

# Mostrar los resultados de métricas con nombres de meses

print(metricas)

         Mes     RMSE        R2
1      Marzo 2.135780 0.8332633
2      Junio 1.987863 0.9345702
3 Septiembre 1.866087 0.9518754
4  Diciembre 2.372220 0.8833696

Gráfico para comparar RMSE entre meses

ggplot(metricas, aes(x = 'Krige', 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áfico para comparar R² entre meses

ggplot(metricas, aes(x = 'Krige', 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))

Visualización de histogramas de residuales por mes

# Convertir la lista a un dataframe para cada mes

loocv_data_list <- lapply(resultados, function(x) x$loocv$modelo)

# Unir los datos en un solo dataframe

loocv_data <- bind_rows(loocv_data_list, .id = "mes")

# Histograma de residuales

ggplot(loocv_data, aes(x = residual)) +
  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 histogramas 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 = var1.pred, 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)

Cartografías predictoras de temperatura con su varianza para los meses de marzo y junio

# Crear una lista para almacenar los mapas

mapas <- list()

# Generar mapas para cada mes y almacenarlos en la lista
for (mes_nombre in nombre_meses) {
  temp_raster <- resultados[[mes_nombre]]$prediccion_raster
  
  # Crear el mapa para la predicción de temperatura promedio
  mapa_pred <- tm_shape(temp_raster[[1]]) +  # Añadir el ráster de predicción
    tm_raster(style = 'quantile', 
              palette = viridis(20), 
              title = "Temperatura Promedio (°C)", 
              colorNA = "grey90",  # Color para valores NA
              alpha = 0.9) +  # Ajustar la transparencia del ráster
    tm_shape(region) +  # Añadir la capa de la región
    tm_borders(lwd = 3) +  # Dibujar los bordes de la región con mayor grosor
    tm_shape(data_temp_filtered) +  # Añadir los puntos de datos de temperatura
    tm_dots(col = 'temp_promedio', 
            title = "Temperatura Media (°C)", 
            style = 'jenks', 
            palette = '-RdYlBu', 
            size = 0.3, 
            alpha = 1) +  
    tm_compass(position = c("left", "top"), type = "4star", size = 0.7) +  
    tm_grid() +  
    tm_layout(
      main.title = paste("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(temp_raster[[2]]) +  # Añadir el ráster de varianza
    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.3, 
            alpha = 1) +  
    tm_compass(position = c("left", "top"), type = "4star", size = 0.7) +  
    tm_grid() +  
    tm_layout(
      main.title = paste("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
}

# Si quieres mostrar en dos etapas, puedes usar:

# 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)

Cartografías predictoras de temperatura con su varianza para los meses de septiembre y diciembre

# 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)

Resultados comparados entre métodos del vecino más cercano, IDW y regresión líneal para todos los meses y modelos

06. Discusión:

Considerando el objetivo de este segundo taller donde se probaron distintos métodos de interpolación determinísticos y de regresión para distintos meses del año 2023 se lograron obtener una serie de métricas mediante varias pruebas, ensayos y evaluaciones que sirvieron como guía para testear el rendimiento de los distintos modelos generados en cuanto a calidad y a su nivel predictivo potencial.

En este sentido, el RMSE y el R2 cada cual, midiendo aspectos diferentes de la precisión de los modelos entregaron resultados disímiles. Por un lado, el RMSE presentó para la mayoría de los modelos valores superiores a 1.5. En el caso de varios meses para los modelos construidos mediante los métodos de IDW y de regresión lineal esta situación empeora superando la barrera de los 2 puntos. El caso del R2 fue más favorable con valores en su mayoría superiores al 0.90, siendo el valor más bajo un 0.83 para el mes de marzo bajo el método de regresión lineal.

Para el caso concreto del vecino más cercano los valores más altos de RMSE para cada mes se encontraron para el caso de k=10 (calculo efectuado con los 10 vecinos más cercanos). Siendo el caso de septiembre el más pronunciado alcanzando un valor de 1,68. En general para el método en cuestión los valores de RMSE se encuentran entre 1.6 y 1.4 que tomando en cuenta los valores de temperatura que van desde los -9 °C a los 23 °C no es un valor tan alto. En este apartado destaca el mes de junio con un resultado para k=3 de un valor de 1.43. Considerando el R2 el método del vecino más próximo entrega muy buenos resultados en todos los meses y en cada uno de los distintos k.

El método del IDW, aunque muestra valores altos para todos los meses para el indicador del R2, sus resultados para el RMSE no son de muy buena calidad, por lo que podría existir en sus modelos mensuales una diferencia significativa entre los valores predichos (interpolados) y los valores observados (estaciones). El caso de la regresión lineal es parecido al del IDW en cuanto a sus comportamientos para las métricas de RMSE y el R2. Particularmente, en este caso los modelos de los meses de marzo y diciembre presentan los peores resultados. Aunque los valores de RMSE para junio y septiembre también son relativamente altos.

En cuento a la varianza en los modelos de regresión se obtuvo que los valores predichos mas dispersos con respecto a los esperados se encontraron en los modelos de marzo y junio. Por su parte, el mes de septiembre con una varianza de entre 3.1 y 5.7 estaría siendo el modelo más confiable en cuanto a que sus valores predichos se encontrarían en bastante concordancia con los valores esperados.

05. Conclusión:

En conclusión, luego de haber observado y evaluado el rendimiento de los diferentes métodos de interpolación utilizados se puede afirmar que el modelo predictivo más confiable estadísticamente, considerando el conjunto de métricas disponibles, corresponde al del vecino más próximo con un valor de k=3.

Cartografía del mejor modelo predictivo construido con el método del vecino más próximo (k=3)

l’esercizio è finito