Introducción

Este informe analiza la predicción de temperaturas en una región de aguacates utilizando tres métodos de interpolación espacial: Voronoi, IDW (Inverse Distance Weighting) y un modelo de ensamble. Se realiza una validación cruzada de 5 particiones para evaluar la precisión de las predicciones, midiendo el RMSE en cada partición. Los resultados obtenidos se comparan para identificar el modelo más efectivo, y se visualizan mediante mapas y gráficos, mostrando cómo cada técnica captura la variabilidad espacial de la temperatura en el área de estudio.

Carga de Datos y Librerías

En esta sección se realiza la configuración inicial necesaria para el análisis de datos espaciales. Esto incluye la carga de las librerías requeridas y la importación de un conjunto de datos que contiene información de temperaturas asociadas a coordenadas geográficas específicas. Se filtran los datos para una fecha particular (01/10/2020) y se seleccionan las columnas relevantes: latitud, longitud y temperatura. Este paso es esencial para preparar los datos para los análisis espaciales e interpolaciones posteriores.

library(readxl)
library(geoR)
library(ggplot2)
library(sf)
library(sp)
library(gstat)
library(FNN)
library(terra)
library(dismo)
library(tmap)
library(dplyr)

ruta_excel <- "C:/Users/vicod/Documents/Maestria/analisis de datos espaciales y geograficos/Unidad 3/Datos_Completos_Aguacate.xlsx"
datos_aguacate <- read_excel(ruta_excel)
datos_aguacate <- datos_aguacate %>% 
                  filter(grepl("01/10/2020", FORMATTED_DATE_TIME))
datos_aguacate <- datos_aguacate %>% 
                  select(Latitude, Longitude, Temperature)

Exploración y Transformación de Datos Espaciales

Este código explora los datos y los transforma en un formato espacial. Se genera un resumen estadístico, un gráfico de dispersión de las coordenadas, y se convierte el conjunto de datos en un objeto espacial (sf) con un sistema métrico UTM para facilitar el análisis. Además, se visualizan las temperaturas en un mapa temático y se crea un marco de datos con coordenadas reproyectadas para análisis posteriores.

summary(datos_aguacate)
##     Latitude       Longitude       Temperature   
##  Min.   :2.392   Min.   :-76.71   Min.   :22.20  
##  1st Qu.:2.393   1st Qu.:-76.71   1st Qu.:24.50  
##  Median :2.393   Median :-76.71   Median :25.80  
##  Mean   :2.393   Mean   :-76.71   Mean   :25.83  
##  3rd Qu.:2.393   3rd Qu.:-76.71   3rd Qu.:27.18  
##  Max.   :2.394   Max.   :-76.71   Max.   :29.70
plot(
  datos_aguacate$Longitude,  
  datos_aguacate$Latitude,   
  col = "blue",              
  pch = 20,                  
  xlab = "Longitude",        
  ylab = "Latitude",         
)

datos_sf <- st_as_sf(
  datos_aguacate,
  coords = c("Longitude", "Latitude"),
  crs = 4326 #wgs84 
)


datos_sf <- st_transform(datos_sf, crs = 32614)
coords_metros <- st_coordinates(datos_sf)

ggplot(datos_sf) +
  geom_sf(aes(color = Temperature), size = 0.9) +  # Color según temperatura
  scale_color_gradient(low = "blue", high = "orange") +  # Gradiente de color
  theme_bw() +  # Tema limpio para el mapa
  labs(
    title = "Mapa de temperaturas",
    color = "Temperatura"
  )

Análisis de Estructura Espacial y Semivariograma

En esta sección, los datos de temperatura, con coordenadas reproyectadas, se convierten en un objeto de tipo geoR para realizar análisis geoestadísticos. Se examina la distribución espacial mediante un gráfico de dispersión y se calcula el semivariograma empírico, que cuantifica la dependencia espacial de la variable en función de la distancia entre puntos. Para refinar el análisis, se utiliza un semivariograma agrupado (binned), segmentando las distancias en intervalos definidos. Además, se genera un intervalo de confianza mediante simulaciones de Monte Carlo, lo que permite determinar si la estructura espacial observada es estadísticamente significativa. Los resultados indican la existencia de correlación espacial en los datos, ya que los puntos del semivariograma empírico se encuentran fuera del intervalo generado por las simulaciones, confirmando una relación espacial en la variable analizada.

datos_aguacate_metros <- data.frame(
  x = coords_metros[, 1],
  y = coords_metros[, 2],
  Temperature = datos_aguacate$Temperature
)

datos_g <- as.geodata(datos_aguacate_metros, coords.col = 1:2, data.col = 3)
plot(datos_g)

summary(dist(datos_g$coords))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.046  48.562  76.796  81.797 109.955 234.895
plot(variog(coords = datos_g$coords, data = datos_g$data,
            option = "cloud", max.dist = 115))
## variog: computing omnidirectional variogram

variograma <- variog(datos_g, option = "bin", uvec = seq( 0, 115, 5))
## variog: computing omnidirectional variogram
plot(variograma)

random=variog.mc.env(datos_g,obj.variog = variograma)
## variog.env: generating 99 simulations by permutating data values
## variog.env: computing the empirical variogram for the 99 simulations
## variog.env: computing the envelops
plot(variograma,envelope = random )

Ajuste del Modelo de Semivariograma

Se ajustan tres modelos teóricos (exponencial, gaussiano y esférico) al semivariograma empírico mediante mínimos cuadrados ordinarios. Los modelos se visualizan en el gráfico del semivariograma con curvas de diferentes colores: rojo para el modelo exponencial, azul para el gaussiano y verde para el esférico. Al evaluar los resultados mediante la suma de cuadrados de los errores, el modelo exponencial (rojo) es el que mejor se ajusta a los datos, ya que presenta la menor suma de cuadrados, lo que lo convierte en el modelo seleccionado para representar la estructura espacial de los datos.

plot(variograma)
ini.vals = expand.grid(seq(2,4,l=0.3),seq(20,80,l =3))
m_mco_exp = variofit(variograma, ini = ini.vals, cov.model ="exp")
## variofit: covariance model used is exponential 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi   tausq kappa
## initial.value "2"     "20"  "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 166355.481792341
m_mco_gaus = variofit(variograma, ini = ini.vals, cov.model ="gaus")
## variofit: covariance model used is gaussian 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi   tausq kappa
## initial.value "2"     "20"  "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 148966.844805868
m_mco_spe = variofit(variograma, ini = ini.vals, cov.model ="sph")
## variofit: covariance model used is spherical 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi   tausq kappa
## initial.value "2"     "20"  "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 130279.742164586
lines(m_mco_exp,col ="red")
lines(m_mco_gaus,col ="blue")
lines(m_mco_spe,col ="green")

summary(m_mco_exp)$sum.of.squares
##    value 
## 2309.283
summary(m_mco_gaus)$sum.of.squares
##    value 
## 3656.729
summary(m_mco_spe)$sum.of.squares
##    value 
## 3089.164
summary(m_mco_exp)
## $pmethod
## [1] "WLS (weighted least squares)"
## 
## $cov.model
## [1] "exponential"
## 
## $spatial.component
##   sigmasq       phi 
##  1.958089 35.419481 
## 
## $spatial.component.extra
## kappa 
##   0.5 
## 
## $nugget.component
##    tausq 
## 1.478306 
## 
## $fix.nugget
## [1] FALSE
## 
## $fix.kappa
## [1] TRUE
## 
## $practicalRange
## [1] 106.1073
## 
## $sum.of.squares
##    value 
## 2309.283 
## 
## $estimated.pars
##     tausq   sigmasq       phi 
##  1.478306  1.958089 35.419481 
## 
## $weights
## [1] "npairs"
## 
## $call
## variofit(vario = variograma, ini.cov.pars = ini.vals, cov.model = "exp")
## 
## attr(,"class")
## [1] "summary.variomodel"

Generación de la Grilla de Predicción y Kriging

Se crea una grilla de predicción de 100x100 puntos en el área de estudio, basada en el rango de las coordenadas de las muestras. Luego, se ajusta el modelo exponencial al semivariograma y se aplica el método de Kriging para predecir los valores de temperatura en cada punto de la grilla. El resultado se visualiza en un mapa que muestra la distribución espacial de la temperatura interpolada, con contornos que indican áreas de igual valor.

x_range <- range(coords_metros[, 1])
y_range <- range(coords_metros[, 2])

datos_grid <- expand.grid(
  Este = seq(x_range[1], x_range[2], length.out = 100),
  Norte = seq(y_range[1], y_range[2], length.out = 100)
)

m_exp <- variofit(variograma, ini = ini.vals, cov.model = "exp")
## variofit: covariance model used is exponential 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi   tausq kappa
## initial.value "2"     "20"  "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 166355.481792341
kriging_result <- krige.conv(
  datos_g, 
  loc = datos_grid, 
  krige = krige.control(
    trend.d = "cte", 
    trend.l = "cte",
    cov.model = "exponential",
    cov.pars = m_exp$cov.pars,
    nugget = m_exp$nugget
  )
)
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
image(
  kriging_result, 
  main = "Prediccin por Kriging", 
  xlab = "Est", 
  ylab = "Nort", 
  col = terrain.colors(100)
)
contour(
  kriging_result, add = TRUE, col = "black"
)

Creación de Celdas de Voronoi y Visualización de Temperaturas

En este fragmento de código, se crea un objeto espacial vect a partir de las coordenadas de los datos de temperatura. Luego, se define el área de estudio y se generan las celdas de Voronoi utilizando las coordenadas como puntos de referencia. Las celdas de Voronoi dividen el espacio en regiones, donde cada celda corresponde a la zona más cercana a un punto de muestreo. A continuación, se extraen las temperaturas correspondientes a cada celda y se visualizan, con las celdas de Voronoi coloreadas según la temperatura. Además, se superponen los puntos originales sobre el mapa para mejor referencia visual.

d <- data.frame(
  x = datos_aguacate_metros$x,
  y = datos_aguacate_metros$y,
  temperature = datos_aguacate_metros$Temperature
)

d_vect <- terra::vect(d, c("x", "y"))


bnd <- terra::ext(min(d$x), max(d$x), min(d$y), max(d$y))


v <- terra::voronoi(d_vect, bnd = bnd)


v$temperature <- terra::extract(v, d_vect)$temperature


plot(v, col = terrain.colors(100)[as.numeric(cut(v$temperature, breaks = 100))])
points(d_vect, cex = 0.5, col = "black")

Interpolación Espacial mediante IDW (Invers Distance Weighting)

En este fragmento de código, se realiza una interpolación espacial utilizando el método de IDW (Invers Distance Weighting). Primero, se crea un modelo IDW usando el paquete gstat, con la temperatura como variable dependiente y las coordenadas como variables explicativas. Luego, se genera una cuadrícula de predicción y se aplica el modelo IDW para predecir la temperatura en esos puntos. Los resultados se convierten en un SpatRaster para visualización, y se utiliza tmap para mostrar la predicción en un mapa, donde las temperaturas estimadas se visualizan con una paleta de colores.

idw_model <- gstat(
  formula = temperature ~ 1, 
  locations = ~ x + y,       
  data = d,                  
  nmax = 6,            
  set = list(idp = 3)        
)

x_range <- seq(min(d$x), max(d$x), length.out = 100) 
y_range <- seq(min(d$y), max(d$y), length.out = 100)

grid <- expand.grid(x = x_range, y = y_range)

coordinates(grid) <- ~ x + y
gridded(grid) <- TRUE


pred_idw <- predict(idw_model, newdata = grid)
## [inverse distance weighted interpolation]
resp <- as.data.frame(pred_idw)
resp$x <- coordinates(grid)[, 1]
resp$y <- coordinates(grid)[, 2]
resp$pred <- resp$var1.pred

pred_raster <- rast(resp, type = "xyz", crs = st_crs(datos_sf)$proj4string)

terra::nlyr(pred_raster)
## [1] 3
pred_raster_single <- pred_raster[["pred"]]
tm_shape(pred_raster_single) +
  tm_raster(palette = "viridis", alpha = 0.6, title = "Prediccion (IDW)") +
  tm_layout(main.title = "Interpolacion por IDW", legend.outside = TRUE)

Interpolación Espacial por Vecinos más Cercanos (IDW con idp = 0)

En este fragmento de código, se utiliza el método de interpolación IDW (Invers Distance Weighting) con el parámetro idp = 0, lo que indica que la predicción de la temperatura en cada punto se realiza utilizando solo los vecinos más cercanos, sin ponderación por la distancia. Se define un modelo IDW con nmax = 4, limitando el número de vecinos a 4. A continuación, se realiza la predicción sobre una cuadrícula de puntos generada previamente, y los resultados se visualizan en un mapa usando tmap, donde las temperaturas estimadas se representan mediante una paleta de colores.

vecinos <- gstat(
  formula = temperature ~ 1, 
  locations = ~ x + y,       
  data = d,           
  nmax = 4,           
  set = list(idp = 0) 
)

pred_vecinos <- predict(vecinos, newdata = grid)
## [inverse distance weighted interpolation]
resp <- as.data.frame(pred_vecinos)
resp$x <- coordinates(grid)[, 1]
resp$y <- coordinates(grid)[, 2]
resp$pred <- resp$var1.pred

pred_raster <- rast(resp, type = "xyz", crs = st_crs(datos_sf)$proj4string)

terra::nlyr(pred_raster)
## [1] 3
pred_raster_single <- pred_raster[["var1.pred"]]
tm_shape(pred_raster_single) +
  tm_raster(palette = "viridis", alpha = 0.6, title = "Prediccin (vecinos)") +
  tm_layout(main.title = "Interpolacin por vecinos", legend.outside = TRUE)

Validación Cruzada y Comparación de Métodos de Interpolación

En este código se realiza una validación cruzada 5-fold para comparar el rendimiento de diferentes métodos de interpolación (Voronoi, IDW con idp = 3, IDW con idp = 0, y un ensamble de estos). En cada iteración, se calcula el RMSE para evaluar la precisión de las predicciones, comparando las temperaturas estimadas con las observadas. Además, para el ensamble, se combinan las predicciones de los tres métodos, ponderadas por su RMSE, para obtener una estimación final. Esto permite identificar el modelo que ofrece la mejor predicción de la temperatura.

RMSE <- function(observed, predicted) {
  sqrt(mean((observed - predicted)^2))
}

set.seed(123)

kf <- dismo::kfold(nrow(d), k = 5)

rmse_voronoi <- rep(NA, 5) 
rmse_idw_1 <- rep(NA, 5)   
rmse_idw_0 <- rep(NA, 5)  
rmse_ensemble <- rep(NA, 5) 




for (k in 1:5) {
  
  test <- d[kf == k, ]
  train <- d[kf != k, ]
  
  train_vect <- terra::vect(train, c("x", "y"))
  
  v <- terra::voronoi(train_vect, bnd = terra::ext(min(d$x), max(d$x), min(d$y), max(d$y)))
  
  v$temperature <- terra::extract(v, terra::vect(test, c("x", "y"))) 
  p1 <- v$temperature
  

  rmse_voronoi[k] <- RMSE(test$temperature, p1)
  

  idw_model_1 <- gstat(
    formula = temperature ~ 1,
    locations = ~ x + y,
    data = train,
    nmax = 6,
    set = list(idp = 3)
  )
  p2 <- predict(idw_model_1, newdata = test)$var1.pred
  rmse_idw_1[k] <- RMSE(test$temperature, p2)
  

  idw_model_0 <- gstat(
    formula = temperature ~ 1,
    locations = ~ x + y,
    data = train,
    nmax = 6,
    set = list(idp = 0)
  )
  p3 <- predict(idw_model_0, newdata = test)$var1.pred
  rmse_idw_0[k] <- RMSE(test$temperature, p3)
  
  #  Ensamble 
  w <- 1 / c(rmse_voronoi[k], rmse_idw_1[k], rmse_idw_0[k])
  weights <- w / sum(w)
  
  p4 <- p1 * weights[1] + p2 * weights[2] + p3 * weights[3]
  rmse_ensemble[k] <- RMSE(test$temperature, p4)
}

Cálculo del Promedio de RMSE

En esta parte del código, se calcula el promedio del RMSE sobre las 5 particiones de la validación cruzada para cada uno de los métodos de interpolación (Voronoi, IDW con idp = 1, IDW con idp = 0, y el modelo de ensamble). Estos valores proporcionan una evaluación general del desempeño de cada método. Los resultados se presentan en un data frame, mostrando la media del RMSE para cada técnica y permitiendo comparar su efectividad en la predicción de las temperaturas.

# Resultados del RMSE obtenidos para cada uno de los 5 splits
resultados_rmse <- data.frame(
  Voronoi = rmse_voronoi,
  IDW_idp3 = rmse_idw_1,
  IDW_idp0 = rmse_idw_0,
  Ensamble = rmse_ensemble
)

# Promedio de RMSE sobre las 5 particiones
promedio_rmse <- data.frame(
  Voronoi = mean(rmse_voronoi),
  IDW_idp3 = mean(rmse_idw_1),
  IDW_idp0 = mean(rmse_idw_0),
  Ensamble = mean(rmse_ensemble)
)

# Mostrar los resultados
print(resultados_rmse)
##    Voronoi IDW_idp3 IDW_idp0 Ensamble
## 1 41.19481 1.091162 1.136826 1.316229
## 2 41.52005 1.069569 1.127216 1.202038
## 3 40.48358 1.126595 1.115893 1.143574
## 4 41.24370 1.101768 1.056408 1.143588
## 5 41.27278 1.259519 1.304061 1.415685
print(promedio_rmse)
##    Voronoi IDW_idp3 IDW_idp0 Ensamble
## 1 41.14298 1.129723 1.148081 1.244223

Conclusión sobre el Análisis de RMSE

Los resultados del análisis de RMSE sobre las 5 particiones de la validación cruzada muestran cómo se desempeñan los distintos métodos de interpolación. Los valores promedio de RMSE para cada técnica son los siguientes:

  • Voronoi: 41.14
  • IDW con idp = 3: 1.13
  • IDW con idp = 0: 1.15
  • Ensamble: 1.24

De estos resultados, podemos concluir que el método Voronoi presenta un RMSE significativamente más alto, lo que indica que su desempeño en la predicción de temperaturas es inferior en comparación con los métodos IDW y el modelo Ensamble.

Entre los métodos IDW, el que utiliza idp = 3 (distancia ponderada inversamente) tiene un rendimiento ligeramente mejor (RMSE = 1.13) que el que utiliza idp = 0 (vecinos más cercanos), aunque la diferencia es pequeña.

El modelo de ensamble, que combina los tres métodos, presenta un RMSE de 1.24, lo que sugiere que no mejora significativamente la precisión sobre los modelos individuales en este caso.

En resumen, los métodos IDW parecen ser más precisos que el Voronoi, y el uso del Ensamble no proporciona una mejora clara en comparación con la selección de un solo modelo como el IDW con idp = 3, esto puede deberse al bajo desempeño del modelo voronoi que afecta significativamente el ensamble.