CONTEXTO

Este estudio se enfoca en la producción de aguacate de una finca ubicada en el departamento del Cauca, Colombia. Se georreferenciaron 534 árboles con coordenadas de latitud y longitud, y se recopilaron datos transversales el 1 de octubre de 2020. Utilizando sensores portátiles distribuidos por toda la plantación, se midieron diversas variables ambientales, incluida la temperatura.

El objetivo principal es examinar cómo la temperatura influye en el cultivo a gran escala. Para ello, se empleará la geoestadística para interpolar esta variable ambiental clave. Se ajustará un modelo que considere la correlación espacial dentro del área de cultivo y se aplicará kriging para generar un mapa de predicciones espaciales. Este mapa se rasterizará, proporcionando una imagen detallada que facilitará futuras comparaciones y el desarrollo de métricas adicionales. Estas métricas serán fundamentales para comprender el impacto de la temperatura en la productividad del aguacate.

PAQUETES

El análisis se llevará a cabo utilizando los siguientes paquetes estadísticos:

library(tidyverse)
library(readxl)
library(geoR)
library(raster)
library(ggplot2)
library(leaflet)
require(rasterVis)

DATOS

Los datos se importan utilizando la siguiente función:

GeoData <- read_excel("C:/Users/Andre/OneDrive/Escritorio/Datos geograficos y espaciales/Casos/Datos_Completos_Aguacate.xlsx", sheet = "Hoja2")
GeoData
## # A tibble: 534 × 21
##    ID_ARBOL LATITUDE LONGITUDE FORMATTED_DATE_TIME        Psychro_Wet_Bulb_Tem…¹
##       <dbl>    <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
## # ℹ 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>, …

El conjunto de datos consta de 21 variables. Para el propósito de este análisis, las variables Latitud, Longitud y Temperatura son de especial interés debido a su relevancia en el estudio geoespacial del comportamiento de la temperatura en la zona de cultivo.

ESTRUCTURA ESPACIAL

Se utiliza la biblioteca leaflet para crear un mapa interactivo. Este mapa sirve como una herramienta visual clave para representar geográficamente la distribución de los árboles de aguacate dentro de la finca.

leaflet() %>% addTiles() %>% addCircleMarkers(lng = GeoData$LONGITUDE,lat = GeoData$LATITUDE,radius = 0.2,color = "black")

Los marcadores circulares, definidos por las coordenadas de longitud y latitud de cada árbol, se superponen sobre una capa base de mapas de ‘tiles’. El color negro y el radio pequeño de 0.2 se utilizan para marcar cada ubicación de manera discreta pero clara. Esta visualización es fundamental para el análisis geoestadístico, ya que permite una comprensión inmediata de la distribución espacial y facilita la identificación de patrones o anomalías que podrían influir en la temperatura y, por ende, en la producción del cultivo.

geo_temp=as.geodata(GeoData,coords.col = 3:2,data.col = 9)
plot(geo_temp)

El análisis gráfico revela que, en términos espaciales, la finca presenta agrupaciones de temperaturas similares en áreas específicas, identificándose al menos cuatro conglomerados térmicos. Según el histograma, estas agrupaciones varían entre los 22 y 30 grados centígrados. Al examinar la relación entre los datos de temperatura y las coordenadas geográficas, se nota una dispersión de puntos que sugiere una posible ausencia de correlación espacial significativa. En particular, la comparación con la latitud no muestra una concentración notable de puntos, a diferencia de la longitud, donde se percibe una mayor aglomeración. Esto podría indicar diferencias en la influencia de las coordenadas geográficas sobre las variaciones de temperatura registradas en la finca.

summary(dist(GeoData[,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

El resumen estadístico generado por summary(dist(GeoData[,3:2])) proporciona una visión cuantitativa de las distancias entre los puntos georreferenciados de los árboles de aguacate en la finca. Este análisis de distancia es crucial para entender la estructura espacial del cultivo.

Resepcto a la métricas que se obtienen:

La distancia más corta entre dos árboles es aproximadamente \[1.712 \times 10^{-5}\] unidades. Esto indica que los árboles más cercanos están muy próximos entre sí. Luego, el 25% de todas las distancias entre árboles son menores o iguales a \[4.051 \times 10^{-4}\] unidades, lo que sugiere una agrupación relativamente densa de algunos árboles. Máximo (Max.): La distancia más larga entre dos árboles es de \[1.959 \times 10^{-3}\] unidades, lo que muestra la extensión máxima del cultivo. Sin embargo como las medidas de distancias no están en metros si no en grados, su interpretación es más complicada de lo normal.

CORRELACIÓN ESPACIAL

variograma = variog(geo_temp, option = "bin", uvec = seq(0, 9.178e-04, length.out = 30))
## variog: computing omnidirectional variogram
plot(variograma, ylim=c(0,5.5))
variograma_mc=variog.mc.env(geo_temp,obj=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
lines(variograma_mc)

Las bandas que se crean en el gráfico al ejecutar lines(variograma_mc) son conocidas como bandas de confianza del variograma o envolventes de Monte Carlo. Estas bandas proporcionan un rango de variabilidad en el que se espera que los valores del variograma empírico caigan con una cierta probabilidad, basada en la simulación de Monte Carlo. Sirven para evaluar la significancia estadística de la estructura espacial observada en los datos y para ayudar a determinar si la correlación espacial percibida es más fuerte o más débil de lo que se podría esperar por variación aleatoria.

La inspección de la gráfica revela que el semivariograma intersecta las bandas de confianza solo marginalmente, lo que sugiere que, conforme a las predicciones, la correlación espacial de la temperatura en el contexto del cultivo es, en efecto, débil. Este hallazgo es indicativo de que otros factores, posiblemente no espaciales, podrían estar influyendo en las variaciones de temperatura observadas en la plantación.

MODELACIÓN

Para optimizar la convergencia de los modelos geoestadísticos, es esencial definir valores iniciales adecuados. Estos valores actúan como puntos de partida en la estimación y son cruciales para los siguientes modelos teóricos de variogramas:

La selección cuidadosa de estos modelos y sus parámetros iniciales es fundamental para capturar con precisión la estructura espacial inherente a los datos y asegurar la validez de las predicciones espaciales resultantes. Los valores elegidos se muestran a continuación:

ini.vals = expand.grid(seq(2,6,l=30), seq(6e-04,9e-04,l=30))

Tras establecer los valores iniciales, se procede con la calibración de los modelos teóricos. Este paso es crucial para afinar los parámetros que permitirán una representación precisa de la correlación espacial en los datos analizados.

## Modelo Exponencial
model_mco_exp=variofit(variograma, ini=ini.vals, cov.model="exponential",
                       wei="npair", min="optim", nugget = 1.7)
## 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.28"  "0"   "1.7" "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 3451.76632415449
## Modelo Gausiano
model_mco_gaus=variofit(variograma, ini=ini.vals, cov.model="gaussian", wei="npair", min="optim", nugget = 2.38)
## 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"     "0"   "2.38" "0.5"
## status        "est"   "est" "est"  "fix"
## loss value: 5502.66333255951
## Modelo Esférico
model_mco_spe=variofit(variograma, ini=ini.vals, cov.model="spheric",fix.nug=TRUE, wei="npair", min="optim", nugget = 1.28)
## 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"     "0"   "1.28" "0.5"
## status        "est"   "est" "fix"  "fix"
## loss value: 4967.17511414114

Tras la calibración de los modelos, se procede a su representación gráfica junto al semivariograma empírico. Este paso visual es esencial para evaluar la adecuación de los modelos y su capacidad para reflejar la estructura espacial de los datos.

plot(variograma)
lines(model_mco_exp,col="blue")
lines(model_mco_gaus,col="red")
lines(model_mco_spe,col="purple")

Finalmente, se considera el mínimo cuadrado ponderado (minimised weighted sum of squares) como un indicador clave de ajuste. Este valor cuantifica la discrepancia entre el semivariograma empírico y el modelo teórico; cuanto menor sea el valor, mejor será el ajuste del modelo a los datos.

Comparando los valores de mínimos cuadrados ponderados para cada modelo:

El modelo Exponencial tiene el valor más bajo de mínimos cuadrados ponderados, lo que indica que es el que mejor se ajusta al semivariograma empírico de los datos proporcionados. Por lo tanto, basándose en este criterio, el modelo Exponencial sería el más adecuado para representar la correlación espacial de la variable temperatura en el área de estudio.

KRIGING

El proceso comienza con la creación de una malla predictiva que se extiende sobre el área georreferenciada del cultivo de aguacate.

##Predicción Espacial Kriging
geodatos_grid=expand.grid( lon=seq(-76.71018,-76.7119,l=100),lat=seq(2.393640 ,2.39210 ,l=100))
plot(geodatos_grid)
points(GeoData[,3:2],col="red")

Posteriormente, se realiza la predicción utilizando los parámetros del modelo exponencial previamente ajustado.

geodatos_ko=krige.conv(geo_temp, loc=geodatos_grid,
      krige= krige.control(nugget=1.7,trend.d="cte", 
      trend.l="cte",cov.pars=c(sigmasq=2.2759, phi=0.0006 )))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood

A partir de los cálculos de kriging, se genera una representación visual que predice las temperaturas en la zona del cultivo.

image(geodatos_ko, main="Kriging Predict", xlab="Este", ylab="Norte")

Además, se examina una imagen que ilustra la desviación estándar de las temperaturas registradas en el cultivo, proporcionando una medida de incertidumbre en las predicciones.

image(geodatos_ko, main="Kriging StDv Predicted",val=sqrt(geodatos_ko$krige.var), xlab="Este", ylab="Norte")

RASTER

Para aprovechar al máximo las predicciones de kriging, se convierten en una imagen raster. A continuación, se procede a la creación de dicha imagen:

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

La función levelplot se utiliza para interpretar las variaciones de temperatura. En la visualización resultante, las áreas azules indican temperaturas más bajas, mientras que las rojas señalan las zonas donde se detectaron las temperaturas más altas.

levelplot(temp_predict,par.settings =BuRdTheme)

Se observa que la variabilidad de temperatura dentro del cultivo oscila entre 1.4 y 2 grados centígrados, con las temperaturas alrededor de los 24 grados siendo las más frecuentes.

temp_error=rasterFromXYZ(cbind(geodatos_grid,sqrt(geodatos_ko$krige.var)))
levelplot(temp_error,par.settings =BuRdTheme)

Para una comprensión más profunda de la variación térmica, se superpone la imagen raster sobre los puntos georreferenciados del cultivo y se añade una capa de leaflet. Esto no solo enriquece la interactividad del análisis, sino que también permite considerar factores externos que podrían influir en el cultivo debido a su ubicación geográfica.

crs(temp_predict)=crs("+proj=longlat +datum=WGS84")

leaflet() %>% addTiles() %>% addRasterImage(temp_predict,opacity = 0.6) %>% addCircleMarkers(lng = GeoData$LONGITUDE,lat = GeoData$LATITUDE,radius = 0.2,color = "black")

Al contrario de las imagenes anteriores, las zonas rojas en realidad son las más frías y las zonas azules las más calurosas.

VALIDACIÓN CRUZADA

valida=xvalid(geodata = geo_temp,model = model_mco_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-validation

Para la evaluación del modelo se realizá una validación cruzada con el modelo exponencial sobre los datos geoespaciales abarcando un total de 534 ubicaciones tanto para los datos como para las localizaciones de validación.

MAE=mean(abs(valida$error))
MAE
## [1] 0.9863183

El resultado de la validación cruzada se cuantifica a través del Error Absoluto Medio (MAE), que en este caso es de aproximadamente 0.9863 grados centígrados. Este valor representa la diferencia media entre los valores observados y los valores predichos por el modelo. Un MAE más bajo indica un mejor ajuste del modelo a los datos, y aunque este valor no es particularmente bajo, proporciona una medida de la precisión del modelo en la predicción de la temperatura en el cultivo.

CONCLUSIONES

El estudio de los árboles de aguacate ha revelado patrones espaciales significativos en la distribución de la temperatura, destacando la existencia de zonas con temperaturas agrupadas que podrían tener implicaciones directas en la productividad del cultivo. A través de técnicas de kriging y la creación de imágenes raster, se ha logrado una representación detallada de la variabilidad térmica, lo que permite una mejor comprensión de las condiciones ambientales que enfrentan los árboles. El Error Absoluto Medio (MAE) obtenido sugiere que, aunque hay espacio para mejorar la precisión del modelo, las predicciones actuales ofrecen una base sólida para futuras investigaciones y la toma de decisiones agronómicas, siendo un modelo con un margen de error de aproximadamente un grado centígrado aproximadamente.

Es importante señalar que se han registrado temperaturas cercanas a los 30 grados centígrados, las cuales impactan solo una sección del cultivo. Para comprender mejor las causas subyacentes de estas variaciones térmicas, sería prudente examinar otros factores, como la altitud o la altura de los árboles, que podrían influir en las condiciones climáticas de estas áreas específicas del cultivo.