Taller 2: Métodos de predicción espacial determinísticos y de regresión en Región de Coquimbo
Author
Victor Contreras, Cristian Aguirre y Christian Chacon
Published
October 16, 2023
Taller 2: Métodos de predicción espacial determinísticos y de regresión para predicción de temperatura en la Región de Coquimbo
1. Introducción
En este taller, se abordará la predicción espacial de la temperatura en la Región de Coquimbo, utilizando datos raster y vectoriales previamente preparados. Los predictores utilizados para este análisis incluyen la elevación, la pendiente (slope), la distancia a la costa, la orientación de la ladera (aspect), el índice de vegetación de diferencia normalizada (NDVI) y la temperatura superficial del suelo (LST). Por otro lado, los datos de temperatura observada provienen de las estaciones meteorológicas de la red Agromet, correspondientes a los meses de enero, junio, octubre y diciembre de 2023.
El objetivo de este taller es aplicar tres métodos de interpolación determinísticos para generar predicciones de la temperatura a nivel regional: el método del vecino más cercano (KNN), el inverso a la distancia pondera (IDW) y la regresión lineal (Trend Surface Fitting - TSF). Estos métodos se implementarán en R utilizando la librería {gstat}, y su desempeño será evaluado mediante las métricas de error cuadrático medio (RMSE) y coeficiente de determinación (R²), a fin de determinar cuál es el método más adecuado para la predicción de la temperatura en la región.
2. Objetivos
Objetivo General:
Realizar la predicción espacial de la temperatura en los meses de enero, junio, octubre y diciembre de 2023, utilizando métodos de interpolación determinísticos implementados en R con la librería {gstat}, y determinar cuál es el método más preciso a partir de las métricas RMSE y R².
Objetivos específicos:
1.- Aplicar el método del vecino más cercano (KNN) con tres valores de k diferentes y evaluar el modelo mediante validación cruzada Leave-One-Out, obteniendo las métricas correspondientes.
2.- Implementar el método de interpolación inversa ponderada por la distancia (IDW) y evaluar su desempeño utilizando el mismo proceso de validación cruzada.
3.- Aplicar un modelo de regresión lineal utilizando los seis predictores raster y evaluar su rendimiento mediante validación cruzada Leave-One-Out.
4.- Generar mapas espaciales de predicción de la temperatura para cada uno de los métodos y para los cuatro meses seleccionados.
5.- Comparar las métricas obtenidas para cada método utilizando herramientas gráficas en R y seleccionar el método de interpolación más adecuado para la predicción de la temperatura en la región.
3. Desarrollo
3.1 Librería y preparación de datos
Se preparan las librerias a utilizar.
library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.1 ✔ tibble 3.2.1
✔ lubridate 1.9.3 ✔ tidyr 1.3.1
✔ purrr 1.0.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(terra)
terra 1.7.78
Adjuntando el paquete: 'terra'
The following object is masked from 'package:tidyr':
extract
library(tmap)
Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
remotes::install_github('r-tmap/tmap')
library(sf)
Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
library(viridis)
Cargando paquete requerido: viridisLite
library(viridisLite)library(spData)
To access larger datasets in this package, install the spDataLarge
package with: `install.packages('spDataLarge',
repos='https://nowosad.github.io/drat/', type='source')`
library(raster)
Cargando paquete requerido: sp
Adjuntando el paquete: 'raster'
The following object is masked from 'package:dplyr':
select
library(purrr)library(gstat)library(stars)
Cargando paquete requerido: abind
library(raster)library(ggplot2)library(patchwork)
Adjuntando el paquete: 'patchwork'
The following object is masked from 'package:raster':
area
The following object is masked from 'package:terra':
area
Se preparan los datos geoespaciales necesarios para el análisis, incluyendo el polígono que delimita la Región de Coquimbo, los datos puntuales correspondientes a las estaciones meteorológicas con sus seis predictores asociados, y los datos ráster de estos predictores. Se garantiza que todos los datos estén en el mismo Sistema de Coordenadas de Referencia (SCR) adecuado para el área de estudio. Posteriormente, se crean subconjuntos de los datos puntuales y ráster, diferenciados por mes, para su uso en el proceso de interpolación.
Se busca visualizar la distribución espacial de las estaciones meteorológicas y las temperaturas registradas en cada mes. Para asegurar una visualización coherente y comparable entre los meses, se define un rango de temperatura mínimo y máximo, lo que permite normalizar los valores de temperatura. De esta forma, se asegura que las visualizaciones de cada mes se realicen bajo la misma escala, facilitando su interpretación.
# Definir los límites y los intervalos para la paletatemp_min <--4temp_max <-24breaks <-seq(temp_min, temp_max, length.out =7) # 5 intervalos => 6 puntos de ruptura
Se confeciona una función para optimizar el proceso de generación de mapas
# Lista con los datos por mes y sus nombresmeses <-list("Enero"= data_ene_df,"Junio"= data_jun_df,"Octubre"= data_oct_df,"Diciembre"= data_dic_df)# Generar los mapasmapas <-lapply(names(meses), function(mes) {stations_maps(meses[[mes]], mes, reg_coq_utm, breaks)})
Y luego se obtienen los mapas para los meses deseados
# Mostrar los mapastmap_arrange(mapas[[1]], mapas[[2]], ncol =2)
Warning: Breaks contains positive and negative values. Better is to use
diverging scale instead, or set auto.palette.mapping to FALSE.
Warning in CPL_crs_from_input(x): GDAL Message 1: +init=epsg:XXXX syntax is
deprecated. It might return a CRS with a non-EPSG compliant axis order.
Some legend labels were too wide. These labels have been resized to 0.63, 0.57, 0.57, 0.57. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
Warning: Breaks contains positive and negative values. Better is to use
diverging scale instead, or set auto.palette.mapping to FALSE.
Some legend labels were too wide. These labels have been resized to 0.63, 0.57, 0.57, 0.57. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
Enero:Se observa una tendencia general de temperaturas cálidas en los valles centrales, con una orientación norte-sur, registrando valores entre 19 y 24 °C. En la costa, las temperaturas son intermedias, fluctuando entre 10 y 19 °C, mientras que en la cordillera principal se aprecian temperaturas significativamente más bajas, entre 0 y 10 °C. Esta distribución está fuertemente influenciada por factores topográficos y la proximidad al mar, ya que las zonas costeras suelen ser más frías que los valles centrales. Esto es esperable durante el pleno verano, cuando predominan las temperaturas altas, aunque se observan algunos valores anómalamente bajos en los valles centrales.
Junio: Durante esta estación más fría, en pleno otoño, predominan las temperaturas bajas en la costa, entre 10 y 14 °C. En los valles centrales, algunas zonas aisladas alcanzan temperaturas de entre 14 y 19 °C, mientras que en la alta cordillera se registran temperaturas frías, como es característico en esa región.
# Mostrar los mapas tmap_arrange(mapas[[3]], mapas[[4]], ncol =2)
Warning: Breaks contains positive and negative values. Better is to use
diverging scale instead, or set auto.palette.mapping to FALSE.
Some legend labels were too wide. These labels have been resized to 0.63, 0.57, 0.57, 0.57. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
Warning: Breaks contains positive and negative values. Better is to use
diverging scale instead, or set auto.palette.mapping to FALSE.
Some legend labels were too wide. These labels have been resized to 0.63, 0.57, 0.57, 0.57. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
Octubre:En primavera, las temperaturas son intermedias en la costa, de 10 a 14 °C, mientras que en los valles centrales se registran temperaturas intermedia-altas, entre 14 y 19 °C. En la alta cordillera, las temperaturas siguen siendo bajas.
Diciembre: Con el inicio del verano, se observa un aumento general de las temperaturas. Las zonas costeras presentan temperaturas intermedias, mientras que los valles centrales alcanzan temperaturas altas. En la alta cordillera, las temperaturas se mantienen bajas.
En cuanto a la distribución de las estaciones meteorológicas, se aprecia una marcada ausencia de estaciones en el norte de la Región de Coquimbo, además de un agrupamiento en ciertos sectores. Los valles centrales concentran una cantidad importante de estaciones con mediciones de temperatura, al igual que algunas zonas costeras. Sin embargo, en la alta cordillera, las estaciones son escasas y, en algunos meses como diciembre, no se dispone de datos. Esta distribución no uniforme podría impactar en la precisión de algunos métodos de interpolación, aunque a simple vista la temperatura tiene una fuerte correlación con estructuras geomorfológica como lo son la cordillera de la costa, valles centrales y cordillera principal, como también longitudinal relacionado con la distancia de la costa hacia el continente.
Como parte de un análisis exploratorio de los datos, se contruye histogramas de temperatura para cada mes con el fin de ver la distribución de la temperatura.
# Gráfica para Eneroplot_ene <-ggplot(data_ene_df, aes(x = temp)) +geom_histogram(aes(y =after_stat(density)), binwidth =2, fill ="skyblue", color ="black", alpha =0.7) +geom_density(color ="blue", size =1) +labs(title ="Enero", x ="Temperatura (°C)", y ="Densidad") +theme_minimal()
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
# Gráfica para Junioplot_jun <-ggplot(data_jun_df, aes(x = temp)) +geom_histogram(aes(y =after_stat(density)), binwidth =2, fill ="lightgreen", color ="black", alpha =0.7) +geom_density(color ="green", size =1) +labs(title ="Junio", x ="Temperatura (°C)", y ="Densidad") +theme_minimal()# Gráfica para Octubreplot_oct <-ggplot(data_oct_df, aes(x = temp)) +geom_histogram(aes(y =after_stat(density)), binwidth =2, fill ="lightpink", color ="black", alpha =0.7) +geom_density(color ="red", size =1) +labs(title ="Octubre", x ="Temperatura (°C)", y ="Densidad") +theme_minimal()# Gráfica para Diciembreplot_dic <-ggplot(data_dic_df, aes(x = temp)) +geom_histogram(aes(y =after_stat(density)), binwidth =2, fill ="purple", color ="black", alpha =0.7) +geom_density(color ="purple", size =1) +labs(title ="Diciembre", x ="Temperatura (°C)", y ="Densidad") +theme_minimal()# Combinar las gráficas en una sola figura, 2x2(combined_plot <- (plot_ene | plot_jun) / (plot_oct | plot_dic))
En los meses de enero y diciembre se concentran las temperaturas más altas. Sin embargo, enero presenta una cola hacia valores negativos, debido a las mediciones realizadas en la zona de alta cordillera. En diciembre, no se disponen de datos para esta estación de alta montaña. En cuanto a la morfología de los histogramas, enero presenta una distribución bimodal, con un pico en torno a los 15 °C y otro cerca de los 20 °C. Esta bimodalidad podría estar relacionada con cambios bruscos de temperatura provocados por las condiciones climáticas y geográficas de la región. En contraste, diciembre tiene una distribución más uniforme, con un pico significativo alrededor de los 18 °C y menor variabilidad extrema.
Por su parte, junio y octubre muestran distribuciones asimétricas hacia la derecha. En junio, la mayoría de las temperaturas se concentran alrededor de los 10 °C, mientras que en octubre el pico principal se sitúa en torno a los 15 °C.
En efecto se puede decir que la asimetría observada en algunos meses sugiere una variabilidad climática y geográfica, mientras que la bimodalidad en enero podría reflejar condiciones climáticas muy contrastadas, como las diferencias entre la zona costera y los valles centrales. Finalmente, la uniformidad de diciembre podría indicar una mayor estabilidad en las temperaturas, aunque es importante señalar que faltan mediciones en las áreas de alta montaña, lo que podría influir en la interpretación de los resultados.
A continuación, se desea modelar la tempertura en la región de Coquimbo a partir de tres métodos no geoestadísticos.
3.3. Vecino más cercano (K-NN)
El algoritmo del vecino más cercano (K-Nearest Neighbor - KNN) es un método de clasificación que asigna una categoría a un objeto basado en sus vecinos más cercanos en el espacio de características. Para cada nuevo punto de datos, el algoritmo identifica los k puntos más cercanos (vecinos) del conjunto de datos de entrenamiento y lo clasifica en función de la clase mayoritaria entre estos vecinos. La distancia entre los puntos puede calcularse utilizando diferentes métricas, siendo la más común la distancia Euclidiana (Adrian et al., 2019).
Este método tiene como característica que la elección del valor k es crucial: un valor de k bajo puede hacer que el modelo sea muy sensible al ruido, mientras que un valor de k alto puede causar una sobre-generalización, donde se pierde detalle.
Para la interpolación mediante este método se crear una grilla base para generar los mapas de interpolación con una resolución espacial de 500 metros y con el sistsema de coordenadas de referencia EPSG:32719. Esta grilla será reutilizada para los otros métodos de interpolación.
Se establece una semilla para garantizar que los resultados sean reproducibles y se asegura que la grilla tenga el mismo sistema de referencia de coordenadas que los datos de las estaciones. Además, se seleccionan los valores de K para la interpolación, eligiéndose en este caso 5, 10 y 15 como valores representativos.
# Lanzar semillaset.seed(432)#forzar que los src sean igualst_crs(grilla_st) <-crs(data_ene_df)# Los k a uilizarkas <-c(5,10,15)
El desempeño de las intepolaciones serán evaluados con las métricas del error cuadrático medio (RMSE) y coeficiente de determinación (R²).
El RMSE (Raíz del Error Cuadrático Medio) es una métrica que compara los valores observados (obs) con los predichos (pred) en un mismo punto o registro. La diferencia entre estos valores se denomina residuo. Para calcular el RMSE, se elevan al cuadrado todos los residuos, se suman, y luego se dividen por el número total de observaciones. Finalmente, se obtiene la raíz cuadrada de este valor. El RMSE proporciona una medida del error promedio en las predicciones, penalizando más los errores grandes. Mientras más pequeño es el RMSE, mejor en la predicción obtenida, y la unidades dependerá de la variable a estudiar.
\[
RMSE = \sqrt{\frac{\sum (obs - pred)^2}{N}}
\]
El R² (coeficiente de determinación) es una métrica sin unidades que mide la proporción de la varianza en los datos observados que puede ser explicada por los datos predichos. Su valor oscila entre 0 y 1, donde un valor cercano a 1 indica que el modelo tiene una mejor capacidad para explicar la varianza del objetivo a partir de las variables predictoras. En otras palabras, un R² más alto refleja una mayor precisión en la predicción del modelo.
Es importante destacar que un R² bajo no necesariamente implica que el modelo sea malo; esto depende del contexto del problema y los datos utilizados. Si un modelo de regresión lineal produce resultados con un R² bajo, mientras que otros algoritmos muestran mejor rendimiento, podría simplemente significar que los datos no se ajustan bien a un modelo lineal y que requieren un enfoque más complejo o no lineal.
Metodo Leave-one-out (LOOCV)
El método Leave-One-Out se utiliza para estimar el rendimiento de los algoritmos de aprendizaje automático cuando se aplican para hacer predicciones sobre datos no utilizados en el entrenamiento. Aunque este método es costoso en términos de rendimiento computacional, ofrece una estimación confiable e imparcial del rendimiento del modelo. Es ideal para conjuntos de datos pequeños.
Este método funciona de la siguiente manera: primero, utiliza todos los puntos de entrada para estimar los parámetros de un modelo de interpolación. Luego, se elimina un único punto de entrada, y el resto de los puntos se utilizan para predecir el valor en la ubicación del punto oculto. El valor predicho se compara con el valor medido. Posteriormente, el punto oculto se reincorpora al conjunto de datos, y se repite el proceso ocultando y prediciendo un punto diferente en cada iteración. Este ciclo se repite para todos los puntos de entrada (ESRI, s.f.).
A continuación, se aplica el método de LOOCV para cada uno de los meses con k = 5, 10 y 15.
Enero
# Calculo predicciones para cada Kas, aplicando la interpolación Vecino más cercano### Enero loocv_nn_ene <-suppressWarnings( purrr::map_df(1:nrow(data_ene_df), \(i) { # corre ciclo desde 1 hasta el número de observaciones purrr::map_df(kas, \(k) { # aplica el inverso a la distancia a todas las observaciones menos a la i temp_idw <-idw(temp ~1, data_ene_df[-i,], grilla_st,nmax=k) temp_idw <- terra::rast(temp_idw)# extrae el valor de la temp interpolada en df <- terra::extract(temp_idw, data_ene_df[i,]) df |>mutate(k = k) }) }))
# Tabla de comparación de métricas# Función para crear gráficos de métricaslibrary(patchwork)crear_metric_plot_knn <-function(data_cv_sum, mes) { data_cv_sum |>pivot_longer(-k) |>ggplot(aes(k, value)) +geom_col(alpha = .7) +coord_flip() +facet_grid(. ~ name, scales ='free') +theme_bw() +ggtitle(paste("Métricas de", mes,"del KNN"))}# Crear gráficos individuales para cada mesplot_ene <-crear_metric_plot_knn(data_nn_ene_cv_sum, "Enero")plot_jun <-crear_metric_plot_knn(data_nn_jun_cv_sum, "Junio")plot_oct <-crear_metric_plot_knn(data_nn_oct_cv_sum, "Octubre")plot_dic <-crear_metric_plot_knn(data_nn_dic_cv_sum, "Diciembre")# Combinar los gráficos en una cuadrículacombined_plot <- plot_ene + plot_jun + plot_oct + plot_dic +plot_layout(ncol =2) # Configura que haya 2 columnas# Mostrar el gráfico combinadocombined_plot
Interpretación de las métricas con KNN
- Enero: se observa que el valor de RMSE es más bajo cuando k = 10 (3.42), lo cual indica un mejor ajuste del modelo, ya que los errores de predicción son menores. Sin embargo, el valor de R² es muy bajo en todos los casos, lo que sugiere que la interpolación tiene una capacidad limitada para explicar la variabilidad de los datos. Un valor de k más alto (15) aumenta el RMSE, lo que puede indicar un sobreajuste o que la interpolación promedia demasiado la información de los vecinos.
- Junio: los valores de RMSE son similares entre k = 5 y k = 15, pero un poco más altos cuando k = 10. Sin embargo, en este caso, k = 10 muestra un R² algo superior (0.0276), lo cual, aunque sigue siendo bajo, sugiere una leve mejora en la capacidad de predicción en comparación con otros valores de k. Aun así, la capacidad explicativa sigue siendo baja, lo que sugiere que la interpolación k-NN no captura bien la variabilidad de los datos.
- Octubre: el valor de RMSE es más bajo con k = 5 (3.47), lo que indica un mejor desempeño en términos de error de predicción. Sin embargo, los valores de R² son extremadamente bajos, lo que implica que la interpolación apenas explica la variabilidad de la temperatura. Al aumentar k, el error aumenta considerablemente, lo cual sugiere que un mayor número de vecinos no ayuda a mejorar el ajuste en este caso.
- Diciembre: el RMSE más bajo se obtiene con k = 5 (3.15), mientras que los valores de R² siguen siendo bajos para todos los valores de k. La ligera diferencia en el RMSE entre k = 5 y k = 10 sugiere que el modelo de interpolación no se ve fuertemente afectado por un aumento en el número de vecinos, pero la capacidad explicativa de la interpolación sigue siendo limitada.
A continuación de presentan los mapas de las interpolaciones para cada mes y con cada valor de k.
meses_map_knn <-list("Enero"=list(data = data_ene_df, data_cv_sum = data_nn_ene_cv_sum),"Junio"=list(data = data_jun_df, data_cv_sum = data_nn_jun_cv_sum),"Octubre"=list(data = data_oct_df, data_cv_sum = data_nn_oct_cv_sum),"Diciembre"=list(data = data_dic_df, data_cv_sum = data_nn_dic_cv_sum))# Generar y mostrar mapas para cada mes y cada valor de Kmaps <-list()for (mes innames(meses_map_knn)) { data_info <- meses_map_knn[[mes]]for (k in kas) { map <-generar_knn_map(data_info$data, grilla_st, reg_coq_utm, k, data_info$data_cv_sum, mes) maps[[paste0(mes, "_K", k)]] <- map }}
Warning in asMethod(object): complete map seems to be NA's -- no selection was
made
[inverse distance weighted interpolation]
Warning in asMethod(object): complete map seems to be NA's -- no selection was
made
Warning in asMethod(object): complete map seems to be NA's -- no selection was
made
[inverse distance weighted interpolation]
Warning in asMethod(object): complete map seems to be NA's -- no selection was
made
Warning in asMethod(object): complete map seems to be NA's -- no selection was
made
[inverse distance weighted interpolation]
Warning in asMethod(object): complete map seems to be NA's -- no selection was
made
Warning in asMethod(object): complete map seems to be NA's -- no selection was
made
[inverse distance weighted interpolation]
Warning in asMethod(object): complete map seems to be NA's -- no selection was
made
Warning in asMethod(object): complete map seems to be NA's -- no selection was
made
[inverse distance weighted interpolation]
Warning in asMethod(object): complete map seems to be NA's -- no selection was
made
Warning in asMethod(object): complete map seems to be NA's -- no selection was
made
[inverse distance weighted interpolation]
Warning in asMethod(object): complete map seems to be NA's -- no selection was
made
Warning in asMethod(object): complete map seems to be NA's -- no selection was
made
[inverse distance weighted interpolation]
Warning in asMethod(object): complete map seems to be NA's -- no selection was
made
Warning in asMethod(object): complete map seems to be NA's -- no selection was
made
[inverse distance weighted interpolation]
Warning in asMethod(object): complete map seems to be NA's -- no selection was
made
Warning in asMethod(object): complete map seems to be NA's -- no selection was
made
[inverse distance weighted interpolation]
Warning in asMethod(object): complete map seems to be NA's -- no selection was
made
Warning in asMethod(object): complete map seems to be NA's -- no selection was
made
[inverse distance weighted interpolation]
Warning in asMethod(object): complete map seems to be NA's -- no selection was
made
Warning in asMethod(object): complete map seems to be NA's -- no selection was
made
[inverse distance weighted interpolation]
Warning in asMethod(object): complete map seems to be NA's -- no selection was
made
Warning in asMethod(object): complete map seems to be NA's -- no selection was
made
[inverse distance weighted interpolation]
Warning in asMethod(object): complete map seems to be NA's -- no selection was
made
Al comparar la distribución espacial de la temperatura, se observa que, a medida que el valor de k aumenta, la interpolación se vuelve más suave y la variabilidad espacial disminuye. Esto se debe probablemente a que, al considerar un mayor número de vecinos, las diferencias locales se suavizan.
Con k = 5, la predicción de la temperatura presenta una mayor variabilidad, con transiciones muy marcadas entre valores extremos. En el sector norte, donde hay una menor densidad de datos, se generan artefactos de discontinuidad que se atenúan conforme se incrementa k. Esto es observado en todas las demás estaciones.
En cuanto a la distribución de la temperatura, se observa que las temperaturas más cálidas se concentran en la parte central de la Región de Coquimbo, siguiendo una tendencia generalizada de norte a sur. En la zona costera, las temperaturas son más templadas, posiblemente debido a la influencia moderadora del océano. En contraste, hacia el este, en la zona de alta cordillera, las temperaturas son más bajas, lo que probablemente esté relacionado con las condiciones topográficas de mayor altitud.
Se observa una disminución generalizada de la temperatura, lo cual está asociado a una estación más fría. En los valles centrales de la Región de Coquimbo, se registran algunos valores anómalos de alta temperatura. Las temperaturas más bajas se concentran en la zona costera, con un descenso aún más marcado en el sector cordillerano.
Se observa una tendencia de temperatura similar a la de los mapas anteriores, pero con un aumento en los valles centrales en comparación con el mes de junio. Esto coincide con el histograma presentado previamente, que muestra una distribución de temperaturas similar, con valores extremos de bajas temperaturas en la cordillera principal. Sin embargo, se nota una mayor concentración de mediciones cercanas a los 15 °C, en contraste con las temperaturas alrededor de los 10 °C observadas en junio.
En general, este mes presenta una distribución de temperatura más uniforme en comparación con los meses anteriores. Las temperaturas mínimas ya no caen por debajo de los 0 °C, situándose en un rango de 10 a 22 °C, lo que dificulta un poco la comparación directa con los mapas previos. Aun así, se observa la misma tendencia de variación de la temperatura a medida que se aleja de la costa: temperaturas medias cerca de la costa, más altas en los valles centrales, y más frías en la zona de la alta cordillera.
3.4. Inverso a la distancia (IDW)
La interpolación por pesos inversos de la distancia (IDW) es un método de interpolación espacial que estima el valor de una variable en ubicaciones no muestreadas a partir de puntos de muestreo cercanos. El principio subyacente es que los puntos más cercanos al sitio de estimación tienen mayor influencia sobre el valor estimado que los puntos más distantes. La función de ponderación es inversamente proporcional a la distancia, lo que permite que los puntos más próximos contribuyan de manera más significativa al valor interpolado (Burrough & McDonnell, 1998).
La fórmula que se presenta a continuación es una expresión para calcular el peso ( lambda(s_0) ) en la interpolación por Pesos Inversos de la Distancia (IDW). Este peso se utiliza para determinar cuánto influye un punto ( s_i ) con un valor conocido ( Z(s_i) ) en la estimación de un valor desconocido en el punto ( s_0 ). La fórmula para ( lambda(s_0) ) es:
Donde: - ( lambda(s_0) ): El peso asignado a la influencia de cada punto ( s_i ) en el valor interpolado en la ubicación ( s_0 ). - ( d(s_0, s_i) ): La distancia entre el punto de estimación ( s_0 ) y el punto con valor conocido ( s_i ). - ( beta ): El parámetro que controla el grado de influencia de la distancia. Valores mayores de ( beta ) reducen drásticamente la influencia de puntos más lejanos (comúnmente ( beta = 2 ) para un peso inverso cuadrático). - ( n ): El número total de puntos con mediciones conocidas que se usan para el cálculo.
Esta fórmula asigna un peso proporcional a cada valor conocido ( Z(s_i) ) en función de su distancia al punto ( s_0 ) que se está estimando. Cuanto más cerca esté ( s_i ) de ( s_0 ), mayor será su contribución a la estimación del valor en ( s_0 ). Todos los pesos se normalizan dividiendo cada uno entre la suma de los pesos totales, de modo que los pesos siempre suman 1, lo que garantiza que la interpolación sea una combinación ponderada de los valores conocidos.
Validacion Cruzada LOOCV para IDW
Enero
# Método Leave-one-out para IDWloocv_idw_ene <-suppressWarnings( # Suprimir la alerta de NA en zonas fuera de interpolación purrr::map_df(1:nrow(data_ene_df), \(i) { purrr::map_df(1, \(k) { temp_idw <-idw(temp ~1, data_ene_df[-i,], grilla_st) temp_idw <- terra::rast(temp_idw) df <- terra::extract(temp_idw, data_ene_df[i,]) df |>mutate(k = k) }) }))
- Enero: El valor del RMSE indica un error promedio de 3.34 unidades de temperatura en las predicciones de enero, lo cual sugiere una variabilidad moderada entre los valores observados y los estimados. El R² de 0.27 indica que el modelo explica el 27% de la variabilidad en los datos, lo cual muestra un ajuste relativamente bajo, reflejando la dificultad de captar patrones espaciales complejos con IDW en este mes, sin embargo, en comparación con KNN para ese mismo mes es muy superior.
- Junio: presenta un menor error promedio (RMSE de 2.92), lo que indica que las predicciones son un poco más precisas que en enero. Además, el R² de 0.36 muestra que el modelo explica un 36% de la variabilidad de los datos, lo que sugiere una mejor captura de los patrones de temperatura en este mes, aunque todavía hay una considerable cantidad de variación no explicada.
- Octubre: tiene el mayor RMSE (3.56), lo que indica un mayor error de predicción comparado con los otros meses, y un R² de 0.28, similar al de enero. Esto sugiere que las características espaciales de las temperaturas en octubre son más difíciles de modelar con IDW, lo cual podría deberse a una mayor variabilidad espacial o a la presencia de valores extremos que el método no logra captar adecuadamente.
- Diciembre: Este mes tiene el menor RMSE (2.48), lo que indica una mejor precisión en las predicciones en comparación con los otros meses. Sin embargo, el R² es el más bajo (0.23), lo que significa que el modelo solo explica el 23% de la variabilidad de los datos. Esto puede indicar que, aunque el error promedio sea menor, el IDW en diciembre no captura bien la estructura espacial de la variabilidad de las temperaturas.
En general, se aprecia que el IDW explica mejora la variabilidad de los datos en comparación al KNN, pero aun sigue siendo bajo.
# Mostrar todos los mapas juntostmap_arrange(map_temp_ene_idw, map_temp_jun_idw, map_temp_oct_idw, map_temp_dic_idw, ncol=2 , nrow =2)
- Enero: La predicción de temperatura muestra valores más cálidos en los valles centrales, con una transición de temperaturas hacia el sur. Las áreas costeras y cordillera princiapl presentan temperaturas más bajas. La baja variabilidad explicada por el modelo (R²) sugiere que hay factores espaciales no capturados adecuadamente por la interpolación.
- Junio: Se observa una distribución más homogénea de temperaturas frías, típica de una estación invernal, con un enfriamiento generalizado en la costa y las áreas de alta montaña. Sin embargo, algunas áreas de los valles centrales presentan valores ligeramente más altos de temperatura.El mejor ajuste de R² entre los meses (0.358) indica que el IDW en junio capta mejor la variabilidad espacial en comparación con otros meses.
- Octubre: La distribución de temperatura en octubre muestra una mayor variabilidad, con un patrón de temperaturas más altas en los valles centrales y un enfriamiento hacia la costa y las zonas más elevadas.Sin embargo, la interpolación presenta un mayor error (RMSE de 3.56) y un R2 moderado, lo que sugiere que hay áreas con variaciones de temperatura difíciles de capturar por el IDW, posiblemente por cambios bruscos entre regiones.
- Diciembre: el mapa muestra un incremento general de las temperaturas, especialmente en las áreas centrales. Las temperaturas más bajas se mantienen en la costa y la alta montaña, aunque ya no se observan temperaturas extremas bajo los 10 °C. A pesar de tener el menor error (RMSE), el R² indica que el modelo sigue explicando una fracción limitada de la variabilidad espacial de la temperatura.
3.5. Modelamiento lineal: Trend Surface Fitting
El Trend Surface Fitting (ajuste de superficie de tendencia) es una técnica de interpolación que se usa para representar patrones generales en datos espaciales. Se basa en crear una ecuación matemática que aproxima los valores de un área de estudio. Esta ecuación puede ser simple, como un plano inclinado, o más compleja, utilizando polinomios de mayor orden para capturar variaciones más detalladas (Adrian et al., 2019).
\[
Z(s) = f(x, y) + \epsilon
\]
El objetivo es ajustar una fórmula que describa cómo varía un fenómeno (como la temperatura, la altitud, etc.) en función de su ubicación espacial. Este ajuste puede ser global, cuando se aplica a toda el área de estudio, o local, cuando se divide el área en partes más pequeñas para capturar mejor las variaciones.
En este caso, se desea modelar la varibale temperatura a partir de un polinomio que considere la altura (DEM), la orientación de las laderas (Aspect), la pendiente (Slope), la distancia de la costa (DC), la temperatura superficial del terreno (LST) y el Indice de diferencia normalizada de vegetación (NDVI)
Es importante tener en cuenta factores como los efectos de borde, que pueden afectar los resultados en los límites del área de estudio, y elegir el tipo de fórmula adecuado según la complejidad de los datos. En este caso, como límite se considera la Región de Coquimbo y las variables. predictorias de la formula anterior.
pred_mes rmse r2
1 Enero ml 2.419045 0.6117440
2 Junio ml 1.587847 0.8030499
3 Octubre ml 2.181483 0.7284522
4 Diciembre ml 1.998125 0.5094585
Gráficos comparativos
data_lm_cv_all |>pivot_longer(-pred_mes) |>mutate(pred_mes =factor(pred_mes, levels =c("Enero ml", "Junio ml", "Octubre ml", "Diciembre ml"))) |>ggplot(aes(pred_mes, value, fill = pred_mes)) +geom_col(alpha = .7) +coord_flip() +facet_grid(. ~ name, scales ='free') +scale_fill_manual(name ="Meses de Predicción", labels =c("Enero ml"="Enero", "Junio ml"="Junio", "Octubre ml"="Octubre", "Diciembre ml"="Diciembre"), # Etiquetas de la leyendavalues =c("Enero ml"="blue", "Junio ml"="green", "Octubre ml"="red", "Diciembre ml"="purple") # Colores específicos ) +labs(title ="Métricas de Modelamiento Lineal por Mes", x ="Meses de Modelamiento", y ="Valor de Métrica" ) +theme_bw() +theme(plot.title =element_text(hjust =0.5, size =14, face ="bold"), legend.position ="right" )
- Enero: Según la métrica R2, este es explica un 61.2% de la variabilidad de la temperatura. El RMSE relativamente bajo (2.42) indica que el ajuste logra predecir la temperatura con una buena precisión, aunque aún hay un margen de error debido a la variabilidad no capturada por el modelo.
- Junio: se observa el mejor rendimiento del ajuste con un R2 de 0.803, lo que significa que el modelo captura el 80.3% de la variabilidad espacial de la temperatura. El RMSE más bajo (1.59) sugiere una excelente precisión en la predicción, lo que puede deberse a una distribución más homogénea de la temperatura en este mes de invierno, facilitando el ajuste del modelo.
- Octubre: el R2 de 0.728 indica que el modelo explica el 72.8% de la variabilidad de la temperatura. El RMSE de 2.18 muestra un buen ajuste, aunque ligeramente inferior al de junio, lo que podría ser resultado de una mayor variabilidad climática en la primavera y cambios más marcados en la temperatura.
- Diciembre: Para este mes, el modelo presenta un R2 de 0.509, lo que significa que solo el 50.9% de la variabilidad de la temperatura es explicada por el ajuste, siendo el mes con el menor coeficiente de determinación. Esto sugiere que en este mes, con el método de Trend Surface Fitting no es posible de captura de forma óptima la variabilidad espacial de la temperatura, posiblemente debido a una mayor complejidad en los patrones de temperatura. Sin embargo, el RMSE de 2.00 sigue siendo relativamente bajo, lo que indica que el ajuste tiene una precisión moderada.
Preparar funcion para interpolación y confeccion de mapas
# Función para generar mapas de predicción y varianzacrear_mapas_lm <-function(pred_ras, df, ndvi, lst, mes, pred_mes_name) {# Crear predictores pred_ras_st <- pred_ras |>st_as_stars() |>split('band')# Ajustar el modelo de regresión lineal kriging formula <-as.formula(paste("temp ~ dem + slope + dist_costa + aspect +", ndvi, "+", lst)) temp_lm <-krige(formula, df, pred_ras_st) temp_lm <- temp_lm |>rast() |>mask(vect(reg_coq_utm)) |>trim()# Extraer capas de predicción y varianza predicciones <-subset(temp_lm, "var1.pred_var1.pred") varianza <-subset(temp_lm, "var1.var_var1.var")# Crear el mapa de predicción map_predicciones <-tm_shape(predicciones) +tm_raster(style ='cont', palette = viridis::viridis(20),title =paste("Temp. Pred.", mes, "°C")) +tm_layout(main.title =paste("Predicción de la Temperatura", mes),main.title.size =1.1, bg.color ="gray90", frame =TRUE, frame.lwd =2, frame.double.line =FALSE) +tm_shape(reg_coq_utm) +tm_borders() +tm_compass(size =3, type ="arrow", position =c("right", "top")) +tm_grid(projection ="+init=epsg:4326", lines =TRUE, col ="gray", alpha =0.5) +tm_legend(outside =TRUE, bg.color ="white", frame =TRUE) +tm_shape(df) +tm_dots() +tm_credits(paste0("RMSE: ", round(data_lm_cv_all |>filter(pred_mes == pred_mes_name) |>pull(rmse), 2), "\nR²: ", round(data_lm_cv_all |>filter(pred_mes == pred_mes_name) |>pull(r2), 3)),position =c("right", "bottom"),bg.color ="white", bg.alpha =1)# Crear el mapa de varianza map_varianza <-tm_shape(varianza) +tm_raster(style ='cont', palette = viridis::viridis(20),title ="Varianza") +tm_layout(main.title =paste("Varianza de la Predicción", mes),main.title.size =1.1, bg.color ="gray90", frame =TRUE, frame.lwd =2, frame.double.line =FALSE) +tm_shape(reg_coq_utm) +tm_borders() +tm_compass(size =3, type ="arrow", position =c("right", "top")) +tm_grid(projection ="+init=epsg:4326", lines =TRUE, col ="gray", alpha =0.5) +tm_legend(outside =TRUE, bg.color ="white", frame =TRUE)# Mostrar ambos mapas en paralelotmap_arrange(map_predicciones, map_varianza, ncol =2)}
Con respecto a la predicción de temperatura para enero muestra una distribución coherente a los métodos de inteporlación anteriores, pero al incorporar otras variables predictorias como características geográficas, se presenta una distribución más coherente de la variable, con temperaturas más cálida en los valles centrales, y más fría en los sectores costero y aun más fria en la alta cordillera. Sin embargo, la variabilidad de la predicción es alta en las zonas con menor cobertura de estaciones, como lo es en el sector norte de la Región de Coquimbo..
Este predicción muestra la misma tendencia de diferenciación en zonas calidad y frias según distribución geográfica en la región y también al comportamiento invernal. Al igual que con el mes de enero, las áreas de mayor incertumbre es donde hay menor densidad de estaciones. Según las métricas, para este mes el modelo sugiere que las predicciones son fiables.
Para este mes, la predicción muestra una transición estacional con un incremento de las temperaturas respecto a los meses invernales. La predicción de temperatura captura adecuadamente la variabilidad regional según el valor de R2, y además se aprecia un contraste térmico muy marcado entre los valles centrales con temperaturas media de 15°C y la cordillera principal con temperaturas bajo 0°C, lo que se relaciona fuertemente con las caracteristicas topográfica de dicho sector.
La predicción de temperaturas para diciembre muestra una transición hacia un clima más cálido, característico del inicio del verano, con una buena representación de la variación espacial de la temperatura en la región. La mayor precisión en las áreas con una mejor cobertura de datos (valles y zonas costeras) contrasta con la mayor incertidumbre en la cordillera. A pesar que este mes tuvo un R2 más bajo que los otros meses, el ajuste del modelo es suficiente para capturas las principales tendencias espaciales de temperatura para este periodo.
4. Discusión y conclusión
La comparación de la predicción de temperatura en la Región de Coquimbo utilizando los métodos del Vecino Más Cercano (KNN), Inverso a la Distancia (IDW) y el Modelamiento Lineal (Trend Surface Fitting - TSF) arroja diversas observaciones tanto a partir de las métricas como de los mapas generados.
El método KNN muestra una notable sensibilidad al número de vecinos considerados (k), lo cual influye en la suavización de la superficie interpolada a medida que aumenta el valor de k. Con valores bajos de k, como k=5, se generan mapas de temperatura con mayor variabilidad espacial y algunos artefactos en las áreas con menos datos. Sin embargo, este enfoque logra capturar con detalle las variaciones locales, especialmente en la zona de la Cordillera Principal. Para los cuatro meses analizados, un valor de k=10 mostró un mejor rendimiento en términos de explicación de la variabilidad y menor error (RMSE), en comparación con k=5. Sin embargo, al aumentar a k=15, no se observaron mejoras significativas, probablemente debido a que un mayor número de vecinos puede llevar a una sobre-generalización y a la pérdida de detalles locales.
Por su parte, el método IDW ofrece un desempeño intermedio entre la capacidad de suavizar la variabilidad local y la precisión de la predicción. Las métricas obtenidas para IDW, con valores de RMSE entre 2.48 y 3.56 y un R² que varía entre 0.23 y 0.36, reflejan una moderada capacidad para explicar la variabilidad en los datos de temperatura. Aunque el IDW logra una distribución espacial más uniforme que el KNN, sigue siendo sensible a la densidad de los puntos de medición, lo que genera variabilidad en las zonas con menos datos. En general, el IDW proporciona una interpolación aceptable, pero no logra capturar de manera eficaz los amplios gradientes de temperatura como lo hace un método de ajuste de tendencia más complejo.
El método TSF, que se basa en la creación de una superficie de tendencia ajustada mediante regresión, ha mostrado el mejor desempeño en términos de métricas, con valores de RMSE más bajos (entre 1.59 y 2.42) y R² más altos (llegando hasta 0.80). Esto indica una mayor capacidad para ajustar y capturar la variabilidad espacial de la temperatura, especialmente en meses como junio, cuando la estacionalidad fría se refleja claramente en los datos. El TSF es particularmente útil para modelar la variabilidad a gran escala y ofrece una representación más estable y continua de la temperatura en toda la región. Los mapas generados con este método presentan un ajuste visualmente más coherente y una menor variabilidad en áreas con poca densidad de datos.
Cada método de interpolación presenta ventajas y limitaciones específicas, dependiendo de la densidad de los datos y las características del área de estudio. KNN es útil para capturar detalles locales en áreas con alta densidad de estaciones, pero su desempeño es limitado cuando los datos son más dispersos. IDW ofrece un equilibrio adecuado, con una buena suavización de la variabilidad local, aunque su capacidad explicativa es más moderada. En cambio, el TSF destaca como el método más robusto en términos de precisión y ajuste, siendo especialmente eficaz para modelar la variabilidad a gran escala, lo que lo hace ideal para análisis regionales.
En conclusión, considerando los resultados obtenidos y el tipo de variable que se está prediciendo, así como la naturaleza de los datos disponibles, el TSF se presenta como la mejor opción de método de interpolación no geoestadístico. Esto se debe a su capacidad para representar con mayor precisión la variabilidad espacial, teniendo en cuenta las características morfoestructurales de la zona de estudio.
5. Referencias
Andrian, R., Utomo, W. H., & Handani, D. I. (2019). Implementation of the K-Nearest Neighbor (KNN) algorithm to recognize different motifs of Batik Lampung. Journal of Physics: Conference Series, 1338, 012061. https://doi.org/10.1088/1742-6596/1338/1/012061
Burrough, P., & McDonnell, R. (1998). Principles of geographical information systems. Oxford University Press.
Esri. (s.f.). Performing cross validation and validation. ArcGIS Pro. https://pro.arcgis.com/es/pro-app/latest/help/analysis/geostatistical-analyst/performing-cross-validation-and-validation.htm
