Módulo 2 Unidad 1
Fabián Castro
2024-12-08
Análisis del clima en finca de aguacate
Preparación de los datos
- Se realiza cargue de las librerías que se utilizarán para el análisis de geoestadística de esta actividad.
require(geoR)## Loading required package: 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
## --------------------------------------------------------------require(rasterVis)## Loading required package: rasterVis## Loading required package: latticelibrary(leaflet)
library(ggplot2)
library(raster)## Loading required package: splibrary(readxl)- Se realiza el cargue de la información habiendo depuerado las fechas y dejando solo en el archivo la que corresponde a 01/10/2020 10:11:12 a, m. Posterior a esto se crea un DataFrame de 4 columnas que son las coordenadas (latitud y longitud), la fecha y la temperatura.
Datos_Completos_Aguacate <- read_excel("C:/Users/ACER/Downloads/Datos_Completos_Aguacate.xlsx")
datos = Datos_Completos_Aguacate
datos <- datos[, c("Latitude", "Longitude","FORMATTED_DATE_TIME", "Temperature")]
colnames(datos) <- c("Latitud", "Longitud", "Fecha", "Temperatura")
datos## # A tibble: 534 × 4
## Latitud Longitud Fecha Temperatura
## <dbl> <dbl> <chr> <dbl>
## 1 2.39 -76.7 01/10/2020 10:11:12 a, m, 23.9
## 2 2.39 -76.7 01/10/2020 10:11:12 a, m, 23.5
## 3 2.39 -76.7 01/10/2020 10:11:12 a, m, 24.5
## 4 2.39 -76.7 01/10/2020 10:11:12 a, m, 25.9
## 5 2.39 -76.7 01/10/2020 10:11:12 a, m, 26
## 6 2.39 -76.7 01/10/2020 10:11:12 a, m, 24.5
## 7 2.39 -76.7 01/10/2020 10:11:12 a, m, 25.5
## 8 2.39 -76.7 01/10/2020 10:11:12 a, m, 25.7
## 9 2.39 -76.7 01/10/2020 10:11:12 a, m, 26
## 10 2.39 -76.7 01/10/2020 10:11:12 a, m, 25.9
## # ℹ 524 more rows- Se crean los marcadores circulares de la zona que se desea explorar de acuerdo con las coordenadas cargadas en el archivo que contiene los datos a analizar.
leaflet() %>% addTiles() %>% addCircleMarkers(lng = datos$Longitud, lat = datos$Latitud, radius = 0.2, color = "black")Creación de datos espaciales
- Se crea un objeto geodata y se genera el gráfico respectivo dónde se puede evidenciar la ubicación que tiene una mayor temperatura (rojo) y los que tienen una menor temperatura (azul), se logra identificar que la distribución se acerca a una distribución simétrica con una ligera concentración en el centro principalmente en las temperaturas entre los 24 y 26 grados.
geodata = as.geodata(datos, coords.col = 2:1, data.col = 4)
plot(geodata)
- Se obtienen las distancias frente a todas las ubicaciones espaciales y con esto poder obtener estadísticas descriptivas básicas que son el insumo para calcular el semivariograma
summary(dist(datos[,2:1]))## 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-03Análisis de semivariograma
- Con la información estadísticas de las distancias pasamos a generar el semivariograma con los siguientes parámetros:
variograma = variog(geodata, option="bin", uvec=seq(0, 9.178e-04,0.650e-4))## variog: computing omnidirectional variogramdatos.env=variog.mc.env(geodata,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 envelopsplot(variograma,pch=16,envelope=datos.env)
- En el semivariograma se evidencia que hay una correlación espacial principalmente en las distancias más cercanas, conforme aumentan hay menos correlación.
- El valor de la semivarianza se estabiliza en un valor cercano a 3 que sería la meseta.
Ajuste del modelo de semivariograma
ini.vals = expand.grid(seq(3,5,l=10), seq(4e-04,6e-04,l=10))
model_mco_exp = variofit(variograma, ini=ini.vals, cov.model = "exp")## 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 "4.11" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 24465.8949033656model_mco_gaus = variofit(variograma, ini=ini.vals, cov.model="gaus")## 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.67" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 74717.8585662642model_mco_spe = variofit(variograma, ini= ini.vals, cov.model="sph")## 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.22" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 10102.484600263plot(variograma,pch=16,envelope=datos.env)
lines(model_mco_exp,col="blue")
lines(model_mco_gaus, col="red")
lines(model_mco_spe, col="purple")
print(model_mco_exp)## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: exponential
## parameter estimates:
## tausq sigmasq phi
## 0.0000 3.1031 0.0001
## Practical Range with cor=0.05 for asymptotic range: 0.0002329142
##
## variofit: minimised weighted sum of squares = 5119.566print(model_mco_gaus)## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: gaussian
## parameter estimates:
## tausq sigmasq phi
## 0.1123 2.9569 0.0001
## Practical Range with cor=0.05 for asymptotic range: 0.000233168
##
## variofit: minimised weighted sum of squares = 6960.626print(model_mco_spe)## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: spherical
## parameter estimates:
## tausq sigmasq phi
## 0.0674 3.0113 0.0002
## Practical Range with cor=0.05 for asymptotic range: 0.0002491095
##
## variofit: minimised weighted sum of squares = 6256.725En el resultado de la suma de los cuadrados de los errores nos da que el valor más bajo es el del modelo exponencial con 5119, por lo tanto este es el que se utilizará para realizar la predicción espacial Kriging.
Predicción espacial Kriging
geodatos_grid = expand.grid(lon=seq(-76.71022,-76.71180,l=100),lat=seq(2.392101 ,2.393634 ,l=100))
plot(geodatos_grid)
points(datos[,2:1], col="red")
datos_ko=krige.conv(geodata, loc=geodatos_grid,
krige=krige.control(nugget=0,trend.d="cte", trend.l="cte",cov.pars=c(sigmasq=3.1031,phi=0.0001)))## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhoodpar(mfrow=c(1,2))
image(datos_ko, main="Predicción espacial Kriging", xlab="Latitud", ylab="Longitud")
contour(datos_ko, main="Kriking predict", add=TRUE, drawlabels=TRUE)
pred = cbind(geodatos_grid, datos_ko$predict)
temp_predict= rasterFromXYZ(cbind(geodatos_grid, datos_ko$predict))
plot(temp_predict)
levelplot(temp_predict,par.settings =BuRdTheme)
temp_error=rasterFromXYZ(cbind(geodatos_grid,sqrt(datos_ko$krige.var)))
levelplot(temp_error,par.settings =BuRdTheme)
Validación cruzada
valida = xvalid(geodata = geodata, model = model_mco_spe)## 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.7932443En la validación cruzada el modelo presenta un error de predicción de 0,79 grados, esto significa que el modelo difiere de las observaciones reales aproximandamente en esta medida.