Cultivo de Aguacate

Los datos corresponden a las mediciones que se han hecho en árboles de producción de aguacate en una finca del departamento del Cauca.

La base de datos tiene la variable temperatura la cual puede influir en el cultivo de aguacate, pues para este cultivo es necesario rangos optimos de temperatura.

1. Análisis exploratorio

Por lo descrito anteriormente para este caso de estudio nos cetraremos en el analisis de la variable Temperatura y se tendra como fecha de referencia 09-10-2019.

Data$fecha <- as.Date(Data$FORMATTED_DATE_TIME, format = "%d/%m/%Y")
Data_20 <- filter(Data, fecha == "2019-08-21")

tabla <- (head(Data_20,5))

tabla %>%
  kbl() %>%
  kable_paper("hover", 
              full_width = F)
id_arbol Latitude Longitude FORMATTED_DATE_TIME Psychro_Wet_Bulb_Temperature Station_Pressure Relative_Humidity Crosswind Temperature Barometric_Pressure Headwind Direction_True Direction_Mag Wind_Speed Heat_Stress_Index Altitude Dew_Point Density_Altitude Wind_Chill Estado_Fenologico_Predominante Frutos_Afectados fecha
1 2.381290 -76.61264 21/08/2019  9:22:57 a, m, 14.8 805.1 33.6 0.2 25.7 805.0 0.7 166 165 0.8 24.1 1896 8.6 2.743 25.7 715 0 2019-08-21
2 2.381320 -76.61259 21/08/2019  9:27:13 a, m, 11.6 805.2 36.8 3.6 20.8 805.2 3.5 314 313 5.1 19.5 1895 5.5 2.570 20.8 715 3 2019-08-21
3 2.381353 -76.61254 21/08/2019  9:36:36 a, m, 12.9 805.8 31.5 0.4 23.7 805.7 0.7 332 331 0.8 22.0 1889 5.8 2.661 23.7 715 0 2019-08-21
4 2.381374 -76.61261 21/08/2019  9:38:02 a, m, 14.1 805.7 33.2 0.6 25.0 805.7 0.7 139 139 0.9 23.2 1890 7.7 2.707 24.9 715 0 2019-08-21
5 2.381335 -76.61266 21/08/2019  9:39:38 a, m, 14.3 805.2 34.3 0.4 25.0 805.2 0.4 129 128 0.6 23.3 1894 8.1 2.714 24.9 715 1 2019-08-21

Localización de los aborles

leaflet() %>% 
  addTiles() %>% 
  addCircleMarkers(lng = Data_20$Longitude, lat = Data_20$Latitude, radius = 0.2, color = "#006400") %>% 
  addControl(html = "<h1>Arboles de Aguacates</h1>", position = "topright")

Comportamiento de la temperatura.

Temp=as.geodata(Data_20,coords.col = 3:2,data.col = 9)
plot(Temp)

Según lo observado en las gráficas, se evidencian áreas donde la temperatura experimenta un notorio aumento. Este cambio podría tener un impacto relevante en el crecimiento de los aguacates, ya que las condiciones térmicas idóneas son importantes en su desarrollo. Cabe destacar que el rango de temperaturas registrado en la finca oscila entre 20 y 34 grados Celsius para ese día en particular.

# Convertir los datos a un objeto sf
Data_20_sf <- st_as_sf(Data_20, coords = c("Longitude", "Latitude"), crs = 4326, remove = FALSE)


ggplot(Data_20_sf) +
  geom_point(aes(x = Longitude, y = Latitude, color = Temperature), size = 2) +
  scale_color_viridis_c() + 
  labs(title = "Temperatura del cutivo de Aguacate",
       x = "Longitud",
       y = "Latitud",
       color = "Temperatura") +
  theme_minimal()

Hay un sector muy especifico con las temperaturas muy bajas, puede haber correlación espacial con esta variable.

2. Autocorrelacion espacial con el semivariograma

dist_summary_df <- summary(dist(Data_20[,3:2]))
tabla1<- as.data.frame(as.matrix(dist_summary_df))

tabla1 %>%
  kbl(caption = "Distancias entre Coordenadas") %>%
  kable_paper("hover", full_width = F)
Distancias entre Coordenadas
V1
Min. 0.0000330
1st Qu. 0.0004017
Median 0.0006544
Mean 0.0007038
3rd Qu. 0.0009703
Max. 0.0019123

Estos valores indican que la mayoría de las distancias entre los puntos de medición son pequeñas, con una media de 0.0006827.

s_variograma = variog(Temp, option = "bin", uvec = seq(0.0004051, 0.0009178, length.out = 20))
## variog: computing omnidirectional variogram
variograma_temp <- variog.mc.env(Temp, obj.variog = s_variograma, nsim = 99)
## variog.env: generating 99 simulations by permutating data values
## variog.env: computing the empirical variogram for the 99 simulations
## variog.env: computing the envelops
# Definir los rangos de los ejes
x_range <- c(0.0003, 0.0009)
y_range <- c(2, 10)

# Graficar con los rangos de los ejes modificados
plot(s_variograma, main = "Semivariograma con Entorno de Monte Carlo", xlim = x_range, ylim = y_range)
lines(variograma_temp, col = "#EE6AA7", lty = 3)

Al observar la forma ascendente indica que la variabilidad entre puntos cambia con la distancia, sugiriendo que los datos tienen una estructura espacial no homogénea.

3. Identificar el mejor modelo teórico

ajuste del modelo exponencial al semivariograma

# Crear una matriz de valores iniciales para el ajuste
ini.vals <- expand.grid(seq(1, 15, length = 10), seq(6e-04, 8e-05, length = 10))



#modelo exponencial al semivariograma
model_exp <- variofit(s_variograma, ini = ini.vals, cov.model = "exponential", wei = "npair", min = "optim")
## 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 "7.22"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 12367.847265673

Ajuste del modelo gaussiano al semivariograma

# Ajustar el modelo gaussiano al semivariograma
model_gaus <- variofit(s_variograma, ini = ini.vals, cov.model = "gaussian", wei = "npair", min = "optim")
## 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 "5.67"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 11437.4720721014

Ajuste del modelo esférico al semivariograma

model_sph <- variofit(s_variograma, ini = ini.vals, cov.model = "spherical", fix.nug = TRUE, wei = "npair", min = "optim")
## 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 "4.11"  "0"   "0"   "0.5"
## status        "est"   "est" "fix" "fix"
## loss value: 38645.2401367575
plot(s_variograma, main = "Semivariograma con Modelos Ajustados")
lines(model_exp, col = "#EEA9B8")
lines(model_gaus, col = "#836FFF")
lines(model_sph, col = "#54FF9F")

# Agregar leyenda
legend("bottomright", legend = c("Modelo Exponencial", "Modelo Gaussiano", "Modelo Esférico"),
       col = c("#EEA9B8", "#836FFF", "#54FF9F"), lty = 1, cex = 0.6)

Después de comparar los modelos, se observa que el modelo Gaussiano es el que mejor se ajusta a los datos. Por lo tanto, se ha seleccionado este modelo para realizar las predicciones.

4. Prediccion espacial con la metodologia Kriging

##Predicción Espacial Kriging
geodatos_grid=expand.grid( lon=seq(-76.71022,-76.7118,l=80),lat=seq(2.392101 ,2.393634 ,l=80))
plot(geodatos_grid)
points(Data_20[,3:2],col="red")

# Realizar la predicción con los mismos parámetros
prediccion <- krige.conv(Temp, loc = geodatos_grid,
                         krige = krige.control(nugget = 0.5, trend.d = "cte",trend.l = "cte",cov.pars = c(sigmasq = 5.67, phi = 0)))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
image(prediccion, main="kriging Predict", xlab="East", ylab="North")

La totalidad de las predicciones permanece constante. Una hipótesis es basada en la suposición de una variabilidad espacial muy cerca a ser constante en todas las direcciones y distancias, sustentada por un modelo de variograma que exhibe una estructura de correlación uniforme en los datos. Además, si el modelo de variograma y los parámetros de ajuste del kriging se optimizaron de manera precisa para adecuarse a un patrón constante, todas las predicciones podrían converger hacia un valor uniforme.

pred=cbind(geodatos_grid,prediccion$predict)
temp_predict=rasterFromXYZ(cbind(geodatos_grid,prediccion$predict))
plot(temp_predict)

5. Validación Cruzada

valida=xvalid(geodata = Temp,model = model_exp)
## xvalid: number of data locations       = 394
## xvalid: number of validation locations = 394
## 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, 
## xvalid: end of cross-validation
MAE=mean(abs(valida$error))
MAE
## [1] 0.9315476

MAE de aproximadamente 0.9315476 indica que, en promedio, las predicciones del modelo se desvían de los valores observados por alrededor de 0.9315476 grados.

Aunque las predicciones se mantienen constantes, el modelo arroja un error absoluto medio (MAE) reducido, lo cual podría atribuirse a la estabilidad observada en los datos de temperatura.