Actividad 3
Edier Guzman
2025-06-01
En el siguiente informe vamos a predecir la humedad de cualquier punto de una finca ubicada en el Cauca mediante la técnica Kriging, tenemos coordenadas de cada punto representando cada uno de los árboles de la zona y tenemos varios datos de cada uno de los árboles. El dato que usaremos en la presente actividad es la humedad
arboles <- read.csv("D:/Users/edier88/Desktop/maestria_ciencia_datos/Maestria Ciencia de Datos/Semestre 2/Analisis Datos Geograficos/Actividad 3/Datos_Filtrados_Aguacate.csv", stringsAsFactors = FALSE)
leaflet() %>% addTiles() %>% addCircleMarkers(lng = arboles$Longitude,lat = arboles$Latitude,radius = 0.2,color = "black")Hay buena aglomeración de los árboles en el espacio, están distribuidos en el área de forma uniforme.
humedades=as.geodata(arboles,coords.col = 3:2,data.col = 7)
plot(humedades)
Se ve que la humedad de los sitios donde se ubican los árboles difiere en el espacio, hay grupos con un mismo rango de humedades por lo que podemos seguir con nuestro análisis al no estar aleatorizados sin patrón aparente
summary(dist(humedades$coords))## 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-03Escogemos hasta el cuartil 3. Y cogemos el minimo (1.7e-05)
semivariograma=variog(humedades,option = "bin",uvec=seq(0,9.17e-04,1.71e-05))## variog: computing omnidirectional variogramsemivariograma_mc=variog.mc.env(humedades,obj=semivariograma)## variog.env: generating 99 simulations by permutating data values
## variog.env: computing the empirical variogram for the 99 simulations
## variog.env: computing the envelopsplot(semivariograma)
lines(semivariograma_mc)
Como se ve en el semivariograma, existe una autocorrelación fuerte entre las humedades de los diferentes árboles y su distancia.
Ajuste del Modelo de Semivarianza
vemos que la semivarianza se estabiliza entre 30 y 35 y entre 8e-4 y 1e-3(10e-4) en el eje X
ini.vals = expand.grid(seq(30,35,l=10), seq(8e-04,1e-03,l=10))modelo_exponencial=variofit(semivariograma, 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 "35" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 12486443.8890203modelo_gaussiano=variofit(semivariograma, 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 "35" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 23107249.1697939modelo_esferico=variofit(semivariograma, 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 "33.33" "0" "0" "0.5"
## status "est" "est" "fix" "fix"
## loss value: 2388467.00151397plot(semivariograma)
lines(modelo_exponencial,col="blue")
lines(modelo_gaussiano,col="red")
lines(modelo_esferico,col="brown")
Siendo la línea azul la del modelo exponencial la línea roja la del modelo gaussiano y la línea café la del modelo esférico la que mejor ejemplifica lo sucedido es el modelo exponencial
Gráficamente nos quedamos con el modelo exponencial
plot(humedades$coords)
head(humedades$coords[,1])## [1] -76.71124 -76.71120 -76.71113 -76.71119 -76.71121 -76.71126maxLongitud = max(humedades$coords[,1])
maxLongitud## [1] -76.71022minLongitud = min(humedades$coords[,1])
minLongitud## [1] -76.7118maxLatitud = max(humedades$coords[,2])
maxLatitud## [1] 2.393634minLatitud = min(humedades$coords[,2])
minLatitud## [1] 2.392101Predicción Espacial Kriging
humedades_grid=expand.grid( lon=seq(-76.7118,-76.71022,l=100),lat=seq(2.392101 ,2.393634 ,l=100))
plot(humedades_grid)
points(humedades$coords,col="red")
geodatos_ko=krige.conv(humedades,
loc=humedades_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)
pred=cbind(humedades_grid,geodatos_ko$predict)
humedad_prediccion=rasterFromXYZ(cbind(humedades_grid,geodatos_ko$predict))par(mfrow = c(1, 1))
plot(humedad_prediccion)
levelplot(humedad_prediccion,par.settings =BuRdTheme)
humedad_error=rasterFromXYZ(cbind(humedades_grid,sqrt(geodatos_ko$krige.var)))
levelplot(humedad_error,par.settings =BuRdTheme)
Ejecutando en el archivo “conver_to_shape” la funcion st_convex_hull obtenemos el shapefile que usamos aquí para predecir la humedad de los puntos dentro del polígono de la finca
finca=shapefile("D:/Users/edier88/Desktop/maestria_ciencia_datos/Maestria Ciencia de Datos/Semestre 2/Analisis Datos Geograficos/Actividad 3/finca cauca shapefile/finca_cauca_v2.shp")plot(humedad_prediccion)
plot(finca,add=T)
humedad_prediccion_finca=mask(humedad_prediccion,finca)plot(humedad_prediccion_finca)
plot(finca,add=T)
crs(humedad_prediccion_finca)=crs("+proj=longlat +datum=WGS84")leaflet() %>% addTiles() %>% addRasterImage(humedad_prediccion_finca,opacity = 0.6)valida=xvalid(geodata = humedades,model = modelo_exponencial)## 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] 2.607193