Instrucciones:

Se empleará el archivo “Datos completos aguacate.xls” para extraer un subconjunto de datos basado en la columna “FORMATTED_DATA_TIME”, donde se seleccionará la fecha “01/10/2020”, lo que resultará en 534 filas de datos. A continuación, se llevarán a cabo las siguientes operaciones centradas en la variable de temperatura:

Autocorrelación espacial: Se calculará el semivariograma.
Identificación del mejor modelo teórico: Se probarán los modelos exponencial, gaussiano y esférico.
Predicción espacial mediante el método Kriging: Se realizará una interpolación.
Generación de una imagen de interpolación o mapa.

En caso de que la temperatura carezca de una autocorrelación espacial significativa, se explorará el uso de otras cuatro variables.

Preparación de los datos:

# Cargar la librería readxl para leer archivos de Excel
library(readxl)

# Suprimir mensajes y advertencias durante la lectura del archivo Excel
suppressMessages({
  suppressWarnings({
    BD_aguacate <- read_excel("Datos_Completos_Aguacate.xlsx")
  })
})

# Filtrar el conjunto de datos para mantener solo las filas con la fecha "01/10/2020"
BD_aguacate_f1 <- subset(BD_aguacate, grepl("01/10/2020", FORMATTED_DATE_TIME), 
                         select = c(id_arbol, Latitude, Longitude,
                                    FORMATTED_DATE_TIME, Relative_Humidity,
                                    Temperature, Altitude, Wind_Speed))

# Mostrar la estructura del nuevo conjunto de datos filtrado
str(BD_aguacate_f1)
## tibble [534 × 8] (S3: tbl_df/tbl/data.frame)
##  $ id_arbol           : chr [1:534] "1" "2" "3" "4" ...
##  $ Latitude           : num [1:534] 2.39 2.39 2.39 2.39 2.39 ...
##  $ Longitude          : num [1:534] -76.7 -76.7 -76.7 -76.7 -76.7 ...
##  $ FORMATTED_DATE_TIME: chr [1:534] "01/10/2020  10:11:12 a, m," "01/10/2020  10:11:12 a, m," "01/10/2020  10:11:12 a, m," "01/10/2020  10:11:12 a, m," ...
##  $ Relative_Humidity  : num [1:534] 85.2 84 79.6 77.6 76.5 77.7 76.5 77.7 78.3 77.8 ...
##  $ Temperature        : num [1:534] 23.9 23.5 24.5 25.9 26 24.5 25.5 25.7 26 25.9 ...
##  $ Altitude           : num [1:534] 1696 1696 1694 1694 1696 ...
##  $ Wind_Speed         : num [1:534] 0 0 0.5 0.5 0 0.3 0 0 0 0.4 ...

Se seleccionaron 4 variables de interés, con el objetivo de comparar la correlación espacial existente en cada uno y poder evidencar las diferencias en la interpolación.

# Cargar la librería geoR para análisis geoespaciales
library(geoR)
## --------------------------------------------------------------
##  Analysis of Geostatistical Data
##  For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
##  geoR version 1.9-4 (built on 2024-02-14) is now loaded
## --------------------------------------------------------------
# Convertir las coordenadas de latitud a segundos
BD_aguacate_f1$Latitude <-  BD_aguacate_f1$Latitude * 3600

# Convertir las coordenadas de longitud a segundos
BD_aguacate_f1$Longitude <-  BD_aguacate_f1$Longitude * 3600

# Crear objetos de datos geoespaciales para la humedad relativa, temperatura, altitud y velocidad del viento
geo_humid <- as.geodata(BD_aguacate_f1, coords.col = c(3,2), data.col = 5)
geo_temp <- as.geodata(BD_aguacate_f1, coords.col = c(3,2), data.col = 6)
geo_altit <- as.geodata(BD_aguacate_f1, coords.col = c(3,2), data.col = 7)
geo_ws <- as.geodata(BD_aguacate_f1, coords.col = c(3,2), data.col = 8)
# Cargar la librería geoR para análisis geoespaciales
library(geoR)

# Graficar los datos de humedad relativa
plot(geo_humid)
# Agregar un título a la gráfica de humedad relativa
title(main = "Geodatos de Humedad", cex.main = 0.9, adj = 0.45, line = -20)

# Graficar los datos de temperatura
plot(geo_temp)
# Agregar un título a la gráfica de temperatura
title(main = "Geodatos de Temperatura", cex.main = 0.9, adj = 0.45, line = -20)

# Graficar los datos de altitud
plot(geo_altit)
# Agregar un título a la gráfica de altitud
title(main = "Geodatos de Altitud", cex.main = 0.9, adj = 0.45, line = -20)

# Graficar los datos de velocidad del viento
plot(geo_ws)
# Agregar un título a la gráfica de velocidad del viento
title(main = "Geodatos de Velocidad del Viento", cex.main = 0.9, adj = 0.45, line = -20)

Los geodatos mostrados en las gráficas anteriores muestran la dispersión de los datos y el histograma para cada variable seleccionada, es de resaltar que en las tres primeras variables, Humedad, Temperatura y Altitud, la dispersión no refleja patrones y los histogramas tienen a ser normales, sin embargo, la última variables, velocidad del viento, muestra patrones y una distribución diferente.

Autocorrelación espacial:

# Cargar la librería geoR para análisis geoespaciales
library(geoR)

# Imprimir un mensaje sobre el cálculo de la distancia entre puntos
print("Distancia entre puntos")
## [1] "Distancia entre puntos"
# Calcular la matriz de distancias entre los puntos geográficos
md <- dist(geo_humid$coords)

# Resumir estadísticas sobre las distancias calculadas
summary(md)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.06162 1.45844 2.30681 2.45763 3.30394 7.05286

Debido las coordenadas de la zona del interés estaban dadas en Grados, la matriz de distancias tenía valores muy pequeños difíciles de manejar, dado lo anterior desde el inicio se convirtieron las coordenadas de grados a segundo para lograr un manejo más sencillo.

Variograma de Humedad Relativa

# Calcular el semivariograma para los datos de humedad
semivariograma_humid <- variog(geodata = geo_humid, uvec = seq(0, 4, 0.07))
## variog: computing omnidirectional variogram
# Calcular el envelope para el semivariograma mediante bootstrap
datos.env <- variog.mc.env(geo_humid, obj = semivariograma_humid)
## variog.env: generating 99 simulations by permutating data values
## variog.env: computing the empirical variogram for the 99 simulations
## variog.env: computing the envelops
# Graficar el semivariograma junto con su envelope
plot(semivariograma_humid, pch = 16, envelope = datos.env)

Variograma de Temperatura

# Calcular el semivariograma para los datos de temperatura
semivariograma_temp <- variog(geodata = geo_temp, uvec = seq(0, 4, 0.07))
## variog: computing omnidirectional variogram
# Calcular el envelope para el semivariograma mediante bootstrap
datos.env <- variog.mc.env(geo_temp, obj = semivariograma_temp)
## variog.env: generating 99 simulations by permutating data values
## variog.env: computing the empirical variogram for the 99 simulations
## variog.env: computing the envelops
# Graficar el semivariograma junto con su envelope
plot(semivariograma_temp, pch = 16, envelope = datos.env)

Variograma de Altitud

# Calcular el semivariograma para los datos de altitud
semivariograma_altit <- variog(geodata = geo_altit, uvec = seq(0, 4, 0.07))
## variog: computing omnidirectional variogram
# Calcular el envelope para el semivariograma mediante bootstrap
datos.env <- variog.mc.env(geo_altit, obj = semivariograma_altit)
## variog.env: generating 99 simulations by permutating data values
## variog.env: computing the empirical variogram for the 99 simulations
## variog.env: computing the envelops
# Graficar el semivariograma junto con su envelope
plot(semivariograma_altit, pch = 16, envelope = datos.env)

Variograma de Velocidad del Viento

# Calcular el semivariograma para los datos de velocidad del viento
semivariograma_ws <- variog(geodata = geo_ws, uvec = seq(0, 4, 0.07))
## variog: computing omnidirectional variogram
# Calcular el envelope para el semivariograma mediante bootstrap
datos.env <- variog.mc.env(geo_ws, obj = semivariograma_ws)
## variog.env: generating 99 simulations by permutating data values
## variog.env: computing the empirical variogram for the 99 simulations
## variog.env: computing the envelops
# Graficar el semivariograma junto con su envelope
plot(semivariograma_ws, pch = 16, envelope = datos.env)

De los semivarigramas anteriores, destaca que las variables Humedad Relativa y Altitud, tiene un rango de semivarianza amplio frente a las variables temperatura y velocidad del viento.

Por otro lado, la cantidad de puntos y cercanía a la zona de discriminación es muy evidente para la variable Velocidad del viendo, indicando así que no existe correlación espacial.

La variable altitud, presenta la mejor correlación espacial, seguida por la humedad relativa y luego la temperatura, las cuales se analizan más a detalle a continuación:

Modelos teóricos

# Definir una cuadrícula de valores iniciales para el ajuste del modelo variográfico
ini.vals1 <- expand.grid(seq(30, 50, length.out = 10), seq(3, 10, length.out = 10))

# Ajustar el modelo variográfico exponencial al semivariograma experimental de humedad relativa
model_mco_exp1 <- variofit(semivariograma_humid, ini = ini.vals1, cov.model = "exponential")
## 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 "50"    "3"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 3756374.08513128
# Ajustar el modelo variográfico gaussiano al semivariograma experimental de humedad relativa
model_mco_gaus1 <- variofit(semivariograma_humid, ini = ini.vals1, cov.model = "gauss")
## 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 "50"    "3"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 14981142.6498284
# Ajustar el modelo variográfico esférico al semivariograma experimental de humedad relativa
model_mco_spe1 <- variofit(semivariograma_humid, ini = ini.vals1, cov.model = "spherical")
## 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 "34.44" "3"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 2693907.4552161
# Graficar el semivariograma experimental y los modelos ajustados
plot(semivariograma_humid, pch = 16, col = "black", main = "Semivariograma Experimental y Modelos Ajustados Humedad R")
lines(model_mco_exp1, col = "blue")
lines(model_mco_gaus1, col = "red")
lines(model_mco_spe1, col = "purple")
legend("topright", legend = c("Exponencial", "Gaussiano", "Esférico"), col = c("blue", "red", "purple"), lty = 1, cex = 0.8)

# Definir una cuadrícula de valores iniciales para el ajuste del modelo variográfico
ini.vals2 <- expand.grid(seq(3, 6, length.out = 10), seq(3, 10, length.out = 10))

# Ajustar el modelo variográfico exponencial al semivariograma experimental de temperatura
model_mco_exp2 <- variofit(semivariograma_temp, ini = ini.vals2, cov.model = "exponential")
## 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 "6"     "3"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 57437.1014324913
# Ajustar el modelo variográfico gaussiano al semivariograma experimental de temperatura
model_mco_gaus2 <- variofit(semivariograma_temp, ini = ini.vals2, cov.model = "gauss")
## 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 "6"     "3"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 214686.967125714
# Ajustar el modelo variográfico esférico al semivariograma experimental de temperatura
model_mco_spe2 <- variofit(semivariograma_temp, ini = ini.vals2, cov.model = "spherical")
## 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 "3.67"  "3"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 44820.5776531328
# Graficar el semivariograma experimental y los modelos ajustados
plot(semivariograma_temp, pch = 16, col = "black", main = "Semivariograma Experimental y Modelos Ajustados Temperatura")
lines(model_mco_exp2, col = "blue")
lines(model_mco_gaus2, col = "red")
lines(model_mco_spe2, col = "purple")
legend("topright", legend = c("Exponencial", "Gaussiano", "Esférico"), col = c("blue", "red", "purple"), lty = 1, cex = 0.8)

# Definir una cuadrícula de valores iniciales para el ajuste del modelo variográfico
ini.vals3 <- expand.grid(seq(100, 200, length.out = 10), seq(3, 10, length.out = 10))

# Ajustar el modelo variográfico exponencial al semivariograma experimental de altitud
model_mco_exp3 <- variofit(semivariograma_altit, ini = ini.vals3, cov.model = "exponential")
## 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 "200"   "3"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 14960733.6484513
# Ajustar el modelo variográfico gaussiano al semivariograma experimental de altitud
model_mco_gaus3 <- variofit(semivariograma_altit, ini = ini.vals3, cov.model = "gauss")
## 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 "200"   "3"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 119057914.072068
# Ajustar el modelo variográfico esférico al semivariograma experimental de altitud
model_mco_spe3 <- variofit(semivariograma_altit, ini = ini.vals3, cov.model = "spherical")
## 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 "144.44" "3.78" "0"   "0.5"
## status        "est"    "est"  "est" "fix"
## loss value: 21901294.3185684
# Graficar el semivariograma experimental y los modelos ajustados
plot(semivariograma_altit, pch = 16, col = "black", main = "Semivariograma Experimental y Modelos Ajustados Altitud")
lines(model_mco_exp3, col = "blue")
lines(model_mco_gaus3, col = "red")
lines(model_mco_spe3, col = "purple")
legend("topright", legend = c("Exponencial", "Gaussiano", "Esférico"), col = c("blue", "red", "purple"), lty = 1, cex = 0.8)

# Definir una cuadrícula de valores iniciales para el ajuste del modelo variográfico
ini.vals4 <- expand.grid(seq(0, 10, length.out = 10), seq(3, 10, length.out = 10))

# Ajustar el modelo variográfico exponencial al semivariograma experimental de velocidad del viento
model_mco_exp4 <- variofit(semivariograma_ws, ini = ini.vals4, cov.model = "exponential")
## 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 "0"     "3"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 1562.282905701
## Warning in variofit(semivariograma_ws, ini = ini.vals4, cov.model =
## "exponential"): unreasonable initial value for sigmasq + nugget (too low)
# Ajustar el modelo variográfico gaussiano al semivariograma experimental de velocidad del viento
model_mco_gaus4 <- variofit(semivariograma_ws, ini = ini.vals4, cov.model = "gauss")
## 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 "1.11"  "9.22" "0"   "0.5"
## status        "est"   "est"  "est" "fix"
## loss value: 602.096829140334
## Warning in variofit(semivariograma_ws, ini = ini.vals4, cov.model = "gauss"):
## unreasonable initial value for sigmasq (too high)
## Warning in variofit(semivariograma_ws, ini = ini.vals4, cov.model = "gauss"):
## unreasonable initial value for sigmasq + nugget (too high)
## Warning in variofit(semivariograma_ws, ini = ini.vals4, cov.model = "gauss"):
## unreasonable initial value for phi (too high)
# Ajustar el modelo variográfico esférico al semivariograma experimental de velocidad del viento
model_mco_spe4 <- variofit(semivariograma_ws, ini = ini.vals4, cov.model = "spherical")
## 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 "0"     "3"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 1562.282905701
## Warning in variofit(semivariograma_ws, ini = ini.vals4, cov.model =
## "spherical"): unreasonable initial value for sigmasq + nugget (too low)
# Graficar el semivariograma experimental y los modelos ajustados
plot(semivariograma_ws, pch = 16, col = "black", main = "Semivariograma Experimental y Modelos Ajustados V. Viendo")
lines(model_mco_exp4, col = "blue")
lines(model_mco_gaus4, col = "red")
lines(model_mco_spe4, col = "purple")
legend("topright", legend = c("Exponencial", "Gaussiano", "Esférico"), col = c("blue", "red", "purple"), lty = 1, cex = 0.8)

Identificación del mejor modelo teórico:

# Crear un data frame para almacenar los resultados de los modelos variográficos
resultados <- data.frame(
  Modelo = c("Exponencial", "Gaussiano", "Esférico"),
  Humedad = NA,
  Temperatura = NA,
  Altitud = NA,
  Viento = NA
)

# Agregar los valores ajustados de los modelos variográficos para la humedad relativa
resultados$Humedad <- c(model_mco_exp1$value, model_mco_gaus1$value, model_mco_spe1$value)

# Agregar los valores ajustados de los modelos variográficos para la temperatura
resultados$Temperatura <- c(model_mco_exp2$value, model_mco_gaus2$value, model_mco_spe2$value)

# Agregar los valores ajustados de los modelos variográficos para la altitud
resultados$Altitud <- c(model_mco_exp3$value, model_mco_gaus3$value, model_mco_spe3$value)

# Agregar los valores ajustados de los modelos variográficos para la velocidad del viento
resultados$Viento <- c(model_mco_exp4$value, model_mco_gaus4$value, model_mco_spe4$value)

# Imprimir los resultados
print(resultados)
##        Modelo  Humedad Temperatura Altitud   Viento
## 1 Exponencial 365359.4    3415.423 4535530 2.945931
## 2   Gaussiano 449110.8    4854.112 6881224 2.945931
## 3    Esférico 365582.9    4041.154 4535396 2.233132

De manera gráfica es un poco difícil discriminar cuál de los modelos se ajusta mejor, por lo cual se construyó una tabla con los valores de los errores medios cuadráticos de cada modelo y su respectiva variable de análisis, logrando así discriminar que los mejores modelos son:

Para la humedad relativa: Modelo exponencial Para la temperatura: Modelo exponencial Para la altitud: Modelo esférico Para la velocidad del viento: Modelo esférico

Predicción espacial mediante el método Kriging

# Crear una cuadrícula de puntos para realizar la interpolación
datos_grid <- expand.grid(Longitude = seq(-276162.5, -276156.8, length = 100),
                          Latitude = seq(8611.564, 8617.082, length = 100))

# Realizar la interpolación de la humedad relativa utilizando el modelo variográfico exponencial
datos_ko1 <- krige.conv(geo_humid, loc = datos_grid,
                        krige = krige.control(nugget = 0, trend.d = "cte", trend.l = "cte", 
                                              cov.pars = c(sigmasq = model_mco_exp1$cov.pars[1],
                                                            phi = model_mco_exp1$cov.pars[2])))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
# Realizar la interpolación de la temperatura utilizando el modelo variográfico exponencial
datos_ko2 <- krige.conv(geo_temp, loc = datos_grid,
                        krige = krige.control(nugget = 0, trend.d = "cte", trend.l = "cte", 
                                              cov.pars = c(sigmasq = model_mco_exp2$cov.pars[1],
                                                            phi = model_mco_exp2$cov.pars[2])))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
# Realizar la interpolación de la altitud utilizando el modelo variográfico esférico
datos_ko3 <- krige.conv(geo_altit, loc = datos_grid,
                        krige = krige.control(nugget = 0, trend.d = "cte", trend.l = "cte", 
                                              cov.pars = c(sigmasq = model_mco_spe3$cov.pars[1],
                                                            phi = model_mco_spe3$cov.pars[2])))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
# Realizar la interpolación de la velocidad del viento utilizando el modelo variográfico esférico
datos_ko4 <- krige.conv(geo_ws, loc = datos_grid,
                        krige = krige.control(nugget = 0, trend.d = "cte", trend.l = "cte", 
                                              cov.pars = c(sigmasq = model_mco_spe4$cov.pars[1],
                                                            phi = model_mco_spe4$cov.pars[2])))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood

Generación de una imagen de interpolación o mapa.

# Graficar el mapa de kriging para la humedad relativa
image(datos_ko1, main = "Kriging Humedad Relativa", xlab = "Longitud (segundos)", ylab = "Latitud (segundos)")

# Graficar el mapa de kriging para la temperatura
image(datos_ko2, main = "Kriging Temperatura", xlab = "Longitud (segundos)", ylab = "Latitud (segundos)")

# Graficar el mapa de kriging para la altitud
image(datos_ko3, main = "Kriging Altitud", xlab = "Longitud (segundos)", ylab = "Latitud (segundos)")

# Graficar el mapa de kriging para la velocidad del viento
image(datos_ko4, main = "Kriging Velocidad del Viento", xlab = "Longitud (segundos)", ylab = "Latitud (segundos)")

Conclusiones

Después de analizar los gráficos de variograma, se observa que la variable de velocidad del viento no muestra correlación espacial evidente. Por lo tanto, al aplicar el método de interpolación Kriging a esta variable, se obtiene un resultado gráfico difícil de interpretar y sin patrones claros.

En contraste, las variables de humedad relativa, temperatura y altitud muestran una correlación espacial significativa, siendo la altitud la más notable. Esta correlación se refleja de manera más adecuada en los gráficos de interpolación, donde se observan patrones más claros.

En conclusión, se identifica que algunas variables pueden tener mediciones en puntos espaciales con coordenadas, pero carecer de correlación espacial. Por lo tanto, la aplicación de métodos de interpolación sobre estas variables puede no producir resultados útiles para el análisis de información espacial. Por otro lado, existen variables con una correlación espacial fuerte o débil, y esta distinción puede determinarse mediante el análisis de variogramas.