Introducción

El presente trabajo aborda el análisis de datos espaciales mediante la aplicación de técnicas geoestadísticas, utilizando el software R y sus librerías especializadas en análisis geoespacial. La geoestadística es particularmente útil cuando se estudian variables regionalizadas y continuas, representadas como procesos estocásticos. Un ejemplo común de este tipo de análisis es el estudio de variables climatológicas como la temperatura, la presión atmosférica, la velocidad del viento y la humedad, entre otras.

En este caso, se emplean datos georreferenciados del departamento del Cauca, centrados en árboles de producción de aguacate y variables climáticas como la humedad relativa, el tiempo y la altitud. El objetivo principal es realizar una interpolación de la temperatura para estimar su comportamiento en áreas donde no se cuenta con mediciones directas. Para ello, se implementan los siguientes pasos metodológicos:

  1. Evaluación de la autocorrelación espacial mediante el semivariograma.
  2. Identificación del modelo teórico que mejor se ajuste a los datos.
  3. Predicción espacial utilizando la metodología de kriging.

Finalmente, se presentan los resultados obtenidos y las conclusiones del análisis.

Desarrollo del Análisis Geoestadístico

Paquete de librerías.

Se consolida los paquetes de librerias a utilizar para el desarrollo del ejercicio.

library(readxl)
library(dplyr)
library(leaflet)
library(ggplot2)
library(geoR)
library(sp)
library(raster)
library(rasterVis)
library(terra)
library(lattice)
library(latticeExtra)
library(kableExtra)

En la carga de la base de datos original se puede observar que cuenta con 20.271 registros en 21 variables, las cuales incluyen datos de la parcela de aguacates en la zona circundante al municipio del Timbio en el departamento del Cauca. Algunas de las variables que incluye la base de datos están relacionadas con la longitud y latitud, temperatura, la presión de la estación, la humedad relativa, la velocidad del viento, altura, densidad, entre otras.

data = read_excel("C:/Users/gustavo.moreno/OneDrive - PUJ Cali/Desktop/Gustavo 2024/Georef/Ejercicio_Aguacate/Datos_Completos_Aguacate.xlsx")
#visualizacion de parcelas originales
leaflet() %>% addTiles() %>% addCircleMarkers(lng = data$Longitude,lat = data$Latitude,radius = 0.2,color = "blue")

El ejercicio nos solicita realizar un filtro en la fecha de análisis, la cual se estipula para el 1 de Octubre del 2020 a las 10:11am.

La nueva base de datos cuenta con un total de 534 registros y conservando las variables.

#realización del filtro x día
data_filtro = data %>% filter(FORMATTED_DATE_TIME == "01/10/2020  10:11:12 a, m,")
dimen=dim(data_filtro)

Se realiza una visualización de la parcela para el analisis. Adicionalmente, se puede observar la existencia de una distribución regular y uniforme de los árboles en el espacio.

#visualización del filtro
leaflet() %>% addTiles() %>% addCircleMarkers(lng = data_filtro$Longitude,lat = data_filtro$Latitude,radius = 0.2,color = "black")

Analisis variable Temperature

En este punto, el ejercicio requiere la elección de una de las variables de análisis geoespacial. Para ello, se propone en primera instancia, utilizar la temperatura. Por ello, es necesario convertirlo en un geodato o variable regionalizad, e imprimimos.

geo_data_temp = as.geodata(obj = data_filtro, coords.col = 3:2, data.col = 9)
plot(geo_data_temp)

Al graficarlo, se puede observar que no hay una tendencia o patrón de la temperatura en el perímetro. Las temperaturas rondan entre los 22° y 30° grados centigrados con pequeñas distribuciones a lo largo de la zona espacial. Por otra parte, se visualiza una distribución normal de los datos con tendencia promedio a los 26° grados centigrados.

ggplot(data = data_filtro, aes(x = Longitude, y = Latitude, color = Temperature)) +
  geom_point(size = 2) +
  scale_color_viridis_c() +
  labs(title = "Distribución de Temperatura",
       x = "Longitud",
       y = "Latitud",
       color = "Temperatura") +
  theme_minimal() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )

1. Autocorrelación espacial [Semivariograma]

A continuación, ajustamos los datos para realizar el semivariograma de la variable temperatura

#semivariograma temperatura
t=dist(geo_data_temp$coords)
summary(t)
##      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-03
semiv_temp = variog(geodata = geo_data_temp,option = "bin", uvec = seq(0, 9.178e-04, 9.178e-05))
## variog: computing omnidirectional variogram
datos_env = variog.mc.env(geo_data_temp, obj.variog = semiv_temp)
## variog.env: generating 99 simulations by permutating data values
## variog.env: computing the empirical variogram for the 99 simulations
## variog.env: computing the envelops

Al realizar el semivariograma, podemos observar que se cumple la correlación espacial, pues entre menor sea la distancia, más se asemejan sus datos, y en ese mismo sentido es mayor la varianza cuando se alejan los puntos. Dado lo anterior, podriamos elegir la temperatura para su estudio, predicción y evaluación de métricas.

Luego, se realizan los modelos teóricos y su respectiva evaluación y análisis.

plot(semiv_temp, pch=16, envelope = datos_env)

ini.vals_temp = expand.grid(seq(2,5,l=10), seq(6e-04, 9e-04, l=10))
model_mco_temp_exp =variofit(semiv_temp, ini=ini.vals_temp, cov.model = "exponential", fix.nug=TRUE, 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 "5"     "0"   "0"   "0.5"
## status        "est"   "est" "fix" "fix"
## loss value: 39269.1275270196
model_mco_temp_gaus =variofit(semiv_temp, ini=ini.vals_temp, cov.model = "gaussian", wei="npair", min="optim")
## 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.67"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 140679.757653415
model_mco_temp_spe =variofit(semiv_temp, ini=ini.vals_temp, cov.model = "spherical", 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 "3.33"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 22729.5850626672
model_mco_temp_mat = variofit(semiv_temp, ini = ini.vals_temp, cov.model = "matern", wei = "npair", min = "optim", nugget = 1, kappa = 5)
## variofit: covariance model used is matern 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi   tausq kappa
## initial.value "5"     "0"   "1"   "5"  
## status        "est"   "est" "est" "fix"
## loss value: 334616.50708926
model_mco_temp_circ = variofit(semiv_temp, ini = ini.vals_temp, cov.model = "circular", wei = "npair", min = "optim")
## variofit: covariance model used is circular 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi   tausq kappa
## initial.value "3.33"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 33372.7189337798
model_mco_temp_cub = variofit(semiv_temp, ini = ini.vals_temp, cov.model = "cubic", wei = "npair", min = "optim")
## variofit: covariance model used is cubic 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi   tausq kappa
## initial.value "3.33"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 27703.1643795236

Elección del mejor modelo teórico

Los modelos “Gaussiano” y “Esférico”, son los que menor error cuadrático medio tienen con 6033.138 y 6686.261, respectivamente. Por su parte los que mayor error presentan en este caso son el “Exponencial” y “Matern”.

Con lo anterior, procederemos a hacer una predicción espacial para esta variable utilizando el método Kriging.

plot(semiv_temp, pch=16, envelope = datos_env)
lines(model_mco_temp_exp, col="blue")
lines(model_mco_temp_gaus, col="red")
lines(model_mco_temp_spe, col="orange")
lines(model_mco_temp_mat, col = "purple")
lines(model_mco_temp_circ, col = "green")
lines(model_mco_temp_cub, col = "black")

legend("bottomright", 
       legend = c("Exponential", "Gaussian", "Spherical", "Matern", "Circular", "Cubic"), 
       col = c("blue", "red", "orange", "purple", "green", "black"), 
       lty = 1, cex = 0.5, bty="n")

error_cuadratico = data.frame(modelo=c("Exponential", "Gaussian", "Spherical", "Matern", "Circular", "Cubic"),
                           Suma_Errores=c(model_mco_temp_exp$value, model_mco_temp_gaus$value,
                                          model_mco_temp_spe$value, model_mco_temp_mat$value,
                                          model_mco_temp_circ$value, model_mco_temp_cub$value))

#Crear tabla para visualizar la suma de errores cuadrados
kbl(error_cuadratico, caption = "<center><b> Error cuadrático medio por Modelo</b></center>", col.names=c("Modelo","Error cuadrático medio"))%>%
  kable_classic(full_width = F)
Error cuadrático medio por Modelo
Modelo Error cuadrático medio
Exponential 39269.128
Gaussian 6033.138
Spherical 6686.261
Matern 334616.507
Circular 6830.621
Cubic 7572.484

A continuación, como parte del ejercicio se propone la predicción de la temperatura en las otras zonas de la parcela. En principio, se realiza el límite de la malla con las coordenadas límite de latitud y longitud.

#Limites de la malla
c(min(data_filtro[,3]),
  max(data_filtro[,3]),
  min(data_filtro[,2]),
  max(data_filtro[,2]))
## [1] -76.711799 -76.710215   2.392101   2.393634

Se crea de la malla o zona a predecir (puntos grises) en relación a la parcela (puntos negros)

geo_datos_grid=expand.grid(lon=seq(-76.711799, -76.710215, l=100),
                           lat=seq(2.392101, 2.393634, l=100))
plot(geo_datos_grid, col="gray")
points(geo_data_temp$coords, col="black")

Predicción Kriggin

Las siguientes imagenes reflejan aquellas zonas donde se realiza la predicción espacial. Se observa algunos patrones y transiciones en los valores predichos, donde un color más claro representa una temperatura más baja y un color más intenso, representa una temperatura más alta teniendo en cuenta el rango evluado entre los 22° y 30° Centigrados.

#Predicción Kriggin
geo_ko_temp = krige.conv(geo_data_temp, loc=geo_datos_grid,
                         krige = krige.control(nugget = 0, trend.d = "cte",
                                               trend.l = "cte", cov.pars = c(sigmaq=4.8289, phi=0.0022)))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
par(mfrow=c(1,2))
contour(geo_ko_temp, main="kriging Predict",  xlab="Longitud", ylab="Latitud")
image(geo_ko_temp, main="Kriggin Predict",  xlab="Longitud", ylab="Latitud")

Kriging StDv Predicted

Los siguientes gráficos permiten identificar áreas donde las predicciones del kriging son menos confiables, las líneas muestran los cambios de la desviación estandar, donde un valor más alto representa mayor incertidumbre en la predicción en determinada área. Las tonalidades más oscuras representan áreas con mayor incertidumbre , mientras que las áreas más claras indican mayor confianza en las predicciones.

par(mfrow=c(1,2))
contour(geo_ko_temp, main="kriging StDv Predict",val=sqrt(geo_ko_temp$krige.var), xlab="Longitud", ylab="Latitud")
image(geo_ko_temp, main="kriging StDv Predicted",val=sqrt(geo_ko_temp$krige.var), xlab="Longitud", ylab="Latitud")

Transformación de la imagen a raster

Ahora, realizamos el ejercicio con la interpretación de imagenes raster. Donde podemos visualizar aquellas áreas con sus respectivas predicciones de temperatura de acuerdo a la escala de color.

pred=cbind(geo_datos_grid, geo_ko_temp$predict)
temp_predict = rasterFromXYZ(cbind(geo_datos_grid, geo_ko_temp$predict))
par(mfrow=c(1,1))
plot(temp_predict)

De igual manera, visualizamos la predicción de la temperatura en otro formato.

levelplot(temp_predict,par.settings =BuRdTheme)

Por último, la predicción del error. Es decir, zonas donde hay alta probabilidad o no, de asertar los valores predichos. Observamos que hacía el centro de la parcela hay menor riesgo de fallo, mientras que a los extremos donde no tenemos referencia, hay una mayor probabilidad de error.

temp_error=rasterFromXYZ(cbind(geo_datos_grid,sqrt(geo_ko_temp$krige.var)))
levelplot(temp_error,par.settings =BuRdTheme)

Validación Cruzada

A continuación, se realiza el método de validación cruzada para evaluar el modelo.

valid = xvalid(geodata = geo_data_temp,model = model_mco_temp_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

El MAE (Mean Absolute Error) es una métrica que mide el error promedio absoluto entre los valores observados y los valores predichos por un modelo. El MAE es útil porque proporciona una medida clara de qué tan lejos, en promedio, están las predicciones del modelo de los valores reales, expresado en las mismas unidades que la variable objetivo. Un valor menor de MAE indica un modelo con mejor desempeño.

El MAE de 1.292256 en la variable temperatura, indica que en promedio, las predicciones están desviadas en aproximadamente 1.29 unidades del valor real. Esto sugiere que el error absoluto promedio representa aproximadamente el 17.2% del rango total , lo que puede considerarse moderado.

Respecto al promedio, el error es cercano al 5% lo que puede ser aceptable en estos contextos.

MAE=mean(abs(valid$error))
MAE
## [1] 1.292256

Analisis variable Station Pressure

Ahora, para efectos del ejercicio se elegirá otra de las variables del dataset y al final poder comparar los resultados de las métricas. En ese sentido, validando el rendimiento de otras variables y su desempeño con el semivariograma, se eligió la variable Station Pressure; la cual indica la presión que ejerce la atmosfera en un punto específico.

Se observan algunos puntos azules son la presión mas baja en la zona inferior izquierda y luego una predominancia de tonalidades verdes y amarillos con presiones medias y medio-altas. No se evidencia un patrón espacial en esa caracteristica, los puntos se encuentran aleatorizados.

#Convertir los dato a variable regionalizada, variable: Presion
geo_data_press = as.geodata(obj = data_filtro, coords.col = 3:2, data.col = 6)
plot(geo_data_press)

ggplot(data = data_filtro, aes(x = Longitude, y = Latitude, color = Station_Pressure)) +
  geom_point(size = 2) +
  scale_color_viridis_c() +
  labs(title = "Distribución de la Presión",
       x = "Longitud",
       y = "Latitud",
       color = "Presión") +
   theme_minimal() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )

Realizamos el mismo proceso que con la variable temperatura, ajuste del variograma.

pr=dist(geo_data_press$coords)
summary(pr)
##      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-03
semiv_press = variog(geodata = geo_data_press,option = "bin", uvec = seq(0, 9.178e-04, 9.178e-05))
## variog: computing omnidirectional variogram
datos_env_press = variog.mc.env(geo_data_press, obj.variog = semiv_press)
## 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(semiv_press, pch=16, envelope = datos_env_press)

Ahora el ajuste de los modelos teoricos:

ini.vals_pres = expand.grid(seq(1, 2, l = 10), seq(5e-04, 8e-04, l = 10))
model_mco_pres_exp = variofit(semiv_press, ini = ini.vals_pres, 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.67"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 587.320366082761
model_mco_pres_gaus = variofit(semiv_press, ini = ini.vals_pres, cov.model = "gaussian", wei = "npair", min = "optim")
## 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.56"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 5360.33442969223
model_mco_pres_spe = variofit(semiv_press, ini = ini.vals_pres, cov.model = "spherical", 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.22"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 1150.39301895435
model_mco_pres_mat = variofit(semiv_press, ini = ini.vals_pres, cov.model = "matern", wei = "npair", min = "optim", nugget = 0.5, kappa = 1)
## variofit: covariance model used is matern 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi   tausq kappa
## initial.value "1.44"  "0"   "0.5" "1"  
## status        "est"   "est" "est" "fix"
## loss value: 127.873593173468
model_mco_pres_circ = variofit(semiv_press, ini = ini.vals_pres, cov.model = "circular", wei = "npair", min = "optim")
## variofit: covariance model used is circular 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi   tausq kappa
## initial.value "1.22"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 1311.95070571413
model_mco_pres_cub = variofit(semiv_press, ini = ini.vals_pres, cov.model = "cubic", wei = "npair", min = "optim")
## variofit: covariance model used is cubic 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi   tausq kappa
## initial.value "1.22"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 2211.8826219878

En este gráfico podemos notar una primera gran diferencia, pues la mayoría de modelos se ajustan casi perfectamente a los puntos del variograma. Lo que indica menor error cuadrático medio y mejor ajuste del modelo.

plot(semiv_press, pch=16, envelope = datos_env_press)
lines(model_mco_pres_exp, col = "blue")
lines(model_mco_pres_gaus, col = "red")
lines(model_mco_pres_spe, col = "orange")
lines(model_mco_pres_mat, col = "purple")
lines(model_mco_pres_circ, col = "green")
lines(model_mco_pres_cub, col = "black")

legend("bottomright", 
       legend = c("Exponential", "Gaussian", "Spherical", "Matern", "Circular", "Cubic"), 
       col = c("blue", "red", "orange", "purple", "green", "black"), 
       lty = 1, cex = 0.5, bty="n")

Para la variable presión de la estación, el modelo que mejor se ajusta resulta siendo el Matern con un valor de error cuadrático medio cercano a los 127. De igual forma, se puede observar en el siguiente cuadro, que los modelos se ajustan mucho mejor que los expuestos para la variable temperatura.

error_cuadratico2 = data.frame(modelo2=c("Exponential", "Gaussian", "Spherical", "Matern", "Circular", "Cubic"),
                           Suma_Errores=c(model_mco_pres_exp$value, model_mco_pres_gaus$value,
                                          model_mco_pres_spe$value, model_mco_pres_mat$value,
                                          model_mco_pres_circ$value, model_mco_pres_cub$value))

#Crear tabla para visualizar la suma de errores cuadrados
kbl(error_cuadratico2, caption = "<center><b> Error cuadrático medio por Modelo</b></center>", col.names=c("Modelo","Error cuadrático medio"))%>%
  kable_classic(full_width = F)
Error cuadrático medio por Modelo
Modelo Error cuadrático medio
Exponential 567.2611
Gaussian 2427.7117
Spherical 1057.5850
Matern 127.8736
Circular 1263.4583
Cubic 1956.1884 
geo_kriggin = krige.conv(geo_data_press, loc=geo_datos_grid,
                         krige = krige.control(nugget = 0, trend.d = "cte",
                                               trend.l = "cte", cov.pars = c(sigmaq=1.6667, phi=0.0005)))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood

Los mapas de predicción muestran cierta similitud a las predicciones de la temperatura en sus diferentes tonalidades.

par(mfrow=c(1,2))
contour(geo_kriggin, main="kriging Predict",  xlab="Longitud", ylab="Latitud")
image(geo_kriggin, main="Kriggin Predict",  xlab="Longitud", ylab="Latitud")

Lo mismo sucede para el patron de la desviación estándar. El gráfico que has generado muestra el desvío estándar (StDv) de las predicciones realizadas mediante kriging, lo cual es una medida de la incertidumbre asociada a esas predicciones

#Kriging StDv Predicted
par(mfrow=c(1,2))
contour(geo_kriggin, main="kriging StDv Predict",val=sqrt(geo_kriggin$krige.var), xlab="Longitud", ylab="Latitud")
image(geo_kriggin, main="kriging StDv Predicted",val=sqrt(geo_kriggin$krige.var), xlab="Longitud", ylab="Latitud")

La imagen raster muestra una predicción de mayor presion de la estación hacia el norte de la parcela, con algunas tendencias mayores en el centro. La media de esta variable estaba entre los valores 826 y 827, por eso la predominancia del color verde y amarillos en el mapa.

#Imagen raster con predicción
pred=cbind(geo_datos_grid, geo_kriggin$predict)
press_predict = rasterFromXYZ(cbind(geo_datos_grid, geo_kriggin$predict))
plot(press_predict)

Se puede observar la misma predominancia en los siguientes gráficos, y en donde valores lejanos a las zonas de medición definida hay un mayor rango de improbabilidad de la precisión de medida.

levelplot(press_predict,par.settings =BuRdTheme)

press_error=rasterFromXYZ(cbind(geo_datos_grid,sqrt(geo_kriggin$krige.var)))
levelplot(press_error,par.settings =BuRdTheme)

Por ultimo, realizamos el proceso de validación cruzada.

valid = xvalid(geodata = geo_data_press,model = model_mco_pres_exp)
## 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

Obtenemos un valor del MAE de 0.42, inclusive mucho menor que el MAE para la variable temperatura (1.29), lo que indica que en promedio las predicciones para está variable están mejor ajustadas por el modelo evaluado. Sin embargo, aún resultan siendo aceptables ambos resultados y se realiza con fines comparativos.

MAE=mean(abs(valid$error))
MAE
## [1] 0.4290692

Resultados y Conclusiones

Como parte inicial del desarrollo del ejercicio, se utilizó la temperatura como variable a predecir. Sin embargo, de forma paralela, se realizó el mismo análisis para otras cuatro variables, identificándose que el modelo teórico se ajustaba mejor a los datos de Station_Pressure. Por esta razón, se continuó el desarrollo con esta variable, logrando un MAE de 0,42, ligeramente mejor que el obtenido para la temperatura. No obstante, de acuerdo con la literatura, ambos valores se encuentran dentro de rangos aceptables para la interpolación o predicción de variables en áreas donde no se dispone de mediciones directas.

1. Selección de modelos teóricos: Para la variable temperatura, los modelos teóricos ajustados muestran diferentes niveles de precisión, con algunos modelos (como el esférico y el cúbico) acercándose más al comportamiento observado. En el caso de la presión (Station_Pressure), el ajuste a los datos es más consistente, sugiriendo que esta variable presenta una estructura espacial más definida y predecible en el área estudiada.

2. Desempeño del modelo: El valor del MAE obtenido para la interpolación de temperatura fue aceptable, pero al cambiar a la variable presión, el modelo mejoró ligeramente, alcanzando un MAE de 0,42. Ambos valores de MAE están dentro de rangos considerados aceptables en la literatura para estudios de interpolación, indicando que los métodos geoestadísticos utilizados son apropiados para la predicción en áreas no muestreadas.

3. Importancia de la autocorrelación espacial: Los semivariogramas confirman que ambas variables tienen una autocorrelación espacial significativa, lo que valida el uso de técnicas como kriging para la interpolación.

Aunque los resultados fueron satisfactorios, los errores residuales podrían reducirse con un mayor número de puntos de muestreo o mediante la incorporación de variables adicionales que expliquen mejor las fluctuaciones locales. La capacidad de estimar la temperatura y la presión en zonas no medidas proporciona herramientas valiosas para planificar actividades como el riego, la fertilización y el manejo del estrés climático.