Introduccion

En este script tiene como propósito explorar datos espaciales aplicando técnicas de geoestadística para analizar la producción de aguacate en un grupo específico de árboles en la región del Cauca. El objetivo es determinar si existe una correlación espacial entre la variable de temperatura y la producción de aguacate.

Datos_Completos_Aguacate <- read_excel("C:/Users/aleja/Downloads/Datos_Completos_Aguacate.xlsx", sheet = "Hoja2")
head(Datos_Completos_Aguacate)
## # A tibble: 6 × 21
##   id_arbol Latitude Longitude FORMATTED_DATE_TIME        Psychro_Wet_Bulb_Temp…¹
##      <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
## # ℹ 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>

Grafico de la ubicacion de los arboles

leaflet() %>% addTiles() %>% addCircleMarkers(lng = Datos_Completos_Aguacate$Longitude,lat = Datos_Completos_Aguacate$Latitude,radius = 0.2,color = "red")

Geoestadística

Datos_Completos_Aguacate_temp=as.geodata(Datos_Completos_Aguacate,coords.col = 3:2,data.col = 9)
plot(Datos_Completos_Aguacate_temp)

summary(dist(Datos_Completos_Aguacate[,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
variograma=variog(Datos_Completos_Aguacate_temp,option = "bin",uvec=seq(0,0.001,0.00002))
## variog: computing omnidirectional variogram
datos.env=variog.mc.env(Datos_Completos_Aguacate_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
plot(variograma)
lines(datos.env)

Al observar las imágenes anteriores, se puede notar que la temperatura varía dentro de un rango específico, generalmente entre 22 y 30 grados Celsius. Además, se puede apreciar cómo los árboles se agrupan según sus temperaturas similares, lo que sugiere una distribución espacial de las temperaturas en el entorno forestal. La agrupación de los puntos o árboles con temperaturas similares en los gráficos de dispersión sugiere que existe una correlación entre la temperatura y la ubicación geográfica de los árboles. Es posible que los árboles con temperaturas cercanas tiendan a agruparse en áreas específicas del territorio, lo que podría indicar preferencias climáticas o condiciones ambientales particulares para ciertas especies de árboles.

Ajuste del modelo teorico

ini.vals = expand.grid(seq(1.2,1.5,l=10), seq(0.0001,0.0008,l=10))

Modelo exponencial

model_mco_exp=variofit(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 "1.5"   "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 296807.231725087
model_mco_exp
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: exponential
## parameter estimates:
##   tausq sigmasq     phi 
##  0.7969  2.2860  0.0001 
## Practical Range with cor=0.05 for asymptotic range: 0.000232806
## 
## variofit: minimised weighted sum of squares = 9158.212

Modelo Gaussiano

model_mco_gaus=variofit(variograma, ini=ini.vals, cov.model="gaussian", wei="npair", min="optim",nugget = 0)
## 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.5"   "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 288205.86784414
model_mco_gaus
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: gaussian
## parameter estimates:
##   tausq sigmasq     phi 
##  0.7902  2.2842  0.0001 
## Practical Range with cor=0.05 for asymptotic range: 0.0002330932
## 
## variofit: minimised weighted sum of squares = 8947.68

Modelo esferico

model_mco_spe=variofit(variograma, ini=ini.vals, cov.model="spheric", 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 "1.5"   "0"   "0"   "0.5"
## status        "est"   "est" "fix" "fix"
## loss value: 283743.808438713
model_mco_spe
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: spherical
## fixed value for tausq =  0 
## parameter estimates:
## sigmasq     phi 
##  3.0186  0.0000 
## Practical Range with cor=0.05 for asymptotic range: 0
## 
## variofit: minimised weighted sum of squares = 17731.76
plot(variograma)

lines(model_mco_exp,col="blue")
lines(model_mco_gaus,col="red")
lines(model_mco_spe,col="green")

Al analizar tanto visualmente las gráficas como el valor del Error Cuadrático Estándar (SCE), se concluye que el modelo gaussiano ofrece el mejor ajuste a los puntos del semivariograma muestral. A pesar de que visualmente el modelo gaussiano no sigue exactamente la tendencia de los puntos de datos, el valor del SCE indica que es el modelo que mejor se ajusta a los datos en términos de minimizar los errores cuadráticos.

Al evaluar los tres modelos considerados, se observa que el modelo esférico exhibe la menor suma de cuadrados de los errores (SCE), con un valor específico de 17731. Esto significa que, a pesar de que el modelo gaussiano se ajusta visualmente mejor, el modelo esférico logra minimizar de manera más efectiva los errores cuadráticos en la estimación de la semivarianza a diferentes distancias

Interpolacion (prediccion espacial)

c(min(Datos_Completos_Aguacate[,3]), max(Datos_Completos_Aguacate[,3]), min(Datos_Completos_Aguacate[,2]), max(Datos_Completos_Aguacate[,2]))
## [1] -76.711799 -76.710215   2.392101   2.393634

Prediccion espacial Kriking

aguacate_grid=expand.grid( lon=seq(-76.710215,-76.711799,l=100),lat=seq(2.392101 ,2.393634 ,l=100))
plot(aguacate_grid)
points(Datos_Completos_Aguacate_temp$coords,col="red")

geodatos_ko=krige.conv(Datos_Completos_Aguacate_temp, loc=aguacate_grid,
                       krige= krige.control(nugget=0,trend.d="cte", 
                                            trend.l="cte",cov.pars=c(sigmasq=3.0186, phi=0.0001 )))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
par(mfrow=c(1,2))
image(geodatos_ko, main="kriging Predict", xlab="East", ylab="North")
contour(geodatos_ko,main="kriging Predict", add=TRUE, drawlabels=TRUE)

Transformacion de la imagen en raster

pred=cbind(aguacate_grid,geodatos_ko$predict)
temp_predict=rasterFromXYZ(cbind(aguacate_grid,geodatos_ko$predict))
levelplot(temp_predict,par.settings =BuRdTheme)

Realizacion de la validacion cruzada

valida=xvalid(geodata = Datos_Completos_Aguacate_temp,model = model_mco_gaus)
## 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
MAE=mean(abs(valida$error))
MAE
## [1] 0.8033426

El modelo presenta un error de predicción de acuerdo con la validación cruzada de 0.8033426 grados.