Geoestadística — Interpolación de Temperatura (Kriging)

1 Introducción

En este trabajo realizaremos un análisis geoestadístico sobre el primer periodo (01/10/2020) de las mediciones de los árboles de aguacate en Cauca. El objetivo principal es interpolar la variable temperatura usando variogramas y kriging (si la variable muestra autocorrelación espacial). Si temperatura no presenta autocorrelación, probaremos con otras variables (humedad relativa, altura, velocidad del viento, etc.). El flujo de trabajo sigue los pasos vistos en clase: EDA, creación de geodato, detección de tendencia, variograma experimental, ajuste de modelos teóricos (exponencial, gaussiano, esférico) y predicción espacial por kriging.

1.1 Preparacion de los datos

# --- Cargar librerías necesarias ---
if (!require(readxl)) install.packages("readxl")
library(readxl)
if (!require(geoR)) install.packages("geoR")
library(geoR)
if (!require(ggplot2)) install.packages("ggplot2")
library(ggplot2)


ruta <- "C:/Users/angie/OneDrive/ANGIE COLLAZOS/MSESTRIA CIENCIA DE DATOS/SIG/SIG/aguacate/Datos_Completos_Aguacate.xlsx"

datos <- read_excel(ruta)
datos <- as.data.frame(datos)

names(datos)
head(datos, 5)

library(dplyr)

# ajustamos el formato de fecha
datos$FORMATTED_DATE_TIME <- gsub(",", "", datos$FORMATTED_DATE_TIME)
datos$FORMATTED_DATE_TIME <- trimws(datos$FORMATTED_DATE_TIME)

datos$FECHA <- sub(" .*", "", datos$FORMATTED_DATE_TIME)

# formato fecha
datos$FECHA <- as.Date(datos$FECHA, format = "%d/%m/%Y")
unique(datos$FECHA)

##  [1] "2019-08-21" "2019-09-04" "2019-09-19" "2019-10-01" "2019-10-18"
##  [6] "2019-11-13" "2020-07-16" "2020-07-30" "2020-08-11" "2020-09-30"
## [11] "2020-10-14" "2019-08-14" "2019-08-15" "2019-08-28" "2019-09-12"
## [16] "2019-09-25" "2019-10-09" "2020-07-07" "2020-07-21" "2020-08-06"
## [21] "2020-08-20" "2020-10-01" "2020-07-09" "2020-07-23" "2020-08-03"
## [26] "2020-08-18" "2020-08-26" "2020-09-18" "2020-09-23" "2020-10-05"
## [31] "2019-08-08" "2019-08-27" "2019-09-11" "2019-09-24" "2019-10-08"
## [36] "2019-11-05" "2020-07-14" "2020-07-28" "2020-08-13" "2020-08-27"
## [41] "2020-09-09" "2020-09-22" "2020-10-06" "2020-10-23" "2019-07-30"
## [46] "2019-08-22" "2019-09-03" "2019-09-20" "2019-10-02" "2019-10-16"

# filtramos segun la actividad
datos_2020 <- filter(datos, FECHA == as.Date("2020-10-01"))

# Revisamos cuántos registros quedaron
nrow(datos_2020)

## [1] 534

# Vemos las primeras filas
head(datos_2020, 6)

##   id_arbol Latitude Longitude      FORMATTED_DATE_TIME
## 1        1 2.393549 -76.71124 01/10/2020  10:11:12 a m
## 2        2 2.393573 -76.71120 01/10/2020  10:11:12 a m
## 3        3 2.393541 -76.71113 01/10/2020  10:11:12 a m
## 4        4 2.393503 -76.71119 01/10/2020  10:11:12 a m
## 5        5 2.393486 -76.71121 01/10/2020  10:11:12 a m
## 6        6 2.393441 -76.71126 01/10/2020  10:11:12 a m
##   Psychro_Wet_Bulb_Temperature Station_Pressure Relative_Humidity Crosswind
## 1                         22.0            825.1              85.2       0.0
## 2                         21.4            825.3              84.0       0.0
## 3                         21.8            825.5              79.6       0.2
## 4                         22.8            825.4              77.6       0.4
## 5                         22.6            825.2              76.5       0.0
## 6                         21.5            825.0              77.7       0.3
##   Temperature Barometric_Pressure Headwind Direction_True Direction_Mag
## 1        23.9               825.2      0.0            313           312
## 2        23.5               825.2      0.0            317           317
## 3        24.5               825.5      0.4            338           337
## 4        25.9               825.4      0.2            299           299
## 5        26.0               825.2      0.0            265           264
## 6        24.5               825.0      0.0            265           264
##   Wind_Speed Heat_Stress_Index Altitude Dew_Point Density_Altitude Wind_Chill
## 1        0.0              25.3     1696      21.3            2.504       23.9
## 2        0.0              24.8     1696      20.7            2.485       23.5
## 3        0.5              25.7     1694      20.8            2.518       24.5
## 4        0.5              28.1     1694      21.7            2.572       25.9
## 5        0.0              28.0     1696      21.5            2.575       25.9
## 6        0.3              25.5     1698      20.4            2.522       24.5
##   Estado_Fenologico_Predominante Frutos_Afectados      FECHA
## 1                            717                0 2020-10-01
## 2                            717                0 2020-10-01
## 3                            717                0 2020-10-01
## 4                            717                0 2020-10-01
## 5                            717                0 2020-10-01
## 6                            717                0 2020-10-01

1.2 Analisis Exploratorio

# resumen estadístico de la temperatura
summary(datos_2020$Temperature)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   22.20   24.50   25.80   25.83   27.18   29.70

# Histograma 
hist(datos_2020$Temperature,
     main = "Histograma de la Temperatura (01/10/2020)",
     xlab = "Temperatura (°C)",
     col = "lightblue")

# Boxplot simple
boxplot(datos_2020$Temperature,
        main = "Boxplot de la Temperatura (°C)",
        horizontal = TRUE,
        col = "lightgreen")

# Gráfico de dispersión (coordenadas de los árboles)
plot(datos_2020$Longitude, datos_2020$Latitude,
     main = "Distribución espacial de los árboles (01/10/2020)",
     xlab = "Longitud", ylab = "Latitud",
     pch = 20, col = "darkgreen")

En la primera exploración de los datos del 1 de octubre de 2020, se observaron 534 registros correspondientes a árboles de aguacate georreferenciados. La temperatura varía entre 22.2 °C y 29.7 °C, con un promedio cercano a 25.8 °C. El histograma muestra una distribución aproximadamente simétrica, aunque con una ligera concentración de valores entre 24 °C y 27 °C, lo que sugiere condiciones térmicas relativamente homogéneas en la zona. El boxplot no evidencia valores atípicos extremos, indicando que las mediciones son consistentes.

En el mapa de dispersión se observa una disposición ordenada de los puntos, que corresponde a la ubicación espacial de los árboles muestreados en el predio o zona de estudio. Esta distribución es fundamental para aplicar los métodos de geoestadística, ya que muestra una cobertura regular del área.

1.3 Creación del Geodato

# Creamos el objeto tipo geoR
geodatos_temp <- as.geodata(datos_2020,
                            coords.col = c("Longitude", "Latitude"),
                            data.col = "Temperature")

# Resumen del objeto
summary(geodatos_temp)

## Number of data points: 534 
## 
## Coordinates summary
##     Longitude Latitude
## min -76.71180 2.392101
## max -76.71022 2.393634
## 
## Distance summary
##          min          max 
## 1.711724e-05 1.959127e-03 
## 
## Data summary
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 22.20000 24.50000 25.80000 25.82903 27.17500 29.70000

# Visualización básica del objeto geodato
plot(geodatos_temp)

# Modelo lineal simple de temperatura respecto a coordenadas
tendencia <- lm(geodatos_temp$data ~ geodatos_temp$coords[,1] + geodatos_temp$coords[,2])
summary(tendencia)

## 
## Call:
## lm(formula = geodatos_temp$data ~ geodatos_temp$coords[, 1] + 
##     geodatos_temp$coords[, 2])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0761 -1.1765 -0.0286  1.1375  4.1201 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -59583.4    19009.2  -3.134 0.001817 ** 
## geodatos_temp$coords[, 1]   -799.3      244.1  -3.274 0.001130 ** 
## geodatos_temp$coords[, 2]   -712.4      213.0  -3.345 0.000882 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.7 on 531 degrees of freedom
## Multiple R-squared:  0.08257,    Adjusted R-squared:  0.07911 
## F-statistic: 23.89 on 2 and 531 DF,  p-value: 1.157e-10

Con la función as.geodata() se construyó el objeto geoespacial que combina las coordenadas de los árboles y sus valores de temperatura. En el resumen del objeto se observa que el área de estudio se encuentra entre las coordenadas −76.7118 y −76.7102 de longitud y 2.3921 a 2.3936 de latitud, con 534 mediciones válidas.

Posteriormente se evaluó la posible tendencia espacial de la temperatura mediante un modelo lineal simple, considerando la longitud y la latitud como variables explicativas. Los resultados muestran que ambos coeficientes son estadísticamente significativos (p < 0.01), lo que indica que existe una leve tendencia espacial: la temperatura varía ligeramente con la ubicación.

Sin embargo, el valor de R² ≈ 0.08 sugiere que esta tendencia explica una proporción reducida de la variación total, por lo que la estructura espacial de la variable deberá analizarse con mayor detalle a través del semivariograma.

1.4 Semivariograma experimental

# Calcular distancias (para referencia)
summary(dist(geodatos_temp$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

# Crear el semivariograma experimental
variograma_temp <- variog(geodatos_temp, max.dist = 0.0015)  # el valor puede ajustarse según la escala

## variog: computing omnidirectional variogram

plot(variograma_temp, main = "Semivariograma experimental - Temperatura")

# Crear una envolvente por simulación (para ver significancia de autocorrelación)
variograma_env <- variog.mc.env(geodatos_temp, obj.var = variograma_temp, nsim = 99)

## 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_temp)
lines(variograma_env$u, variograma_env$upper, lty = 2)
lines(variograma_env$u, variograma_env$lower, lty = 2)

A partir del conjunto de datos del 1 de octubre de 2020 se construyó el semivariograma experimental de la temperatura. En el gráfico se observa que la semivarianza aumenta a medida que lo hace la distancia entre los puntos, alcanzando una meseta aproximadamente a partir de los 0.001 grados de separación.

Este comportamiento indica la presencia de autocorrelación espacial positiva, es decir, las mediciones de temperatura tomadas en árboles cercanos tienden a ser más similares que aquellas más distantes. La forma del variograma sugiere la existencia de una estructura espacial definida, lo cual permite continuar con el ajuste de un modelo teórico de variograma (exponencial, esférico o gaussiano) para posteriormente realizar la interpolación espacial mediante kriging.

1.5 Ajuste de los modelos

# Creamos una grilla de valores iniciales
ini.vals <- expand.grid(seq(8, 12, length = 10),
                        seq(6e-04, 8e-04, length = 10))

# Visualizamos el variograma experimental
plot(variograma_temp)

# Modelo Exponencial
model_exp <- variofit(variograma_temp,
                      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 "8"     "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 331571.360898001

# Modelo Gaussiano
model_gaus <- variofit(variograma_temp,
                       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 "8"     "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 557739.825951394

# Modelo Esférico
model_esf <- variofit(variograma_temp,
                      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 "8"     "0"   "0"   "0.5"
## status        "est"   "est" "fix" "fix"
## loss value: 1942538.50492307

# Agregamos las líneas de los tres modelos al gráfico
lines(model_exp, col = "blue", lwd = 2)
lines(model_gaus, col = "red", lwd = 2)
lines(model_esf, col = "purple", lwd = 2)

# Agregamos una leyenda
legend("bottomright",
       legend = c("Exponencial", "Gaussiano", "Esférico"),
       col = c("blue", "red", "purple"),
       lwd = 2)

Se realizó el ajuste de tres modelos teóricos de semivarianza: exponencial, gaussiano y esférico, utilizando el método de mínimos cuadrados ponderados con valores iniciales definidos en una grilla. En el gráfico se observa que el modelo exponencial (línea azul) sigue más de cerca la tendencia del semivariograma experimental, especialmente en las distancias cortas, donde la autocorrelación espacial es más fuerte.

Los resultados numéricos también confirman este comportamiento: el modelo exponencial presentó el menor valor del error del ajuste (loss value = 331 571), mientras que los modelos gaussiano y esférico mostraron errores considerablemente mayores.

Por lo tanto, el modelo exponencial fue seleccionado como el que mejor representa la estructura espacial de la variable temperatura, y se utilizará como base para la interpolación espacial mediante el método de kriging.

1.6 Predicción Kriging

# Revisamos los rangos de coordenadas
range_long <- range(geodatos_temp$coords[,1])
range_lat  <- range(geodatos_temp$coords[,2])
range_long

## [1] -76.71180 -76.71022

range_lat

## [1] 2.392101 2.393634

# Creamos la cuadrícula de predicción (100 x 100 puntos)
geodatos_grid <- expand.grid(
  lon = seq(range_long[1], range_long[2], length = 100),
  lat = seq(range_lat[1], range_lat[2], length = 100)
)

# Visualizamos la malla 
plot(geodatos_grid, pch = ".", main = "Malla de predicción (100x100)")
points(geodatos_temp$coords, col = "red", pch = 20)

# modelo exp
model_exp$cov.pars

## [1] 8.000000202 0.001755716

model_exp$nugget

## [1] 1.123523e-06

geodatos_ko <- krige.conv(
  geodata = geodatos_temp,
  loc = geodatos_grid,
  krige = krige.control(
    nugget = model_exp$nugget,
    trend.d = "cte",
    trend.l = "cte",
    cov.pars = model_exp$cov.pars,   # parámetros del modelo
    cov.model = "exponential"        # modelo elegido
  )
)

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

Para la interpolación se generó una malla regular de predicción de 100×100 puntos que cubre la extensión de las coordenadas observadas (longitud entre −76.71180 y −76.71022, latitud entre 2.392101 y 2.393634). La malla permite estimar la temperatura en una superficie continua; la densidad de puntos muestreados (puntos rojos) es mayor en el interior del dominio, por lo que las predicciones son más confiables en esa zona que en los bordes.

#Visualización
par(mfrow = c(1, 2))

# ️Mapa de predicción
image(geodatos_ko,
      main = "Predicción de Temperatura-K",
      xlab = "Longitud", ylab = "Latitud")
contour(geodatos_ko, add = TRUE, drawlabels = TRUE)

#  Mapa de desviación estándar (incertidumbre)
image(geodatos_ko,
      main = "Desv.estándar de la predicción-K",
      val = sqrt(geodatos_ko$krige.var),
      xlab = "Longitud", ylab = "Latitud")
contour(geodatos_ko, val = sqrt(geodatos_ko$krige.var), add = TRUE, drawlabels = TRUE)

El mapa de predicción revela un patrón espacial no homogéneo: se observan áreas con temperaturas ligeramente más altas y otras más bajas, organizadas en una estructura continua que sugiere gradientes locales. El modelo exponencial utilizado (sigmasq ≈ 8, rango ≈ 0.00175°) indica que la dependencia espacial tiene alcance en torno a ~200 m; por tanto, las variaciones observadas son consistentes con procesos que actúan a esa escala.

El mapa de desviación estándar muestra que la incertidumbre de la estimación es menor en las áreas con mayor densidad de muestreo y aumenta hacia los bordes del dominio. Esto es consistente con la naturaleza del kriging: la precisión es mayor donde hay puntos observados cercanos y decrece en zonas poco muestreadas. El patrón de incertidumbre sugiere que se debe tener más cautela al interpretar las predicciones en los límites del área de estudio.