knitr::opts_chunk$set(echo = TRUE)

library(readxl)
## Warning: package 'readxl' was built under R version 4.2.3
library(geoR)
## Warning: package 'geoR' was built under R version 4.2.3
## --------------------------------------------------------------
##  Analysis of Geostatistical Data
##  For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
##  geoR version 1.9-2 (built on 2022-08-09) is now loaded
## --------------------------------------------------------------
library(raster)
## Warning: package 'raster' was built under R version 4.2.3
## Loading required package: sp
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.2.3
library(leaflet)
## Warning: package 'leaflet' was built under R version 4.2.3
require(rasterVis)
## Loading required package: rasterVis
## Warning: package 'rasterVis' was built under R version 4.2.3
## Loading required package: lattice

1. Carga de datos y visualización inicial

Se carga la base de corresponde a mediciones climáticas de arboles de aguacate tomadas en el Cauca.

geo_datos =read_excel("C:/Users/ADMIN/Desktop/Maestria Ciencia de Datos/SEMESTRE III - 2023.2/Analisis de Inf Geografica y Espacial/M2U1 - Exploracion de Datos con Geoestadistica/Datos_Completos_Aguacate.xlsx")
leaflet() %>% addTiles() %>% addCircleMarkers(lng = geo_datos$Longitude,lat = geo_datos$Latitude,radius = 0.2,color = "blue")

2. Geoestadística

El objetivo del análisis estadística evaluar si existe una correlación espacial en las mediciones de temperatura, usando técnicas geoestadísticas en R. Inicialmente se crea el objeto geodata el cual es fundamental para el análisis.

2.1. Crear geodata

El objeto se crea a partir de la base cargada anteriormente indicando las columnas de longitud y latitud de las mediciones y la columna de la variable de interés (Temperatura) que en este caso es la 5.

geo_temp=as.geodata(geo_datos,coords.col = 3:2,data.col = 5)
plot(geo_temp)

En este reporte inicial se ve que las temperaturas se encuentran entre los 19 y 25 °C y que hay ciertos lugares donde se concentran en su mayoria temperaturas altas o bajas.

2.2. Variograma

Primero se hace el análisis de las distancias entre los puntos, para evaluar el rango que tendrá la secuencia en la creación del variograma. Los rangos de distancia estarán entre 0 y 0.001, separados por 0.00002.

options(scipen=999)
summary(dist(geo_datos[,3:2]))
##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## 0.00001712 0.00040512 0.00064078 0.00068267 0.00091776 0.00195913
variograma=variog(geo_temp,option = "bin",uvec=seq(0,0.00091776,0.00001712))
## variog: computing omnidirectional variogram
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 envelops
plot(variograma)
lines(variograma_mc)

El semivariograma indica que posiblemente existe una correlación espacial ya que para distancias más pequeñas la semivarianza es más baja. Para las distancias más pequeñas la semivarianza inicia en 0.7 aproximadamente y se estabiliza en 1.3 a una distancia de 0.00025 aproximandamente.

2.3. Ajuste del modelo teórico

Se ajustan tres modelos teóricos: - Exponencial - Gaussiano - Esférico

ini.vals = expand.grid(seq(1.5,1.5,l=10), seq(0.0003,0.0008,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 "1.5"   "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 6508.77757895642
model_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 "1.5"   "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 11216.7379645166
model_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 "1.5"   "0"   "0"   "0.5"
## status        "est"   "est" "fix" "fix"
## loss value: 3695.31528075193
lines(model_mco_exp,col="yellow")
lines(model_mco_gaus,col="blue")
lines(model_mco_spe,col="red")

model_mco_exp
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: exponential
## parameter estimates:
##   tausq sigmasq     phi 
##  0.0011  1.3126  0.0001 
## Practical Range with cor=0.05 for asymptotic range: 0.0002325445
## 
## variofit: minimised weighted sum of squares = 298.3581
model_mco_gaus
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: gaussian
## parameter estimates:
##   tausq sigmasq     phi 
##  0.0387  1.2706  0.0001 
## Practical Range with cor=0.05 for asymptotic range: 0.0001163598
## 
## variofit: minimised weighted sum of squares = 438.3308
model_mco_spe
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: spherical
## fixed value for tausq =  0 
## parameter estimates:
## sigmasq     phi 
##  1.3118  0.0002 
## Practical Range with cor=0.05 for asymptotic range: 0.0001643328
## 
## variofit: minimised weighted sum of squares = 389.7993

Según el gráfico y la suma de cuadrados el modelo que mejor se ajusta es el exponencial, con un valor de sigmasq = 1.3126 y phi = 0.0001 y un minimised weighted sum of squares = 298.3581

2.4. Interpolación (Predicción espacial)

Para la interpolación se crea la malla con la extensión de las mediciones, la cual sirve como entrada del kriggin con los datos y los parámetros estimados en el modelo exponencial.

c(min(geo_datos[,3]),
max(geo_datos[,3]),
min(geo_datos[,2]),
max(geo_datos[,2]))
## [1] -76.711799 -76.710215   2.392101   2.393634
geodatos_grid=expand.grid( lon=seq(-76.711799,-76.710215,l=100),lat=seq(2.392101 ,2.393634 ,l=100))
plot(geodatos_grid)
points(geo_datos[,3:2],col="yellow")

geodatos_ko=krige.conv(geo_temp, loc=geodatos_grid, krige= krige.control(nugget=0,trend.d="cte", trend.l="cte",cov.pars=c(sigmasq=1.3126, phi=0.0001 )))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood

El resultado de la predicción espacial es el siguiente.

par(mfrow=c(1,2))

contour(geodatos_ko,main="kriging Predict", drawlabels=TRUE)
image(geodatos_ko, main="kriging Predict", xlab="East", ylab="North")

Y la desviación estandar espacial:

par(mfrow=c(1,2))

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

2.5. Transformar imagen a raster

Ahora vamos a transformar la imagen en 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)

valida=xvalid(geodata = geo_temp,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.7010667

El modelo presenta un error de predicción de acuerdo con la validación cruzada de 0.7 grados.