Caso Predicciones de temperatura en cultivos de aguacate
Oscar Eduardo Ramirez Zambrano
2025-05-18
Introducción
A partir de información sobre las variables climatológicas de una finca en el departamento de Cauca (Colombia), se realizará un análisis geoestadístico con el fin de realizar predicciones sobre estas variables, para esto se:
Evaluará la autocorrelación espacial con el semivariograma. Identificaremos el mejor modelo teórico. Realizaremos una predicción espacial usando la metodología krigin. Generamos imagen de la predicción espacial.
Reto
El reto principal consiste en determinar el modelo que presente el mejor ajuste y en base a este realizar predicciones de temperatura usando el interpolador de temperatura Kriging.
1. Análisis explotarorio de los datos
A partir de la base de datos “Datos_Completos_Aguacate.xlsx” que contiene información tomaremos el filtro del periodo“01/10/2020 10:11:12 a.m.” usando la variable “Temperature”.
Carga de librerias necesarias
require (readxl)
require(geoR)
require(raster)
require(rasterVis)
require(leaflet)
require(dplyr)Carga del archivo de la base de datos
Cargamos la base de datos filtrando los datos necesarios.
datos = read_excel("C:/Users/ORAM/OneDrive/Desktop/Datos_Completos_Aguacate.xlsx")
head(datos)# A tibble: 6 × 21
id_arbol Latitude Longitude FORMATTED_DATE_TIME Psychro_Wet_Bulb_Tempe…¹
<chr> <dbl> <dbl> <chr> <dbl>
1 1 2.38 -76.6 21/08/2019 9:22:57 a, m, 14.8
2 2 2.38 -76.6 21/08/2019 9:27:13 a, m, 11.6
3 3 2.38 -76.6 21/08/2019 9:36:36 a, m, 12.9
4 4 2.38 -76.6 21/08/2019 9:38:02 a, m, 14.1
5 5 2.38 -76.6 21/08/2019 9:39:38 a, m, 14.3
6 6 2.38 -76.6 21/08/2019 9:42:02 a, m, 14.2
# ℹ 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>En los registros se observan las variables climatológicas de cada árbol: humedad relativa, temperatura, dirección del viento, entre otras.
leaflet() %>% addTiles() %>% addCircleMarkers(lng = datos$Longitude,lat = datos$Latitude,radius = 0.2,color = "red")De acuerdo con la visualización anterior, existen 4 lotes donde se encuentran los cultivos de aguacate cercanos al municipio de Timbio, Cauca
a continuación, filtramos para el periodo 01/10/2020 10:11:12 a.m.
datos$fecha <- as.Date(datos$FORMATTED_DATE_TIME, format = "%d/%m/%Y")
# Filtrar los datos por la fecha específica
datos_2020 <- datos %>%
filter(FORMATTED_DATE_TIME == "01/10/2020 10:11:12 a, m,")
# Muestra los datos filtrados
head(datos_2020)# A tibble: 0 × 22
# ℹ 22 variables: id_arbol <chr>, Latitude <dbl>, Longitude <dbl>,
# FORMATTED_DATE_TIME <chr>, Psychro_Wet_Bulb_Temperature <dbl>,
# 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>, …Observamos que el filtro no nos ha arrojado ningún resultado, por lo tanto procedemos a verificar.
# Calculando datos faltantes
Datos_faltantes <- data.frame(colnames(datos), sapply(datos, function(x) sum(is.na(x))))
Datos_faltantes <- Datos_faltantes[c(2:4, 9), ]
# Limpiando Dataframe
rownames(Datos_faltantes) <- NULL
colnames(Datos_faltantes) <- c("Variable", "Datos Faltantes")
# Presentando información en formato tabla
library(kableExtra) # Asegúrate de tener cargada esta librería
kable(Datos_faltantes, caption = "<center><b>Tabla 1. Datos Faltantes por Variable</b></center>") %>%
kable_classic(full_width = FALSE)| Variable | Datos Faltantes |
|---|---|
| Latitude | 0 |
| Longitude | 0 |
| FORMATTED_DATE_TIME | 0 |
| Temperature | 0 |
Observamos que ninguna de las variables tiene datos faltantes
# Encontramos registros que contengan el dato
print(head(datos[grepl(
"01/10/2020",datos$FORMATTED_DATE_TIME), ][c(1:4)], 5))# A tibble: 5 × 4
id_arbol Latitude Longitude FORMATTED_DATE_TIME
<chr> <dbl> <dbl> <chr>
1 1 2.39 -76.7 01/10/2020 10:11:12 a, m,
2 2 2.39 -76.7 01/10/2020 10:11:12 a, m,
3 3 2.39 -76.7 01/10/2020 10:11:12 a, m,
4 4 2.39 -76.7 01/10/2020 10:11:12 a, m,
5 5 2.39 -76.7 01/10/2020 10:11:12 a, m,Se observa que solo 5 registros cumplen la condición, por lo tanto procedemos a validar las diferencias entre las dos cadenas de caracteres.
# Verificamos los valores hexadecimales de ambas cadenas de texto
charToRaw(datos$FORMATTED_DATE_TIME[9140]) [1] 30 31 2f 31 30 2f 32 30 32 30 c2 a0 20 31 30 3a 31 31 3a 31 32 20 61 2c 20
[26] 6d 2ccharToRaw("01/10/2020 10:11:12 a, m,") [1] 30 31 2f 31 30 2f 32 30 32 30 20 20 31 30 3a 31 31 3a 31 32 20 61 2c 20 6d
[26] 2cObservamos que se tienen difercias entre las dos cadenas de caracteres, procedemos a la limpieza:
# Eliminar espacios
datos_fix <- datos
datos_fix$FORMATTED_DATE_TIME <- gsub(
"\u00A0", " ", datos_fix$FORMATTED_DATE_TIME)# Filtramos nuevamente la cadena de texto solicitada
datos_2020 <- datos_fix[
datos_fix$FORMATTED_DATE_TIME == "01/10/2020 10:11:12 a, m,",]
print(datos_2020)# A tibble: 534 × 22
id_arbol Latitude Longitude FORMATTED_DATE_TIME Psychro_Wet_Bulb_Tem…¹
<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
7 7 2.39 -76.7 01/10/2020 10:11:12 a, m, 22.2
8 8 2.39 -76.7 01/10/2020 10:11:12 a, m, 22.6
9 9 2.39 -76.7 01/10/2020 10:11:12 a, m, 23
10 10 2.39 -76.7 01/10/2020 10:11:12 a, m, 22.7
# ℹ 524 more rows
# ℹ abbreviated name: ¹Psychro_Wet_Bulb_Temperature
# ℹ 17 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>, …Obtenemos de esta manera los datos limpios, procedemos a realizar la conversión a geodata.
# Convertir dataframe en geodata
datos_geo <- as.geodata(
datos_2020,coords.col = 3:2,data.col = 9)
# Gráficar geodata
par(mfrow = c(2, 1), mar = c(5, 4, 1, 4)) # Ajustar márgenes: 5 en la parte inferior, 1 en la superior
# Graficar el objeto geodata
plot(datos_geo)
# Título: Añadir texto sobre el gráfico
mtext("Gráfico 1. Objeto geodata de la\nvariable Temperatura",
side = 3, line = 0.1, cex = 1 )
2. Evalución del Semivariograma
Determinamos determinar las distancias a graficar, normalmente lo recomendado es hasta el tercer cuartil de las distancias entre los puntos a evaluar, así que calculamos las distancias y verificamos.
# Calculamos las estadísticas de los puntos de estudio
summary(dist(datos_2020[3:2])) Min. 1st Qu. Median Mean 3rd Qu. Max.
1.712e-05 4.051e-04 6.408e-04 6.827e-04 9.178e-04 1.959e-03 Tomando el dato del tercer cuartil.
# Creamos el variograma
datos_variog <- variog(datos_geo,option = "bin",
uvec=seq(0,9.178e-04,9.178e-05), messages = FALSE)
# Calcular margen de correlación
datos_marks <- variog.mc.env(datos_geo,
obj=datos_variog, messages = FALSE)
# Graficar semivariograma
plot(datos_variog, main = "Gráfico 2. Semivariograma de Temperatura",
xlab = "Distancia", ylab = "Semivarianza")
lines(datos_marks)
Con la gráfica del semivariograma podemos observar que temperatura tiene una correlación espacial, sin embargo analizaremos otras variables como son la humedad relativa y la velocidad del viento.
# Humedad relativa
# Convertir dataframe en geodata
datos_geo2 <- as.geodata(
datos_2020,coords.col = 3:2,data.col = 7)
# Crear variograma
datos_variog2 <- variog(datos_geo2,option = "bin",
uvec=seq(0,9.178e-04,9.178e-05), messages = FALSE)
# Calcular margen de correlación
datos_marks2 <- variog.mc.env(datos_geo2,
obj=datos_variog2, messages = FALSE)
# Graficar semivariograma
plot(datos_variog2, main = "Gráfico 3. Semivariograma de Humedad Relativa",
xlab = "Distancia", ylab = "Semivarianza")
lines(datos_marks2)
# Velocidad del viento
# Convertir dataframe en geodata
datos_geo2 <- as.geodata(
datos_2020,coords.col = 3:2,data.col = 14)
# Crear variograma
datos_variog2 <- variog(datos_geo2,option = "bin",
uvec=seq(0,9.178e-04,9.178e-05), messages = FALSE)
# Calcular margen de correlación
datos_marks2 <- variog.mc.env(datos_geo2,
obj=datos_variog2, messages = FALSE)
# Graficar semivariograma
plot(datos_variog2, main = "Gráfico 4. Semivariograma de Velocidad del Viento",
xlab = "Distancia", ylab = "Semivarianza")
lines(datos_marks2)
Se observa que la velocidad del viento no tiene una correlación espacial ya que la mayoría de sus puntos caen en la zona de no correlación.
3. Identificar el Modelo
Escogiendo la variable de temperatura, verificamoslos modelos para determinar el grado de ajuste con los datos.
#Fijar rango de evaluación
val_ini <- expand.grid(seq(1.5,3,l=10), seq(4e-04,6e-04,l=10))
# Gráficar semiovariograma junto a las curvas teóricas
plot(datos_variog, main = "Gráfico 6. Semivariograma de Temperatura\ncon modelos teóricos",
xlab = "Distancia", ylab = "Semivarianza")
datos_exp <- variofit(datos_variog, ini=val_ini,
cov.model="exponential", wei="npair",
min="optim", messages = FALSE)
datos_gauss <- variofit(datos_variog, ini=val_ini,
cov.model="gaussian", wei="npair",
min="optim",nugget = 0, messages = FALSE)
datos_spe <- variofit(datos_variog, ini=val_ini,
cov.model="spheric",fix.nug=TRUE,
wei="npair", min="optim", messages = FALSE)
lines(datos_exp,col="blue")
lines(datos_gauss,col="red")
lines(datos_spe,col="yellow")
Se observa que el modelo exponencial se ajusta mejor a los datos, verificamos la suma de errores al cuadrado para tener un indicador numérico de ajuste.
#Crear dataframe de la summa de errores cuadrados por modelo
Errores_cuadrado <- data.frame(Modelo=c("Exponencial", "Gaussiano", "Esférico"),
Suma_Errores=c(datos_exp$value, datos_gauss$value, datos_spe$value))
#Crear tabla para visualizar la suma de errores cuadrados
kbl(Errores_cuadrado, caption = "<center><b>Tabla 2. Suma de Errores al cuadrado por Modelo</b></center>", col.names=c("Modelo","Suma de errores al cuadrado"))%>%
kable_classic(full_width = F)| Modelo | Suma de errores al cuadrado |
|---|---|
| Exponencial | 4857.673 |
| Gaussiano | 14947.769 |
| Esférico | 6583.649 |
El modelo exponencial tiene un mejor desempeño, seguido del esférico, el modelo gaussiano no se comporta bien para esta variable, es posible que esto se deba a la baja correlación espacial del modelo.
4. Realizar Predicción espacial
Necesitamos obtener el área donde va a ser graficada la predicción, para ello obtenemos el contorno obtenidos por los puntos extremos de latitud y longitud en la base de datos.
#Crear dataframe del área a gráficar
area_val <- data.frame(Modelo=c("Min. Long", "Max Long.", "Min. Lat", "Max. Lat"),
Suma_Errores=c(min(datos_2020[3]), max(datos_2020[3]), min(datos_2020[2]), max(datos_2020[2])))
#Crear tabla para visualizar los valores del área
kbl(area_val, caption = "<center><b>Tabla 3. Rango de área a graficar</b></center>", col.names=c("Línea","Valor"))%>%
kable_classic(full_width = F)| Línea | Valor |
|---|---|
| Min. Long | -76.711799 |
| Max Long. | -76.710215 |
| Min. Lat | 2.392101 |
| Max. Lat | 2.393634 |
Graficamos la “grilla”, donde se realizará la predicción espacial.
# Graficar grilla de terreno
datos_grilla <- expand.grid(lon=seq(-76.7118, -76.71022, l=100),
lat=seq(2.392101, 2.393634, l=100))
plot(datos_grilla, main = "Gráfico 7. Grilla del terreno",
xlab = "Longitud", ylab = "Latitud")
points(datos_2020[3:2],col="cyan")
Procedemos a realizar la predicción espacial de la temperatura del terreno.
# Crear predicción, pasando modelo
salida <- output.control(messages=FALSE)
datos_pred <- krige.conv(
datos_geo, loc=datos_grilla, output = salida,
krige= krige.control(obj.model = datos_exp))
# Gráficar predicción
par(mfrow=c(1,2))
image(datos_pred, main="Gráfico 8. Predicción Espacial",
xlab="Longitud", ylab="Latitud")
par(new=TRUE)
contour(datos_pred, xlab="Longitud", ylab="Latitud", drawlabels=TRUE)
Generamos una visualización gráfica espacial del error de la predicción de los datos.
# Gráficar desviación del error predicción
par(mfrow=c(1,2))
image(datos_pred, main="Gráfico 9. Error de la Predicción\nEspacial",
val=sqrt(datos_pred$krige.var), xlab="Longitud", ylab="Latitud")
par(new=TRUE)
contour(datos_pred, xlab="Longitud",
val=sqrt(datos_pred$krige.var), ylab="Latitud", drawlabels=TRUE)
5. Generar Raster
Generamos imágenes raster a partir de las predicciones realizadas.
# Creando mapa raster a partir de los datos
datos_tempPred <- rasterFromXYZ(cbind(datos_grilla, datos_pred$predict))
plot(datos_tempPred, main="Gráfico 10. Mapa raster de temperatura",
xlab="Longitud", ylab="Latitud")
Generamos imágenes raster a partir de los errores de las predicciones realizadas.
# Creando mapa raster del error a partir de los datos
datos_errorPred <- rasterFromXYZ(cbind(datos_grilla, val=sqrt(datos_pred$krige.var)))
plot(datos_errorPred, main="Gráfico 11. Mapa raster del error de la predicción",
xlab="Longitud", ylab="Latitud")
Realizamos una validación del modelo exponencial mediante validación cruzada.
datos_valid=xvalid(geodata = datos_geo,
model = datos_exp)xvalid: number of data locations = 534
xvalid: number of validation locations = 534
xvalid: performing cross-validation at location ... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260, 261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275, 276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 336, 337, 338, 339, 340, 341, 342, 343, 344, 345, 346, 347, 348, 349, 350, 351, 352, 353, 354, 355, 356, 357, 358, 359, 360, 361, 362, 363, 364, 365, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, 376, 377, 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389, 390, 391, 392, 393, 394, 395, 396, 397, 398, 399, 400, 401, 402, 403, 404, 405, 406, 407, 408, 409, 410, 411, 412, 413, 414, 415, 416, 417, 418, 419, 420, 421, 422, 423, 424, 425, 426, 427, 428, 429, 430, 431, 432, 433, 434, 435, 436, 437, 438, 439, 440, 441, 442, 443, 444, 445, 446, 447, 448, 449, 450, 451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 461, 462, 463, 464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 474, 475, 476, 477, 478, 479, 480, 481, 482, 483, 484, 485, 486, 487, 488, 489, 490, 491, 492, 493, 494, 495, 496, 497, 498, 499, 500, 501, 502, 503, 504, 505, 506, 507, 508, 509, 510, 511, 512, 513, 514, 515, 516, 517, 518, 519, 520, 521, 522, 523, 524, 525, 526, 527, 528, 529, 530, 531, 532, 533, 534,
xvalid: end of cross-validationCalculmos el error absoluto medio para dar un indicador general del error del modelo.
MAE=mean(abs(datos_valid$error))
print(MAE)[1] 0.7909316El modelo presenta un error de predicciónde 0.79 grados de acuerdo con la validación cruzada, se puede decir solo considerando este indicador que el modelo tuvo un buen ajuste.
6. Conclusiones
A partir de los resultados obtenidos podemos señalar las siguientes conclusiones:
la temperatura no presenta un patrón estacionario, no se distribuye en forma de gradiente uniforme.
La correlación de temperatura encontrada no es muy fuerte, nos se observa una buenacorrelación espacial, sin embargo en comparación con la humedad relativa y la velocidad del viento es la que mejor correlaciona, esto puede deberse a otros factores no evaluados.
El modelo que presentó el mejor ajuste fue el modelo exponencial, al final en la validación cruzada su error absoluto medio en la temperatura con respecto del valor encontrado en cada punto evaluado fue de 0.79 grados.