Analisis del cultivo de aguacate y sus variables

El presente trabajo se desarrolla teniendo presente la importancia de analizar la variable de temperatura en el proceso de cultivo de aguacate, teniendo en cuenta información de una finca ubicada en el departamento del Cauca

Inicialmente, se hace la importación del dataset a partir de lo que se genera una geodata que permita realizar una valoración desde la estadística geoespcial.

Importación dataset

library(readxl)
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
## --------------------------------------------------------------
library(raster)
## Loading required package: sp
library(ggplot2)
library(leaflet)

geo_datos<-read_excel("/Users/juan/Desktop/MAESTRIA/Estadistica espacial/ejercicio agucate.xlsx")
knitr::kable(head(geo_datos, 10), caption = "Base de Datos Cultivo de Aguacate")
Base de Datos Cultivo de Aguacate
id_arbol Latitude Longitude Psychro_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
1 2.393549 -76.71124 22.0 825.1 85.2 0.0 23.9 825.2 0.0 313 312 0.0 25.3 1696 21.3 2.504 23.9 717
2 2.393573 -76.71120 21.4 825.3 84.0 0.0 23.5 825.2 0.0 317 317 0.0 24.8 1696 20.7 2.485 23.5 717
3 2.393541 -76.71113 21.8 825.5 79.6 0.2 24.5 825.5 0.4 338 337 0.5 25.7 1694 20.8 2.518 24.5 717
4 2.393503 -76.71119 22.8 825.4 77.6 0.4 25.9 825.4 0.2 299 299 0.5 28.1 1694 21.7 2.572 25.9 717
5 2.393486 -76.71121 22.6 825.2 76.5 0.0 26.0 825.2 0.0 265 264 0.0 28.0 1696 21.5 2.575 25.9 717
6 2.393441 -76.71126 21.5 825.0 77.7 0.3 24.5 825.0 0.0 265 264 0.3 25.5 1698 20.4 2.522 24.5 717
7 2.393311 -76.71128 22.2 825.0 76.5 0.0 25.5 825.0 0.0 210 209 0.0 27.3 1698 21.1 2.558 25.5 717
8 2.393382 -76.71122 22.6 825.1 77.7 0.0 25.7 825.0 0.0 304 303 0.0 27.5 1698 21.5 2.566 25.6 717
9 2.393425 -76.71120 23.0 825.2 78.3 0.0 26.0 825.2 0.0 284 284 0.0 28.4 1696 22.0 2.579 26.0 717
10 2.393467 -76.71116 22.7 825.4 77.8 0.3 25.9 825.4 0.3 315 315 0.4 28.1 1694 21.7 2.572 25.9 717 
leaflet() %>% addTiles() %>% addCircleMarkers(lng = geo_datos$Longitude,lat = geo_datos$Latitude,radius = 0.2,color = "brown")

Teniendo en cuenta la concentración de la zona de cultivo, es posible establecer que el mismo se ubica en el Municipio de Timbio.

Análisis geoestadístico

A continuación se procede a crear la variable regionalizada

topo_geo= as.geodata(geo_datos,coords.col = 3:2, data.col = 8)
plot(topo_geo)

Semivariograma

Con el propósito de corroborar si existe o no autocorrelación espacial, se efectua el semivariograma para la variable objeto de analisis y con ello se efectua el gráfico correspondiente

summary(dist(geo_datos[,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(topo_geo, option = "bin", uvec = seq(0, 9.178e-04, 8.712e-05))
## variog: computing omnidirectional variogram
plot(variograma)

variograma_mc=variog.mc.env(topo_geo,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(variograma_mc)

Al observar el gráfico, se nota que la semivarianza aumenta rápidamente al inicio, lo que sugiere una alta variabilidad entre puntos muy cercanos. Después de este incremento inicial, la semivarianza se estabiliza, alcanzando una meseta conocida como “sill”. Esta meseta indica el nivel máximo de variabilidad presente en el conjunto de datos. El punto en el que la semivarianza se estabiliza se denomina “rango”, y más allá de esta distancia, los puntos no presentan correlación espacial significativa. En algunos variogramas, puede observarse una discontinuidad en el origen llamada “nugget”, que representa variabilidad a una escala menor que la mínima distancia de muestreo o errores de medición.

Este variograma sugiere que existe una rápida pérdida de correlación espacial a distancias cortas, con la semivarianza alcanzando rápidamente un nivel estable. Esto implica que la mayoría de la variabilidad espacial se encuentra en distancias muy cortas y que, a partir de cierta distancia, la variabilidad ya no depende de la distancia entre puntos. Comprender este patrón es esencial para modelar estructuras espaciales y realizar interpolaciones precisas, como en el método de kriging. Ajustar un modelo teórico de variograma a estos datos permitiría obtener parámetros precisos para técnicas de análisis espacial avanzadas.

Ajuste del modelo teórico

Posteriormente, se efectúa el ajuste del modelo teórico considerando tres opciones: el modelo Exponencial, Gaussiano y Esférico.

ini.vals = expand.grid(seq(2.0, 3.0,l=10), seq(4e-04,6e-04,l=10))
plot(variograma)

modelo_exponencial=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 "3"     "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 101067.987408248
modelo_gaussiano=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 "3"     "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 110561.009962844
modelo_esfer <- 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 "3"     "0"   "0"   "0.5"
## status        "est"   "est" "fix" "fix"
## loss value: 12131.5152418142
plot(variograma)

# Calcular el variograma
variograma <- variog(topo_geo, option = "bin", uvec = seq(0, 9.178e-04, 8.712e-05))
## variog: computing omnidirectional variogram
# Ajustar el modelo exponencial
ini.vals <- expand.grid(seq(2.0, 3.0, l = 10), seq(4e-04, 8e-04, l = 10))
modelo_expon <- 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 "3"     "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 101067.987408248
# Ajustar el modelo gaussiano
modelo_gaussi <- 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 "3"     "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 110561.009962844
# Ajustar el modelo esférico
modelo_esfer <- 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 "3"     "0"   "0"   "0.5"
## status        "est"   "est" "fix" "fix"
## loss value: 12131.5152418142
# Graficar el variograma
plot(variograma)

# Agregar líneas de los modelos al gráfico
lines(modelo_expon, col = "green")
lines(modelo_gaussi, col = "red")
lines(modelo_esfer, col = "purple")

El modelo exponencial generalmente muestra un aumento rápido de la semivarianza a distancias cortas, seguido por un aplanamiento gradual. En el gráfico, la línea azul oscura debería reflejar este comportamiento.

El modelo gaussiano tiende a aumentar más suavemente al principio y luego se aplana. La línea naranja en el gráfico representa este ajuste.

El modelo esférico aumenta linealmente al principio y luego se aplana a un valor constante (sill). La línea púrpura representa este ajuste.

El gráfico y los resultados indican que el modelo esférico es el que mejor describe la variabilidad espacial de los datos, ya que minimiza la función de pérdida más eficazmente que los modelos exponencial y gaussiano. Este modelo refleja que la variabilidad espacial se estabiliza rápidamente a una cierta distancia, después de la cual los puntos ya no están correlacionados espacialmente. Este conocimiento es crucial para técnicas de interpolación como el kriging, donde un buen ajuste del variograma es fundamental para obtener predicciones precisas.

modelo_expon
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: exponential
## parameter estimates:
##   tausq sigmasq     phi 
##  0.0538  3.0539  0.0001 
## Practical Range with cor=0.05 for asymptotic range: 0.0003535736
## 
## variofit: minimised weighted sum of squares = 3730.522
modelo_gaussi
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: gaussian
## parameter estimates:
##   tausq sigmasq     phi 
##  1.2802  1.6990  0.0000 
## Practical Range with cor=0.05 for asymptotic range: 0.0001159668
## 
## variofit: minimised weighted sum of squares = 13512.47
modelo_esfer
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: spherical
## fixed value for tausq =  0 
## parameter estimates:
## sigmasq     phi 
##  3.0644  0.0002 
## Practical Range with cor=0.05 for asymptotic range: 0.000248055
## 
## variofit: minimised weighted sum of squares = 5586.095

Entre los tres modelos evaluados, el modelo esférico presenta el valor de pérdida más bajo (9.83), lo que sugiere que es el modelo que mejor se ajusta a los datos, seguido por el modelo exponencial (44.92) y finalmente el modelo gaussiano (56.16). Los parámetros estimados son similares entre los modelos, con una varianza

##Interpolación (predicción espacial)

c(min(geo_datos[,3]),
max(geo_datos[,3]),
min(geo_datos[,2]),
max(geo_datos[,2]))
## [1] -76.711799 -76.710215   2.392101   2.393634
geodatos_grid=expand.grid( lon=seq(-76.71022,-76.7118,l=100),lat=seq(2.392101 ,2.393634 ,l=100))
plot(geodatos_grid)
points(geo_datos[,3:2],col="red")

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

par(mfrow=c(1,2))

contour(geodatos_ko,main="kriging StDv Predict",val=sqrt(geodatos_ko$krige.var),  drawlabels=TRUE)
image(geodatos_ko, main="kriging StDv Predicted",val=sqrt(geodatos_ko$krige.var), xlab="East", ylab="North")

library(rasterVis)
## Loading required package: lattice
library(lattice)
# Cargar el paquete 'raster'
library(raster)

# Tus datos
pred <- cbind(geodatos_grid, geodatos_ko$predict)

# Crear un objeto raster desde los datos XYZ
temp_predict=rasterFromXYZ(cbind(geodatos_grid,geodatos_ko$predict))

# Graficar el objeto raster
plot(temp_predict)


```r
levelplot(temp_predict)

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

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

library(geoR)
valida=xvalid(geodata = topo_geo, model = modelo_gaussi)
## 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] 1.442553

Con los resultados anteriores, es posible determinar que el modelo presenta un error de predicción con validación cruzada de 1.44

Conclusiones

El analisis presentado incluye los resultados de tres modelos de covarianza diferentes: exponencial, gaussiano y esférico. Cada modelo tiene sus propias características y se ajusta a los datos de manera distinta. Los parámetros y las líneas de ajuste de estos modelos se utilizan para describir la estructura espacial de los datos.