Repoorte De Aguacates Actividad 3
Carlos Camilo Cuchumbe Escandon
2024-12-09
Librerias Utilizadas
library(readxl)
library(geoR)
library(ggplot2)
library(sf)
library(sp)
library(gstat)
library(FNN)
library(terra)
library(dismo)
library(tmap)
library(dplyr)
library(readxl)Lectura de los datos
datos_aguacate <- read_excel("C:/Users/camil/Downloads/Datos_Completos_Aguacate.xlsx")
datos_aguacate <- datos_aguacate %>% filter(grepl("01/10/2020", FORMATTED_DATE_TIME))
head(datos_aguacate)## # A tibble: 6 × 21
## id_arbol Latitude Longitude FORMATTED_DATE_TIME Psychro_Wet_Bulb_Temp…¹
## <chr> <dbl> <dbl> <chr> <dbl>
## 1 1 2.39 -76.7 01/10/2020 10:11:12 a, m, 22
## 2 2 2.39 -76.7 01/10/2020 10:11:12 a, m, 21.4
## 3 3 2.39 -76.7 01/10/2020 10:11:12 a, m, 21.8
## 4 4 2.39 -76.7 01/10/2020 10:11:12 a, m, 22.8
## 5 5 2.39 -76.7 01/10/2020 10:11:12 a, m, 22.6
## 6 6 2.39 -76.7 01/10/2020 10:11:12 a, m, 21.5
## # ℹ abbreviated name: ¹Psychro_Wet_Bulb_Temperature
## # ℹ 16 more variables: Station_Pressure <dbl>, Relative_Humidity <dbl>,
## # Crosswind <dbl>, Temperature <dbl>, Barometric_Pressure <dbl>,
## # Headwind <dbl>, Direction_True <dbl>, Direction_Mag <dbl>,
## # Wind_Speed <dbl>, Heat_Stress_Index <dbl>, Altitude <dbl>, Dew_Point <dbl>,
## # Density_Altitude <dbl>, Wind_Chill <dbl>,
## # Estado_Fenologico_Predominante <dbl>, Frutos_Afectados <dbl>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.
Resumen de los datos
datos_aguacate <- datos_aguacate %>% select(Latitude, Longitude, Temperature)
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.70Este 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.
Grafico De Cultivos
plot(
datos_aguacate$Longitude, # Eje X: longitud
datos_aguacate$Latitude, # Eje Y: latitud
col = "blue", # Color de los puntos
pch = 20, # Tipo de marcador (punto sólido)
xlab = "Longitude", # Etiqueta para el eje X
ylab = "Latitude", # Etiqueta para el eje Y
)
Grafico por medio de Temperaturas
datos_sf <- sf::st_as_sf( datos_aguacate,
coords = c("Longitude", "Latitude"),
crs = 4326 #wgs84
)
datos_sf <- sf::st_transform(datos_sf, crs = 32614)
coords_metros <- sf::st_coordinates(datos_sf)
ggplot2::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.
Grafico Geo referenciado
datos_aguacate_metros <- data.frame(
x = coords_metros[, 1],
y = coords_metros[, 2],
Temperature = datos_aguacate$Temperature
)
# geoR
datos_g <- as.geodata(datos_aguacate_metros, coords.col = 1:2, data.col = 3)
plot(datos_g)
plot(variog(coords = datos_g$coords, data = datos_g$data,
option = "cloud", max.dist = 115))## variog: computing omnidirectional variogram
Resumen de los datos
summary(dist(datos_g$coords))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.046 48.562 76.796 81.797 109.955 234.895Semivariograma con distancias corregidas
variograma <- variog(datos_g, option = "bin", uvec = seq( 0, 115, 5))## variog: computing omnidirectional variogramplot(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 envelopsplot(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.
Grafico con las distribuciones Ajustadas
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.481792341m_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.844805868m_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.742164586plot(variograma)
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.283summary(m_mco_gaus)$sum.of.squares## value
## 3656.729summary(m_mco_spe)$sum.of.squares## value
## 3089.164summary(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.
Generar grilla de predicción
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.481792341Kriging con el modelo ajustado
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 neighbourhoodVisualización del resultado kriging_result
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.
# Crear un objeto de tipo vect con las coordenadas x, y (suponiendo que tienes 'x' y 'y' de tus datos)
d <- data.frame(
x = datos_aguacate_metros$x,
y = datos_aguacate_metros$y,
temperature = datos_aguacate_metros$Temperature
)
# Convertir el data frame a un objeto de tipo vect (terreno espacial)
d_vect <- terra::vect(d, c("x", "y"))
# Definir los límites de la región para la creación de las celdas de Voronoi
bnd <- terra::ext(min(d$x), max(d$x), min(d$y), max(d$y))
# Crear las celdas de Voronoi
v <- terra::voronoi(d_vect, bnd = bnd)
# Extraer la temperatura más cercana para cada celda de Voronoi
v$temperature <- terra::extract(v, d_vect)$temperature
vector_1=v$temperature
# Visualizar las celdas de Voronoi coloreadas por la temperatura
plot(v, col = terrain.colors(100)[as.numeric(cut(vector_1, breaks = 100))])
points(d_vect, cex = 0.5, col = "black") # Agregar los puntos originales
### 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.
Crear el objeto ‘gstat’ para IDW
idw_model <- gstat::gstat(
formula = temperature ~ 1, # Fórmula simple, predice solo usando el valor de la temperatura
locations = ~ x + y, # Coordenadas
data = d, # Datos
nmax = 6, # Usar todos los vecinos
set = list(idp = 3) # IDW con beta = 1 (ponderación inversa a la distancia)
)
# Crear una cuadrícula de predicción (raster)
x_range <- seq(min(d$x), max(d$x), length.out = 100) # Ajustar la resolución según sea necesario
y_range <- seq(min(d$y), max(d$y), length.out = 100)
# Crear un data frame de la grilla
grid <- expand.grid(x = x_range, y = y_range)
head(grid)## x y
## 1 3043166 285884.3
## 2 3043168 285884.3
## 3 3043170 285884.3
## 4 3043172 285884.3
## 5 3043174 285884.3
## 6 3043176 285884.3Metodo por Interpolacin por IDW
# Convertir la grilla a Spatial (requerido para predicción con gstat)
coordinates(grid) <- ~ x + y
gridded(grid) <- TRUE
# Realizar la predicción IDW en la cuadrícula
pred_idw <- predict(idw_model, newdata = grid)## [inverse distance weighted interpolation]# Convertir los resultados de predicción en un data frame para posterior visualización
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] 3pred_raster_single <- pred_raster[["var1.pred"]]
tm_shape(pred_raster_single) +
tm_raster(palette = "viridis", alpha = 0.6, title = "Prediccion (IDW)") +
tm_layout(main.title = "Interpolacin 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.
Metodo por Interpolacin por vecinos
# Crear el objeto 'gstat' para IDW
vecinos <- gstat(
formula = temperature ~ 1, # Fórmula simple, predice solo usando el valor de la temperatura
locations = ~ x + y, # Coordenadas
data = d, # Datos
nmax = 4, # Usar todos los vecinos
set = list(idp = 0) # IDW con beta = 1 (ponderación inversa a la distancia)
)
# Realizar la predicción IDW en la cuadrícula
pred_vecinos <- predict(vecinos, newdata = grid)## [inverse distance weighted interpolation]# Convertir los resultados de predicción en un data frame para posterior visualización
resp <- as.data.frame(pred_vecinos)
resp$x <- coordinates(grid)[, 1]
resp$y <- coordinates(grid)[, 2]
resp$pred <- resp$var1.pred
# Convertir a SpatRaster para visualización con tmap
pred_raster <- rast(resp, type = "xyz", crs = st_crs(datos_sf)$proj4string)
terra::nlyr(pred_raster)## [1] 3pred_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.
Función para calcular el RMSE
RMSE <- function(observed, predicted) {
sqrt(mean((observed - predicted)^2))
}Funcion para el calculo de metricas
# Establecer la semilla para reproducibilidad
set.seed(123)
# Dividir los datos en 5 subconjuntos (validación cruzada 5-fold)
kf <- dismo::kfold(nrow(d), k = 5)
# Vectores para almacenar los RMSE obtenidos con cada método
rmse_voronoi <- rep(NA, 5) # Voronoi
rmse_idw_1 <- rep(NA, 5) # IDW con idp = 1
rmse_idw_0 <- rep(NA, 5) # IDW con idp = 0
rmse_ensemble <- rep(NA, 5) # Ensamble
for (k in 1:5) {
# Dividir los datos en conjuntos de entrenamiento y prueba
test <- d[kf == k, ]
train <- d[kf != k, ]
# *** Voronoi ***
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"))) # Asegúrate de que test también esté en formato espacial
p1 <- v$temperature
rmse_voronoi[k] <- RMSE(test$temperature, p1)
# *** IDW con idp = 1 ***
idw_model_1 <- gstat::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 con idp = 0 ***
idw_model_0 <-gstat:: 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 de métodos ***
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)
}## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]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 Error Cuadratico Medio
resultados_rmse <- data.frame(
Voronoi = rmse_voronoi,
IDW_idp1 = rmse_idw_1,
IDW_idp0 = rmse_idw_0,
Ensamble = rmse_ensemble
)
resultados_rmse## Voronoi IDW_idp1 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.415685Promedio de los resultados obtenidos
promedio_rmse <- data.frame(
Voronoi = mean(rmse_voronoi),
IDW_idp1 = mean(rmse_idw_1),
IDW_idp0 = mean(rmse_idw_0),
Ensamble = mean(rmse_ensemble)
)
promedio_rmse## Voronoi IDW_idp1 IDW_idp0 Ensamble
## 1 41.14298 1.129723 1.148081 1.244223Conclusió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.