Geoestadistica
Julian Diaz
2023-06-10
Cargue de librerÃas y datos
Mapa de la finca
leaflet() %>% addTiles() %>% addCircleMarkers(lng = aguaca_data$Longitude,lat = aguaca_data$Latitude,radius = 0.2,color = "black")Geoestadistica Seleccionamos la variable altitud como objeto de análisis.
En este grafico podemos observar en azul la altitud más baja de la zona y en rojo las más altas. Amarillo y verde serian las altitudes intermedias.
geo_temp=as.geodata(aguaca_data,coords.col = 3:2,data.col = 16)
plot(geo_temp)
a_Análisis descriptivo b.Variograma
Se construye el semivariograma para entender si existe una correlación espacial. Para el caso de la altitud se encuentra que hay una correlación significativa, dado que la mayor parte de los puntos observados en el grafico (puntos negros) se encuentran fuera de la banda
summary(dist(aguaca_data[,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-03b.Variograma
variograma=variog(geodata =geo_temp ,uvec = seq(0,0.0009,0.00005), option = "bin")## variog: computing omnidirectional variogramdatos_env = variog.mc.env(geo_temp, obj.variog =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, envelop = datos_env)
A partir del comportamiento observado, procederemos a realizar el ajuste del modelo a los datos observados, evaluando 3 modelos: Exponencial, gaussiaano y esferico.
plot(variograma, pch=16, envelop = datos_env)
ini.vals = expand.grid( seq(2,4,l=10), seq(0.0017,0.0172, l = 10))
modelo_mco_exp = variofit(variograma, ini.vals, cov.model = "exponential",wei="npair", min="optim") # Modleo exponenci## 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" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 970754245.337251## Warning in variofit(variograma, ini.vals, cov.model = "exponential", wei =
## "npair", : unreasonable initial value for sigmasq + nugget (too low)## Warning in variofit(variograma, ini.vals, cov.model = "exponential", wei =
## "npair", : unreasonable initial value for phi (too high)modelo_mco_gau = variofit(variograma, ini.vals, cov.model = "gaussian",wei="npair", min="optim",nugget = 0) # Modelo gausianico## 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 "4" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 983675060.324604## Warning in variofit(variograma, ini.vals, cov.model = "gaussian", wei =
## "npair", : unreasonable initial value for sigmasq + nugget (too low)## Warning in variofit(variograma, ini.vals, cov.model = "gaussian", wei =
## "npair", : unreasonable initial value for phi (too high)modelo_mco_shp = variofit(variograma, ini.vals, cov.model = "spherical",fix.nug=TRUE, wei="npair", min="optim") # Modelo Esferico## 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 "4" "0" "0" "0.5"
## status "est" "est" "fix" "fix"
## loss value: 955046276.281372## Warning in variofit(variograma, ini.vals, cov.model = "spherical", fix.nug =
## TRUE, : unreasonable initial value for sigmasq + nugget (too low)## Warning in variofit(variograma, ini.vals, cov.model = "spherical", fix.nug =
## TRUE, : unreasonable initial value for phi (too high)lines(modelo_mco_exp, col = "blue")
lines(modelo_mco_gau, col = "red")
lines(modelo_mco_shp, col = "green")
legend("topright", legend = c("Exponential", "Gaussian", "Spherical"), col = c("blue", "red", "green"), lty = 1)
Selección mejor modelo De acuerdo a la gráfica y al MSE
de cada modelo, se oobserva que el modelo esférico es el que mas se
ajusta, teniendo un comportamiento o tendencia lineal que alcanza un
rango (hacia el final)donde se estabiliza.
print(paste("MSE modelo Exponencial =", as.character(modelo_mco_exp$value)))## [1] "MSE modelo Exponencial = 38669955.0717345"print(paste("MSE modelo Gaussiano =", as.character(modelo_mco_gau$value)))## [1] "MSE modelo Gaussiano = 42596031.9652708"print(paste("MSE modelo Esferico =", as.character(modelo_mco_shp$value)))## [1] "MSE modelo Esferico = 9299965.62692245"Predicción
Habiendo seleccionado el modelo tomaremos los parametros del mismo para realizar la predicción.
max(aguaca_data[,3])## [1] -76.71022Realizamos la predicción espacial y los diversos gráficos
##Predicción Espacial Kriging#
aguacadatos_grid=expand.grid(Longitud = seq(-76.7118,-76.7102, l=100),Latitude = seq(2.392,2.394, l=100))
plot(aguacadatos_grid)
points(aguaca_data[,3:2],col="blue")
Se observa que las areas más oscuras corresponden al area de mayor
altitud de acuerdo a lo evidenciado con el análisis exploratorio
aguacadatos_ko <- krige.conv(geo_temp, loc = aguacadatos_grid,
krige = krige.control(nugget = 0, trend.d = "cte",
trend.l = "cte",
cov.pars = c(sigmasq=4, phi=0.04)))## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhoodimage(aguacadatos_ko, main = "Kriging Predict - Altitude", xlab ="Latitude", ylab ="Longitude")
image(aguacadatos_ko, main="kriging StDv Predicted",val=sqrt(aguacadatos_ko$krige.var), xlab="East", ylab="North")
contour(aguacadatos_ko,main="kriging StDv Predict",val=sqrt(aguacadatos_ko$krige.var), add=TRUE, drawlabels=TRUE)
Transformación imagen en raster donde se observa que las zonas en verde
tienen la mayor altitud y en el segundo mapa corresponde a las areas en
rojo
require(rasterVis)## Loading required package: rasterVis## Warning: package 'rasterVis' was built under R version 4.1.3## Loading required package: lattice## Warning: package 'lattice' was built under R version 4.1.3pred=cbind(aguacadatos_grid,aguacadatos_ko$predict)
BaromPres_predict=rasterFromXYZ(cbind(aguacadatos_grid,aguacadatos_ko$predict))
plot(BaromPres_predict)
levelplot(BaromPres_predict,par.settings =BuRdTheme)
Validación cruzada y estimación del error
valida=xvalid(geodata=geo_temp,model = modelo_mco_shp)## 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] 4.134694El modelo presenta un error de predicción de acuerdo con la validación cruzada de 4.1 metros.
finca=shapefile("C:/Users/Julian/Downloads/Municipios202305_shp/Munpio.shp")
plot(BaromPres_predict)
plot(finca,add=T)