Geoestadística

Planteamiento del Problema

Se cuenta con información de producción de aguacate para unos árboles de una finca en el departamento del Cauca, se cuenta con variables climáticas complementarias.

Para las mediciones realizadas el 01/10/2020 (534 registros) se desea realizar el análisis geoestadístico que permita interpolar la variable Temperatura.

#Carga de librerías
library(readxl)
library(summarytools)
library(dplyr)
library(lubridate)
library(leaflet)
library(geoR)
library(RColorBrewer)

Análisis exploratorio.

#Cargue de datos
Datos_Aguacate <- read_excel("D:/MaestriaCienciaDatos/S2_AnalisisInfoGeografica/M2U1_Geoestadistica/CasoAccion/Datos_Completos_Aguacate.xlsx")

#Adecuación de fecha
Datos_Aguacate <- Datos_Aguacate %>%
  mutate(
    # Eliminar las comas y el "a m" / "p m"
    fecha_limpia = gsub(",", "", FORMATTED_DATE_TIME),
    fecha_limpia = gsub(" a m", " AM", fecha_limpia),
    fecha_limpia = gsub(" p m", " PM", fecha_limpia),

    # Convertir a fecha POSIXct
    fecha = dmy_hms(fecha_limpia)
  )

# Selección de fecha de interés 01 de octubre de 2020
df <- Datos_Aguacate %>%
  filter(as.Date(fecha) == as.Date("2020-10-01"))

df %>%
  select(Relative_Humidity, Crosswind, Temperature, Barometric_Pressure, Headwind, Wind_Speed, Altitude) %>%
  dfSummary() %>%
  print(method = "render")

Data Frame Summary

df

Dimensions: 534 x 7
Duplicates: 0

No	Variable	Stats / Values	Freqs (% of Valid)	Graph	Valid	Missing
1	Relative_Humidity [numeric]	Mean (sd) : 71.2 (5.3) min ≤ med ≤ max: 59.5 ≤ 71.2 ≤ 91.6 IQR (CV) : 8 (0.1)	200 distinct values		534 (100.0%)	0 (0.0%)
2	Crosswind [numeric]	Mean (sd) : 0.2 (0.2) min ≤ med ≤ max: 0 ≤ 0.1 ≤ 1.5 IQR (CV) : 0.4 (1.2)	13 distinct values		534 (100.0%)	0 (0.0%)
3	Temperature [numeric]	Mean (sd) : 25.8 (1.8) min ≤ med ≤ max: 22.2 ≤ 25.8 ≤ 29.7 IQR (CV) : 2.7 (0.1)	75 distinct values		534 (100.0%)	0 (0.0%)
4	Barometric_Pressure [numeric]	Mean (sd) : 825.5 (1.1) min ≤ med ≤ max: 822.4 ≤ 825.7 ≤ 827.4 IQR (CV) : 1.8 (0)	31 distinct values		534 (100.0%)	0 (0.0%)
5	Headwind [numeric]	Mean (sd) : 0.2 (0.3) min ≤ med ≤ max: -0.7 ≤ 0.1 ≤ 1.3 IQR (CV) : 0.4 (1.6)	20 distinct values		534 (100.0%)	0 (0.0%)
6	Wind_Speed [numeric]	Mean (sd) : 0.3 (0.3) min ≤ med ≤ max: 0 ≤ 0.4 ≤ 1.8 IQR (CV) : 0.5 (1)	15 distinct values		534 (100.0%)	0 (0.0%)
7	Altitude [numeric]	Mean (sd) : 1693 (11) min ≤ med ≤ max: 1675 ≤ 1692 ≤ 1724 IQR (CV) : 15.8 (0)	47 distinct values		534 (100.0%)	0 (0.0%)

Generated by summarytools 1.1.4 (R version 4.5.1)
2025-11-20

Distribución de la temperatura

hist(df$Temperature,
     main = "Histograma de Temperatura",
     xlab = "Temperatura (°C)",
     col = "lightblue",
     breaks = 20)

boxplot(df$Temperature,
        main = "Boxplot de Temperatura",
        ylab = "Temperatura (°C)",
        col = "orange")

Localización de los árboles de aguacate

# Crear paleta continua (azul → rojo)
pal <- colorNumeric(
  palette = "YlOrRd",
  domain = df$Temperature
)

leaflet(df) %>% addTiles() %>% addCircleMarkers(
  lng = df$Longitude,
  lat = df$Latitude, 
  radius = 3, 
  fillColor = ~pal(Temperature),
  fillOpacity = 0.8,
  color = "black",
  weight = 0.3)

Se observa la distribución de los árboles de aguacate con información de temperatura. La finca se encuentra ubicada en el municipio de Timbío en el departamento del Cauca y se observan mayores valores de temperatura (rojos) en la región suroccidental y en el extremo oriental, presentando una variación suave y continua en franjas diagonales de occidente a oriente. Lo cual sugiere un patrón de variación no aleatorio.

Análisis de autocorrelación espacial con el semivariograma.

# Creación de la variable regionalizada
geodatos=as.geodata(df, coords.col=3:2,data.col=9)
plot(geodatos)

Se observa una diferencación en la temperatura de la finca principalmente para los valores más bajos y altos identificados con circulos azules y equis rojas respectivamente.

Resumen estadístico de las distancias entre observaciones

#Estadísticas de matrix de distancias
# Referencia para establecer los intervalos de distancia del semivariograma
summary(dist(geodatos$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-03

# Semivariograma observado
variograma=variog(geodatos,option = "bin",uvec=seq(0,0.0009178,0.00009178))

## variog: computing omnidirectional variogram

plot(variograma,pch=16)

# Semivariogramas aleatorios que definen el intervalo esperado si no hay correlación espacial
variograma_mc=variog.mc.env(geodatos,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

lines(variograma_mc)

Se observa que a medida que aumenta la distancia entre los árboles el valor de la temperatura tiende a presentar una variación. Se observa como el semivariograma observado se diferencia del semivariograma no correlacionado representado por las dos bandas continuas en la gráfica. Por lo cual se puede afirmar que existe una dependencia espacial de la temperatura.

Se observan los siguientes parámetros:

• Alcance: Entre 5e-04 y 7e-04
• Meseta: Cercana a 3
• Efecto pepita: Cercano a 1

Ajuste del mejor modelo teórico.

ini.vals = expand.grid(seq(3,6,l=10), seq(5e-04,7e-04,l=10))

plot(variograma)

# Modelo exponencial
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 "4.67"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 32512.4651218213

# Modelo gaussiano
model_mco_gaus=variofit(variograma, ini=ini.vals, cov.model="gaussian", wei="npair", min="optim",nugget = 1)

## 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 "3"     "0"   "1"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 38189.2213679867

# Modelo esférico
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 "3.33"  "0"   "0"   "0.5"
## status        "est"   "est" "fix" "fix"
## loss value: 15721.417256591

lines(model_mco_exp,col="blue")
lines(model_mco_gaus,col="red")
lines(model_mco_spe,col="purple")

model_mco_exp

## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: exponential
## parameter estimates:
##   tausq sigmasq     phi 
##  0.0000  4.6667  0.0005 
## Practical Range with cor=0.05 for asymptotic range: 0.001471112
## 
## variofit: minimised weighted sum of squares = 32512.47

model_mco_gaus

## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: gaussian
## parameter estimates:
##   tausq sigmasq     phi 
##  0.9408  2.1366  0.0001 
## Practical Range with cor=0.05 for asymptotic range: 0.0002333715
## 
## variofit: minimised weighted sum of squares = 6205.156

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 
##  3.0469  0.0002 
## Practical Range with cor=0.05 for asymptotic range: 0.0001862422
## 
## variofit: minimised weighted sum of squares = 8373.414

Gráficamente y a partir de WSS se observa que el que más se ajusta es el modelo gaussiano.

Predicción espacial con la metodología de kriging.

Definición de la grilla de interés

#max(df[,2])
##Predicción Espacial Kriging
geodatos_grid=expand.grid( lon=seq(-76.7119,-76.7101,l=100),lat=seq(2.3920 ,2.3937 ,l=100))
plot(geodatos_grid)
points(df[,3:2],col="yellow")

geodatos_ko=krige.conv(geodatos, loc=geodatos_grid,
      krige= krige.control(nugget=0.94,trend.d="cte", 
      trend.l="cte",cov.pars=c(sigmasq=2.1366, phi=0.0001 )))

## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood

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

Incertidumbre del Kriging

image(geodatos_ko, main="Prediccion Kriging StDv",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)

Este mapa de incertidumbre del modelo de kriging indica que la zona de cubrimiento de los árboles presenta una predicción del kriging muy confiable al tener una desviación estandar pequeña. Mientras que en las zonas externas se observa una desviación estandar alta esperada debido a la escasez de datos.

Validación cruzada

valida=xvalid(geodata = geodatos,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.8307022

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

Conclusiones

• La finca localizada en el municipio de Timbío (Cauca), presenta valores de temperatura que varían entre 22 y 30 grados Celsius, sobre los cuales no se identifican valores atípicos. Adicionalmente el mapa evidencia una variación térmica no aleatoria dentro de la finca, con temperaturas mayores hacia el suroccidente y extremo oriental, que disminuyen hacia el nororiente.

• A partir del semivariograma, se observa una diferencación en la temperatura de la finca principalmente para los valores más bajos y altos identificados con circulos azules y equis rojas respectivamente. Estos se encuentran predominantemente concentrados y a medida que aumenta la distancia se observa una variación en el valor.

• Tras evaluar los modelos exponencial, gaussiano y esférico, se identifica que el modelo que mejor se ajusta a los datos es el gaussiano con un valor de WSS (suma ponderada de cuadrados) de: 6205.15, en comparación con el valor de 32512 obtenido para el modelo exponencial y 8373 obtenido para el modelo esférico.

• Tras obtener la predicción por kriging se genera el mapa de incertidumbre, el cual asegura una predicción confiable al presentar una desviación estándar pequeña en la zona de interés.

• Se culmina el proceso con la validación cruzada, que indica a partir de el MAE (Error absoluto medio) que en promedio, el modelo se equivoca en ±0.83 grados Celsius en la predicción de temperatura.