datos <- read_excel("C:/Users/andre/OneDrive/Escritorio/Maestria en Ciencia de Datos/Primer Semestre/3_Electiva_Analisis_Infor_Geografia y Esp/Modulo 2 Geoestadistica/Act/datos.xlsx")

Ubicacion Arboles de Aguacate en la finca

plot(datos[,2:3],
     main="Ubicacion de Arboles de Agucate en la Finca",
     col="#00CD00")

Temperatura

datos_geo=as.geodata(datos,coords.col = 2:3,data.col = 8)
plot(datos_geo)

distancia=dist(datos_geo$coords[1:10])
distancia
##           1        2        3        4        5        6        7        8
## 2  0.000024                                                               
## 3  0.000008 0.000032                                                      
## 4  0.000046 0.000070 0.000038                                             
## 5  0.000063 0.000087 0.000055 0.000017                                    
## 6  0.000108 0.000132 0.000100 0.000062 0.000045                           
## 7  0.000238 0.000262 0.000230 0.000192 0.000175 0.000130                  
## 8  0.000167 0.000191 0.000159 0.000121 0.000104 0.000059 0.000071         
## 9  0.000124 0.000148 0.000116 0.000078 0.000061 0.000016 0.000114 0.000043
## 10 0.000082 0.000106 0.000074 0.000036 0.000019 0.000026 0.000156 0.000085
##           9
## 2          
## 3          
## 4          
## 5          
## 6          
## 7          
## 8          
## 9          
## 10 0.000042
resumen=summary(distancia)
print(format(resumen,scientific = FALSE))
##         Min.      1st Qu.       Median         Mean      3rd Qu.         Max. 
## "0.00000800" "0.00004500" "0.00008200" "0.00009582" "0.00013000" "0.00026200"
semiva=variog(datos_geo, option="bin", uvec = seq(0, 0.00023000, 0.0000101))
## variog: computing omnidirectional variogram

Semivariograma

plot(semiva,
     main="Semivariograma Finca Aguacate",
     xlab="Distancia",
     ylab="Semivarianza",
     pch=19,
     col="blue",
     cex=1.1)
grid(nx=NULL, ny=NULL, col="black", lty="dotted")

Interpretación:

Tendencia general: Observamos que la semivarianza tiende a aumentar a medida que la distancia entre los puntos de muestreo también aumenta. Esto sugiere que existe una autocorrelación espacial, es decir, los puntos que están más cerca entre sí tienden a tener valores más similares que los puntos que están más alejados. Meseta: Parece que la semivarianza se estabiliza alrededor de un valor constante a partir de cierta distancia. Este valor se le conoce como meseta o sill y representa la varianza total del conjunto de datos. Lo cual indica que a partir de esa distancia, la influencia espacial deja de ser significativa. Rango: La distancia a la que se alcanza la meseta se denomina rango. En este caso, el rango parece estar alrededor de 0.001 unidades de distancia (el valor exacto dependería de la escala del gráfico). El rango nos da una idea de la extensión de la autocorrelación espacial.

datos.env=variog.mc.env(datos_geo,obj=semiva)
## variog.env: generating 99 simulations by permutating data values
## variog.env: computing the empirical variogram for the 99 simulations
## variog.env: computing the envelops
plot(semiva,pch=16,envelope=datos.env,
     xlab="Distancia",
      ylab= "Semivarianza",
       col="blue",
          cex=1.1)
grid(nx=NULL, ny=NULL, col="black", lty="dotted")

Ajuste del Modelo Teorico

ini.vals = expand.grid(seq(2,4,l=10), seq(0.00008,0.000020,l=10))
# Ajustar un modelo exponencial

model_mco_exp=variofit(semiva, 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 "2.89"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 527.510708434048
# Ajustar un modelo gaussiano

model_mco_gaus=variofit(semiva, 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 "2.67"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 737.757194972905
# Ajustar un modelo esferico

model_mco_spe=variofit(semiva, 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 "2.44"  "0"   "0"   "0.5"
## status        "est"   "est" "fix" "fix"
## loss value: 2611.28606179996
print(model_mco_exp,col="blue")
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: exponential
## parameter estimates:
##   tausq sigmasq     phi 
##  0.0000  2.8590  0.0001 
## Practical Range with cor=0.05 for asymptotic range: 0.0002328517
## 
## variofit: minimised weighted sum of squares = 524.4337
print(model_mco_gaus,col="red")
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: gaussian
## parameter estimates:
##   tausq sigmasq     phi 
##  0.0221  2.5464  0.0001 
## Practical Range with cor=0.05 for asymptotic range: 0.0001164077
## 
## variofit: minimised weighted sum of squares = 723.0815
print(model_mco_spe,col="purple")
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: spherical
## fixed value for tausq =  0 
## parameter estimates:
## sigmasq     phi 
##  2.3336  0.0001 
## Practical Range with cor=0.05 for asymptotic range: 7.616976e-05
## 
## variofit: minimised weighted sum of squares = 2597.765

Modelos teoricos ajustado

# Comparar los variogramas empíricos y ajustados
plot(semiva, main="Variograma Empírico y Modelos Ajustados")
lines(model_mco_exp, col="red", lty=20) # Exponencial
lines(model_mco_gaus, col="blue", lty=20) # Gaussiano
lines(model_mco_spe, col="green", lty=20) # Esférico
legend("bottomright", legend=c("Exponencial", "Gaussiano", "Esférico"), col=c("red", "blue", "green"), lty=6:8)

Interpretación:

Ajuste de los modelos: El gráfico muestra cómo los tres modelos teóricos intentan capturar la tendencia del variograma empírico. A simple vista, parece que el modelo esférico (verde) se ajusta mejor a los datos, especialmente a distancias cortas y medias. El modelo exponencial (rojo) también parece un ajuste razonable, mientras que el gaussiano (azul) se desvía un poco más de los datos.

Comportamiento espacial: El variograma empírico muestra que la semivarianza aumenta con la distancia, lo que indica autocorrelación espacial positiva. Esto significa que los puntos cercanos entre sí tienden a tener valores más similares que los puntos más alejados.

Rango: El rango es la distancia a la que el variograma alcanza una meseta, es decir, donde la semivarianza se estabiliza. En este caso, el rango parece estar alrededor de 0.001 unidades de distancia (el valor exacto dependería de la escala del gráfico). El rango nos da una idea de la extensión de la autocorrelación espacial.

Efecto pepita: El efecto pepita es la semivarianza a una distancia de cero. En este gráfico, el efecto pepita parece ser bajo, lo que sugiere poca variabilidad a pequeña escala o un error de medición pequeño.

Interpolacion ( prediccion espacial)

Prediccion Espacial de Kriging

datos_grid=expand.grid(Latitude=seq(2.392101,2.393634,l=100), Longitude=seq(-76.71022,-76.7118,l=100))
plot(datos_grid)
points(datos[,2:3],col="green")

datos_ko=krige.conv(datos_geo, loc=datos_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 neighbourhood
par(mfrow=c(1,2))
image(datos_ko, main="kriging Predict", xlab="North", ylab="East")
contour(datos_ko,main="kriging Predict", add=TRUE, drawlabels=TRUE)

image(datos_ko, main="kriging StDv Predicted",val=sqrt(datos_ko$krige.var), xlab="North", ylab="East")
contour(datos_ko,main="kriging StDv Predict",val=sqrt(datos_ko$krige.var), add=TRUE, drawlabels=TRUE)

Imagen Raster

require(rasterVis)
require(sp)
require(raster)
pred=cbind(datos_grid,datos_ko$predict)
temp_predict=rasterFromXYZ(cbind(datos_grid,datos_ko$predict))
plot(temp_predict)

require(rasterVis)

crs(temp_predict)=crs("+proj=longlat +datum=WGS84")
leaflet() %>% addTiles() %>% addRasterImage(temp_predict,opacity = 1.8)

Validacion cruzada

valida=xvalid(geodata = datos_geo,model = model_mco_gaus)
## 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-validation
MAE=mean(abs(valida$error))
MAE
## [1] 0.8445374

Teniendo como base el resultado anterior, se puede decir que el modelo presenta un error de prediccion de acuerdo con la validacion cruzada de 0.84 grados.