Geo_Aguacate
2024-05-20
Elaborado por: Juan Carlos Villalba Acevedo
Presentado a: Prof. David Arango Londoño
1. Visualización de Datos
head(geo_datos,10)## # A tibble: 10 × 21
## id_arbol Latitude Longitude FORMATTED_DATE_TIME Psychro_Wet_Bulb_Tem…¹
## <dbl> <dbl> <dbl> <chr> <dbl>
## 1 448 2.39 -76.7 01/10/2020 10:11:12 a, m, 20.3
## 2 457 2.39 -76.7 01/10/2020 10:11:12 a, m, 21.4
## 3 444 2.39 -76.7 01/10/2020 10:11:12 a, m, 21.6
## 4 449 2.39 -76.7 01/10/2020 10:11:12 a, m, 20.2
## 5 470 2.39 -76.7 01/10/2020 10:11:12 a, m, 22.5
## 6 471 2.39 -76.7 01/10/2020 10:11:12 a, m, 22.3
## 7 443 2.39 -76.7 01/10/2020 10:11:12 a, m, 20.4
## 8 447 2.39 -76.7 01/10/2020 10:11:12 a, m, 20.7
## 9 445 2.39 -76.7 01/10/2020 10:11:12 a, m, 21.3
## 10 458 2.39 -76.7 01/10/2020 10:11:12 a, m, 22.3
## # ℹ 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>leaflet() %>% addTiles() %>% addCircleMarkers(lng = geo_datos$Longitude,lat = geo_datos$Latitude,radius = 0.2,color = "black")Los datos se encuentran localizados en una finca del municipio de Timbío (Cauca).
2. Geoestadística (Variable de Temperatura)
A continuación, se procede a crear la variable regionalizada:
geo_temp = as.geodata(geo_datos,coords.col = 3:2,data.col = 9)
plot(geo_temp)
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-03Se genera el semivariograma teniendo en cuenta el valor mínimo y el valor correspondiente al tercer cuartil:
variograma = variog(geo_temp, option = "bin", uvec = seq(0,0.0010000,0.00002000))## variog: computing omnidirectional variogramplot(variograma,pch=16)
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 envelopslines(variograma_mc)
Se observa en el semivariograma una autocorrelación espacial fuerte en los valores iniciales, sin embargo, al aumentar la distancia, los valores de temperatura se acercan a un valor constante, lo que indica que la variable regionalizada a partir de la distancia 2e-04 deja de ser espacialmente correlacionada.
3. Ajuste del Modelo Teórico
A continuación, se procede a ajustar un modelo exponencial, gausiano y esférico al semivariograma de la variable regionalizada temperatura:
ini.vals = expand.grid(seq(2.5,3.5,l=10), seq(2e-04,6e-04,l=10))
plot(variograma)
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 "3.39" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 6972.37513629932model_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 "3.28" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 16787.5617532935model_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 "3.17" "0" "0" "0.5"
## status "est" "est" "fix" "fix"
## loss value: 7876.82779918255lines(model_mco_exp,col="blue")
lines(model_mco_gaus,col="red")
lines(model_mco_spe,col="purple")
Se observa que el modelo exponecial (de color azul) se ajusta más a los puntos del semivariograma muestral. Esto nos indica que la correlacion espacial disminuye de manera gradual, la curva es más suave que la de los demás modelos y se acerca a la meseta de manera asintótica para aplastarse complementamente a partir de ese punto.
4. Predicción Espacial
Se procede a calcular los valores máximos y mínimos para definir el área de interés en el cual se realizará la predición espacial:
summary(geo_datos$Latitude)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.392 2.393 2.393 2.393 2.393 2.394summary(geo_datos$Longitude)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -76.71 -76.71 -76.71 -76.71 -76.71 -76.71geodatos_grid=expand.grid( lon=seq(-76.712,-76.710,l=100),lat=seq(2.3920,2.3937,l=100))
plot(geodatos_grid)
points(geo_datos[,3:2],col="red")
Para la predicción espacial se usará el modelo espacial Kriging:
geodatos_ko=krige.conv(geo_temp, loc=geodatos_grid,
krige= krige.control(nugget=0,trend.d="cte",
trend.l="cte",cov.pars=c(sigmasq=12.7080, phi=0.0007 )))## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhoodpar(mfrow=c(1,2))
image(geodatos_ko, main="kriging Predict", xlab="East", ylab="North")
contour(geodatos_ko, main="kriging Predict", add=TRUE, drawlabels=TRUE)
image(geodatos_ko, main="kriging StDv Predicted",val=sqrt(geodatos_ko$krige.var), xlab="East", ylab="North")
contour(geodatos_ko,main="kriging StDv Predict",val=sqrt(geodatos_ko$krige.var), add=TRUE, drawlabels=TRUE)
A continuación, se procede a transformar la imagen obtenida del modelo
Kriging en forma tipo raster:
pred=cbind(geodatos_grid,geodatos_ko$predict)
temp_predict=rasterFromXYZ(cbind(geodatos_grid,geodatos_ko$predict))
plot(temp_predict)
levelplot(temp_predict, par.settings = BuRdTheme)
temp_error=rasterFromXYZ(cbind(geodatos_grid,sqrt(geodatos_ko$krige.var)))
levelplot(temp_error,par.settings =BuRdTheme)
En el siguiente mapa se despliega la imagen raster producto de la predicción del modelo:
crs(temp_predict)=crs("+proj=longlat +datum=WGS84")
leaflet() %>% addTiles() %>% addRasterImage(temp_predict,opacity = 0.6)La siguiente validación cruzada se realiza para evaluar la capacidad predictiva del modelo de interpolación Kriging:
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-validationFinalmente, se calcula la métrica (Mean Absolute Error) utilizada para evaluar la precisión del modelo predictivo:
MAE=mean(abs(valida$error))
MAE## [1] 0.7913021Si el valor del MAE es bajo nos indica que las prediciones estan cerca de los valores observados, si es muy alto indica que la predicciones estan lejos. En este caso, el valor de 0.79 nos indica que el modelo es aceptable debido a que las predicciones no se alejan tanto de los valores observados, en otras palabras, el modelo es preciso.