Ejercicio Unidad 2.1
Juan Sebastián Meza Buchely
2024-05-22
Análisis de variables para el cultivo de agucate
En el presente trabajo, se procede a analizar la variable de temperatura en el cultivo de aguacate para una finca en el departamento del Cauca. El análisis del trabajo incluye una valoración desde la estadística espacial para identificar si existe correlación espacial e identificar modelos teóricos con el mejor ajuste a los datos, para posteriormente, realizar una predicción e identificar su nivel de error.
En ese orden de ideas, para realizar el procedimiento es necesario comenzar con una exploración inicial de los datos
Exploración inicial del dataset
En un primer momento, se procece a importar la data con la cual se va a trabajar y luego se creará la geodata necesaria para solucionar el ejercicio.
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: splibrary(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")| 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 |
Como se puede observar, la base de datos contiene información acerca de las coordenadas donde se encuentran los árboles de agucate, así como otras características relacionadas con las condiciones climáticas y fenotipicas del cultivo.
A continuación, se procede a identificar el espacio geográfico donde se ubican los cultivos.
leaflet() %>% addTiles() %>% addCircleMarkers(lng = geo_datos$Longitude,lat = geo_datos$Latitude,radius = 0.2,color = "darkgreen")Tal y como se observa en el anterior gráfico, la zona del cultivo se ubica en el Municipio de Timbio, cerca a la capital del departamento del Cauca, Popayán. Es importante reconocer, que en cuanto a su distribución espacial está tiene una forma regular.
Análisis geoestadístico
Luego de realizar, el análisis inicial, procedemos a crear la variable regionalizada. Este proceso se muestra a continuación:
topo_geo= as.geodata(geo_datos,coords.col = 3:2, data.col = 8)
plot(topo_geo)
En el reporte inicial, se puede observar que la temperatura en la zona de cultivo se encuentra entre los 22 a 30 grados centigrados, presentando una mayor frecuencia en una temperatura de 26 grados. De igual manera, se observa en la primera figura que la distribución de las temperatura en la zona parecieran ser aleatorias, es decir, parecen no tener una autocorrelación espacial, sin embargo dicha hipótesis se debe comprobar estadísticamente.
Semivariograma
Para comprobar la existencia o no de autocorrelación espacial, se procede a calcular el semivariograma para la variable de interés. Para ello, se requiere calcular los valores máximos y mínimos de la diferentes puntos para realizar el calculo del gráfico.
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-03Con los datos presentados anteriormente, se procede a calcular el semivariograma
variograma = variog(topo_geo, option = "bin", uvec = seq(0, 9.178e-04, 8.712e-05))## variog: computing omnidirectional variogramplot(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 envelopslines(variograma_mc)
Teniendo en cuenta los resultados presentados previamente, se puede observar que el semivariograma presenta una forma tipica de datos espaciales autocorrelacionados, es decir, que a distancias cortas los datos de la temepratura se parecen y a distancias larga ya no guardan tanta relación. Igualmente, las bandas del semivariograma permiten establecer que en dichos limites se encuentran los semivariograma donde no se presentan autocorrelación.
En consecuencia, se logra reconocer que los datos de la temperatura son espacialmente correlacionadas. Este resultado permite continuar con el análisis para la modelación y estimación de puntos de temperatura.
Ajuste del modelo teórico
Luego de realizar y analizar el semivariograma, se procede a realizar el ajuste del modelo teórico considerando tres opciones: el modelo Exponencial, Gaussiano y Esférico.
ini.vals = expand.grid(seq(2.5,3.5,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.5" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 50467.5575967037modelo_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.5" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 78351.288113765modelo_esferico=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.17" "0" "0" "0.5"
## status "est" "est" "fix" "fix"
## loss value: 9401.47805910977lines(modelo_exponencial,col="darkblue")
lines(modelo_gaussiano,col="orange")
lines(modelo_esferico,col="purple")
Luego de observa la gráfica, se reconoce qu el modelo exponencial parace tener el mejor ajuste frente a los datos obervados. Sin embargo, para aseguranos de la decisión, se procede a ver los valores del SCE para comprobar tanto visualmente como estadísticamente el mejor ajuste del modelo teórico.
modelo_exponencial## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: exponential
## parameter estimates:
## tausq sigmasq phi
## 0.1834 2.9856 0.0001
## Practical Range with cor=0.05 for asymptotic range: 0.0004057251
##
## variofit: minimised weighted sum of squares = 2903.744modelo_gaussiano## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: gaussian
## parameter estimates:
## tausq sigmasq phi
## 0.0885 2.9640 0.0001
## Practical Range with cor=0.05 for asymptotic range: 0.0002331344
##
## variofit: minimised weighted sum of squares = 6077.238modelo_esferico## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: spherical
## fixed value for tausq = 0
## parameter estimates:
## sigmasq phi
## 3.0630 0.0002
## Practical Range with cor=0.05 for asymptotic range: 0.000248561
##
## variofit: minimised weighted sum of squares = 5584.077Considerando los resultados expuestos tanto en la gráfica como en los valores del SCE se logra reconocer que el modelo teórico con el mejor ajuste es el exponencial donde la suma de los errores al cuadrado fue de 2903.744 y su gráfica es la que mejor se ajusta al semivariograma.
Interpolación (predicción espacial)
Luego de realizar el ajuste del modelo, se procede a realizar la interpolación o predicción espacial mediante el método de Kriging
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.393634geodatos_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")
Ahora procedemos a realizar el modelo de predicción espacial
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 neighbourhoodpar(mfrow=c(1,2))
contour(geodatos_ko,main="kriging Predict", drawlabels=TRUE)
image(geodatos_ko, main="kriging Predict", xlab="East", ylab="North")
Ahora se procede a graficar la desviación estandar espacial
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")
Luego de obtener las imágenes previas, procedemos a transformala en un archivo tipo Raster.
require(rasterVis)## Loading required package: rasterVis## Loading required package: latticepred=cbind(geodatos_grid,geodatos_ko$predict)
temp_predict=rasterFromXYZ(cbind(geodatos_grid,geodatos_ko$predict))
plot(temp_predict)
Luego de realizar la transformación al arcihi raster, procedemos a utilizar la funcion Levelplot para identificar en que zonas se presentan las mayores variaciones en relación con el clima. Este procedimiento se realiza a continuación.
library(lattice)
levelplot(temp_predict)
Tal y como se observa en la anterior gráfica, se puede afirmar que las zonas con las mayores temperaturas en las zonas para el cultivo de aguacate se presentan en las latitudes medias. Mientras que en los extremos se evidencia una mayor dispersión de las temperaturas.
Ahora, se continua el análisis con los errores que se puede generar al momento de realizar los modelos estadisticos.
error_temperatura =rasterFromXYZ(cbind(geodatos_grid, sqrt (geodatos_ko$krige.var)))
levelplot(error_temperatura, par.settings=BuRdTheme)
Nuevamente, al analizar el gráfico anterior, se puede reconocer que los errores del modelo seleccionado se encuentran en un rango de 0.2 a 1.8 centigrados, siendo la zona externa la región con las mayores variaciones al momento del cálculo.
Finalmente, procedemos a realizar la validación cruzada del modelo con el objetivo de identificar el error que se generó al momento del cálculo con el modelo.
valida=xvalid(geodata = topo_geo, model = modelo_gaussiano)## 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-validationMAE=mean(abs(valida$error))
MAE## [1] 0.8160037De acuerdo con el anterior resultado, se puede reconocer que el modelo presenta un error de predicción con la validación cruzada de 0.81 grados.
Conclusiones
Con el anterior análisis, se puede reconocer que a través de la modelación estadística espacial es posible reconocer e interpretar datos que permitan facilitar la toma de decisiones alrededor de un conjunto de datos y atributos sobre un fenómeno en particular. A raíz de ello, en el presente caso, el análisis sobre la temperatura en la región donde se cultuva el agucate permitirá que los agricultores tengan mejor información sobre la distribución de temperaturas y predicción de la misma para realizar el cultivo en las zonas más apropiadas para el proceso agricola.
En resumen, se reconoce que la estadistica espacial permitirá mejorar la toma de decisiones de las personas a través de un análisis por regiones y sectores considerando las variables de interés de los individuos.