1 Introducción

El conjunto de datos utilizado en este trabajo corresponde a mediciones ambientales realizadas en una finca productora de aguacate ubicada en el departamento del Cauca. Cada registro representa un árbol individual, georreferenciado mediante sus coordenadas espaciales y asociado a diversas variables ambientales medidas con sensores portátiles, tales como temperatura, humedad relativa, velocidad del viento y altura.

Para este análisis, se trabajará exclusivamente con el primer periodo de muestreo, correspondiente al 01/10/2020, lo que permite contar con una base homogénea y comparable de aproximadamente 534 observaciones. El objetivo principal es estudiar la distribución espacial de la temperatura, variable climática crítica en procesos agrícolas, ya que influye en el crecimiento, productividad y estrés fisiológico de los cultivos.

La metodología empleada se basa en técnicas de geoestadística, que permiten analizar la dependencia espacial entre mediciones y generar predicciones en áreas no muestreadas. En primera instancia, se realizará un análisis exploratorio de datos para comprender la variabilidad y estructura espacial de la temperatura. Luego, se construirá el semivariograma experimental, herramienta fundamental para evaluar la autocorrelación espacial. A partir de este, se ajustarán diferentes modelos teóricos —como el esférico, exponencial o gaussiano— con el fin de identificar el que mejor representa el comportamiento espacial de la variable.

Finalmente, utilizando la metodología de kriging, se generará una predicción espacial continua de la temperatura sobre el área de estudio, obteniendo un mapa interpolado que permitirá visualizar su distribución espacial. En caso de que la temperatura no presente autocorrelación significativa, se considerarán otras variables ambientales disponibles.

Este proceso permitirá no solo caracterizar el comportamiento espacial de la variable seleccionada, sino también aplicar de forma completa el flujo de trabajo propio del análisis geoestadístico. Si deseas, también puedo redactar la sección de metodología o preparar el código en R paso a paso.

library(dplyr)
library(geoR)
library(sp)
#install.packages("gstat")
library(gstat)
library(ggplot2)
library(knitr)
library(readxl)

2 Carga, filtrado y creación del subconjunto datos_2020

##  [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"
## [1] 534

Tras cargar la base de datos original (Datos_Completos_Aguacate), se realizó un proceso de limpieza y preparación de la variable de fecha. La columna FORMATTED_DATE_TIME presentaba caracteres adicionales (comas y espacios), por lo que fue estandarizada y transformada al formato de fecha (Date). Posteriormente, se exploraron todas las fechas disponibles en el conjunto, identificándose múltiples días de muestreo distribuidos entre 2019 y 2020.

El análisis requería trabajar exclusivamente con el primer periodo seleccionado, correspondiente al 01 de octubre de 2020, por lo que se filtraron todas las mediciones asociadas a ese día. El resultado fue un subconjunto de 534 registros, coincidiendo con lo esperado para este periodo de muestreo. Este dataset representa mediciones tomadas a nivel de árbol, cada uno georreferenciado, lo cual lo convierte en un insumo adecuado para la aplicación de métodos geoestadísticos.

Finalmente, se creó la variable Temperatura, equivalente a la columna original Temperature, con el fin de facilitar su uso en los paquetes gstat y geoR. Esta preparación asegura un formato limpio y consistente para el análisis posterior, y establece una base sólida para las etapas de exploración, modelado de dependencia espacial y predicción mediante kriging.

3 Creación del geodato y análisis de tendencia espacial

library(geoR)
library(ggplot2)
library(viridis)   # escala de colores agradable
library(patchwork) # para combinar gráficos

# 3.1 Crear objeto geodatos_temp (geoR)
geodatos_temp <- geoR::as.geodata(datos_2020,
                                  coords.col = c("Longitude", "Latitude"),
                                  data.col   = "Temperature")

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
# 3.2 Modelo lineal simple para tendencia
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
# ---- Estilo base para todos los gráficos ----
estilo_base <- theme_minimal(base_size = 12) +
  theme(
    plot.title    = element_text(face = "bold", hjust = 0.5),
    panel.grid.minor = element_blank()
  )

# 1) Mapa de puntos coloreados por temperatura
p1 <- ggplot(datos_2020,
             aes(x = Longitude, y = Latitude, colour = Temperature)) +
  geom_point(size = 1.5, alpha = 0.8) +
  scale_colour_viridis(option = "plasma", name = "°C") +
  labs(
    title = "Distribución espacial de la temperatura",
    x = "Longitud", y = "Latitud"
  ) +
  coord_equal() +
  estilo_base

# 2) Tendencia según Latitud
p2 <- ggplot(datos_2020,
             aes(x = Latitude, y = Temperature)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", se = FALSE, colour = "red", linewidth = 1) +
  labs(
    title = "Tendencia de temperatura según Latitud",
    x = "Latitud", y = "Temperatura (°C)"
  ) +
  estilo_base

# 3) Tendencia según Longitud
p3 <- ggplot(datos_2020,
             aes(x = Longitude, y = Temperature)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", se = FALSE, colour = "blue", linewidth = 1) +
  labs(
    title = "Tendencia de temperatura según Longitud",
    x = "Longitud", y = "Temperatura (°C)"
  ) +
  estilo_base

# 4) Histograma + densidad
p4 <- ggplot(datos_2020, aes(x = Temperature)) +
  geom_histogram(aes(y = ..density..),
                 bins = 20,
                 fill = "palegreen3",
                 colour = "grey30") +
  geom_density(linewidth = 1.1, colour = "darkgreen") +
  labs(
    title = "Distribución estadística de la temperatura",
    x = "Temperatura (°C)", y = "Densidad"
  ) +
  estilo_base

# Panel 2x2
(p1 | p2) /
(p3 | p4)

El objeto geodatos_temp permitió integrar las coordenadas geográficas y los valores de temperatura del 01/10/2020 en una estructura adecuada para la aplicación de métodos geoestadísticos. Según el resumen del objeto, el conjunto abarca un rango espacial reducido, cercano a 0.0016 grados en longitud y 0.0015 grados en latitud, lo cual es coherente con mediciones realizadas dentro de una parcela agrícola. En total se registraron 534 observaciones válidas, con temperaturas entre 22.2 °C y 29.7 °C, y un valor mediano de 25.8 °C.

Los gráficos de distribución espacial muestran que los árboles están ubicados de manera homogénea dentro del área de estudio, lo cual es favorable para el análisis espacial. Además, el mapa de puntos coloreados por temperatura evidencia la presencia de zonas con valores ligeramente más altos hacia el sector superior derecho del terreno, mientras que las temperaturas más bajas tienden a concentrarse en la parte inferior izquierda. Este patrón sugiere la influencia de gradientes microclimáticos dentro de la finca.

Para evaluar formalmente la presencia de tendencia espacial en la variable, se ajustó un modelo lineal considerando la latitud y la longitud como predictores. Ambos coeficientes resultaron estadísticamente significativos, lo que indica que la temperatura varía sistemáticamente con la posición en el espacio. Esto coincide con los gráficos de tendencia, donde se observa que la temperatura disminuye ligeramente hacia el norte y hacia el este del área de estudio.

El modelo lineal presentó un valor de R2 aproximado de 0.082, lo cual indica que la tendencia espacial explica una fracción pequeña pero consistente de la variación observada. Aunque la magnitud de esta tendencia no es muy elevada, sí evidencia la existencia de una estructura espacial que debe ser considerada en el análisis geoestadístico.

En conjunto, tanto la exploración visual como los resultados del modelo lineal muestran que la temperatura presenta una tendencia espacial moderada, lo que justifica continuar con la construcción del semivariograma y la posterior aplicación del método de kriging para la estimación espacial.

4 Semivariograma experimental

# Distancias entre puntos (referencia en consola)
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
# Semivariograma experimental
variograma_temp <- geoR::variog(
  geodatos_temp,
  max.dist = 0.0015
)
## variog: computing omnidirectional variogram
# Envolvente por simulación
set.seed(123)
variograma_env <- geoR::variog.mc.env(
  geodatos_temp,
  obj.variog = 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
# Pasar a data.frame para ggplot
df_vario <- data.frame(
  dist  = variograma_temp$u,
  gamma = variograma_temp$v
)

df_env <- data.frame(
  dist  = variograma_env$u,
  lower = variograma_env$v.lower,
  upper = variograma_env$v.upper
)

# Gráfico mejorado
ggplot(df_vario, aes(x = dist, y = gamma)) +
  # banda de la envolvente (sin heredar el y = gamma)
  geom_ribbon(
    data = df_env,
    inherit.aes = FALSE,
    aes(x = dist, ymin = lower, ymax = upper),
    fill = "grey80",
    alpha = 0.6
  ) +
  # puntos + línea del semivariograma
  geom_point(size = 2, colour = "black") +
  geom_line(linewidth = 0.7, colour = "darkgreen") +
  labs(
    title = "Semivariograma experimental de Temperatura",
    x = "Distancia",
    y = expression(gamma(h))
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5),
    panel.grid.minor = element_blank()
  )

El semivariograma experimental permite evaluar la dependencia espacial de la temperatura dentro del área de estudio. En este caso, el semivariograma presenta un comportamiento ascendente claro: a medida que aumenta la distancia entre los árboles, también aumenta la semivarianza. Este patrón sugiere que existe autocorrelación espacial positiva, es decir, que los valores de temperatura tienden a ser más similares entre puntos cercanos que entre puntos alejados.

Se observa un crecimiento rápido en las primeras distancias, seguido de una tendencia hacia la estabilización alrededor de valores de gamma cercanos a 3.5. Esta meseta indica la aproximación al sill, lo que sugiere que más allá de cierta distancia la temperatura deja de mostrar similitud espacial significativa. Este comportamiento es típico en variables ambientales afectadas por gradientes locales como exposición solar o circulación de aire.

La banda gris correspondiente a la envolvente de simulación (generada mediante 99 permutaciones aleatorias) permite evaluar la significancia de la autocorrelación. La mayor parte del semivariograma observado se encuentra por encima de esta banda, lo que confirma que la estructura espacial detectada no es producto del azar, sino que refleja un patrón real en los datos.

En conjunto, el semivariograma experimental muestra una estructura espacial bien definida, lo que justifica continuar con la etapa de ajuste del modelo teórico y la aplicación del método de kriging para generar la predicción espacial de la temperatura.

5 Ajuste de modelos teóricos (geoR)

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

# 2. Ajuste de modelos (usa variograma_temp ya calculado antes)
model_exp <- geoR::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
model_gaus <- geoR::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
model_esf <- geoR::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
# 3. Semivariograma experimental a data.frame
df_vario <- data.frame(
  dist  = variograma_temp$u,
  gamma = variograma_temp$v
)

# 4. Secuencia de distancias para dibujar los modelos
dist_seq <- seq(0, max(df_vario$dist), length.out = 200)

# 5. Funciones teóricas de semivariograma
gamma_exp <- function(h, sigma2, phi, tau2) {
  tau2 + sigma2 * (1 - exp(-h / phi))
}

gamma_gaus <- function(h, sigma2, phi, tau2) {
  tau2 + sigma2 * (1 - exp(-(h^2) / (phi^2)))
}

gamma_esf <- function(h, sigma2, phi, tau2) {
  # modelo esférico clásico
  g <- ifelse(
    h <= phi,
    tau2 + sigma2 * (1.5 * (h / phi) - 0.5 * (h / phi)^3),
    tau2 + sigma2
  )
  g
}

# 6. Data frame con las curvas de los tres modelos
df_modelos <- rbind(
  data.frame(
    dist  = dist_seq,
    gamma = gamma_exp(dist_seq,
                      sigma2 = model_exp$cov.pars[1],
                      phi    = model_exp$cov.pars[2],
                      tau2   = model_exp$nugget),
    modelo = "Exponencial"
  ),
  data.frame(
    dist  = dist_seq,
    gamma = gamma_gaus(dist_seq,
                       sigma2 = model_gaus$cov.pars[1],
                       phi    = model_gaus$cov.pars[2],
                       tau2   = model_gaus$nugget),
    modelo = "Gaussiano"
  ),
  data.frame(
    dist  = dist_seq,
    gamma = gamma_esf(dist_seq,
                      sigma2 = model_esf$cov.pars[1],
                      phi    = model_esf$cov.pars[2],
                      tau2   = model_esf$nugget),
    modelo = "Esférico"
  )
)

# 7. Gráfico final: experimental + modelos
ggplot() +
  # puntos del semivariograma experimental
  geom_point(data = df_vario,
             aes(x = dist, y = gamma),
             size = 2, alpha = 0.8) +
  # líneas de los modelos
  geom_line(data = df_modelos,
            aes(x = dist, y = gamma, colour = modelo),
            linewidth = 1.1) +
  scale_colour_manual(
    values = c(
      "Exponencial" = "#1f77b4",
      "Gaussiano"   = "#d62728",
      "Esférico"    = "#9467bd"
    )
  ) +
  labs(
    title    = "Ajuste de modelos de semivariograma",
    subtitle = "Comparación de modelos exponencial, gaussiano y esférico",
    x        = "Distancia",
    y        = expression(gamma(h)),
    colour   = "Modelo"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title    = element_text(face = "bold", hjust = 0.5),
    plot.subtitle = element_text(hjust = 0.5),
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )

A partir del semivariograma experimental de la temperatura se procedió al ajuste de tres modelos teóricos de semivarianza: exponencial, gaussiano y esférico. El ajuste se realizó mediante la función variofit del paquete geoR, utilizando mínimos cuadrados ponderados y una grilla de valores iniciales para los parámetros de cada modelo.

En el gráfico se observa la comparación entre los puntos del semivariograma experimental y las curvas teóricas ajustadas. El modelo exponencial es el que sigue de forma más cercana la tendencia de los datos, especialmente en el rango de distancias cortas e intermedias, donde la autocorrelación espacial es más relevante para el kriging. El modelo gaussiano presenta una curva más suave que subestima la semivarianza en las primeras distancias y tiende a sobrestimarla en las distancias más grandes. Por su parte, el modelo esférico se separa con mayor rapidez del semivariograma experimental, en particular a partir de la mitad del rango de distancias, mostrando un ajuste menos consistente.

Los valores de la función de pérdida (loss value) reportados por variofit confirman esta apreciación visual: el modelo exponencial presenta el menor error de ajuste, seguido por el modelo gaussiano, mientras que el esférico muestra el error más alto. Aunque se emitieron advertencias sobre valores iniciales elevados del parámetro de varianza, el procedimiento de optimización logró converger en soluciones razonables que describen adecuadamente la estructura espacial observada.

En síntesis, la combinación de la evidencia gráfica y numérica indica que el modelo exponencial es el que mejor representa la dependencia espacial de la temperatura para este conjunto de datos. Por esta razón, dicho modelo se selecciona como base para la etapa de interpolación mediante kriging en la sección siguiente.

6 Interpolación espacial mediante Kriging y generación de mapas

library(ggplot2)
library(viridis)
library(patchwork)

# 1. Rangos de coordenadas
range_long <- range(geodatos_temp$coords[, 1])
range_lat  <- range(geodatos_temp$coords[, 2])

# 2. Malla 100x100
geodatos_grid <- expand.grid(
  lon = seq(range_long[1], range_long[2], length = 100),
  lat = seq(range_lat[1],  range_lat[2],  length = 100)
)

# (opcional) Visualización rápida de la malla + puntos de muestreo
coords_df <- data.frame(
  lon = geodatos_temp$coords[, 1],
  lat = geodatos_temp$coords[, 2]
)

p_grid <- ggplot() +
  geom_point(data = geodatos_grid,
             aes(x = lon, y = lat),
             size = 0.2, alpha = 0.5) +
  geom_point(data = coords_df,
             aes(x = lon, y = lat),
             colour = "red", size = 1.5, alpha = 0.9) +
  labs(
    title = "Malla de predicción (100x100)",
    x = "Longitud", y = "Latitud"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5),
    panel.grid.minor = element_blank()
  ) +
  coord_equal()

# 3. Kriging con el modelo exponencial ajustado antes (model_exp)
geodatos_ko <- geoR::krige.conv(
  geodata = geodatos_temp,
  loc     = geodatos_grid,
  krige   = geoR::krige.control(
    nugget    = model_exp$nugget,
    trend.d   = "cte",
    trend.l   = "cte",
    cov.pars  = model_exp$cov.pars,
    cov.model = "exponential"
  )
)
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
# 4. Data frame con predicción y desviación estándar
df_ko <- data.frame(
  lon  = geodatos_grid$lon,
  lat  = geodatos_grid$lat,
  pred = geodatos_ko$predict,
  sd   = sqrt(geodatos_ko$krige.var)
)

theme_mapa <- theme_minimal(base_size = 12) +
  theme(
    plot.title    = element_text(face = "bold", hjust = 0.5),
    plot.subtitle = element_text(hjust = 0.5),
    panel.grid.minor = element_blank(),
    legend.position  = "right"
  )

# 5. Mapa de predicción
p_pred <- ggplot(df_ko, aes(x = lon, y = lat)) +
  geom_raster(aes(fill = pred)) +
  geom_point(data = coords_df,
             aes(x = lon, y = lat),
             colour = "black", size = 0.6, alpha = 0.9) +
  scale_fill_viridis(
    option = "plasma",
    name   = "Temperatura"
  ) +
  labs(
    title = "Predicción de Temperatura (Kriging)",
    x = "Longitud", y = "Latitud"
  ) +
  coord_equal() +
  theme_mapa

# 6. Mapa de desviación estándar
p_sd <- ggplot(df_ko, aes(x = lon, y = lat)) +
  geom_raster(aes(fill = sd)) +
  geom_point(data = coords_df,
             aes(x = lon, y = lat),
             colour = "black", size = 0.6, alpha = 0.9) +
  scale_fill_viridis(
    option = "magma",
    name   = "Desv. estándar"
  ) +
  labs(
    title = "Incertidumbre de la predicción",
    x = "Longitud", y = "Latitud"
  ) +
  coord_equal() +
  theme_mapa

# 7. Panel final: malla (opcional) + mapas
(p_grid | p_pred | p_sd)

La interpolación espacial se realizó utilizando el método de kriging ordinario, empleando el modelo exponencial previamente identificado como el mejor ajuste para describir la dependencia espacial de la temperatura. Para ello, se creó una malla regular de 100 por 100 puntos que cubre completamente el rango de coordenadas de las observaciones originales. Esta malla sirve como base para estimar la temperatura en ubicaciones no muestreadas dentro del área de estudio.

El primer mapa muestra la distribución de la malla junto con los puntos de muestreo. Se observa una cobertura espacial densa y homogénea de los árboles, lo que garantiza estabilidad en la interpolación. La malla regular, por su parte, permite una representación continua del espacio interpolado.

El mapa resultante del kriging revela un patrón térmico no uniforme dentro de la finca. Se identifican zonas con temperaturas más elevadas en sectores específicos, mientras que otras áreas presentan valores más bajos. Estos gradientes de temperatura son consistentes con lo observado en la exploración inicial y pueden estar asociados a microvariaciones ambientales dentro del terreno, como diferencias en exposición solar, relieve o circulación de aire.

El tercer mapa corresponde a la desviación estándar de la predicción, la cual representa la incertidumbre del kriging. Como es esperable, la incertidumbre es menor en las zonas con mayor densidad de puntos de muestreo y aumenta hacia los extremos de la malla, donde la distancia a las observaciones es mayor. En general, los niveles de incertidumbre son bajos, lo que indica que la predicción es confiable en buena parte del dominio espacial.

En conjunto, los mapas generados demuestran que el kriging ordinario fue capaz de capturar la estructura espacial de la temperatura y producir una superficie continua que refleja adecuadamente los patrones térmicos del área. Además, la incertidumbre obtenida confirma que las estimaciones son robustas dentro del rango cubierto por los datos originales.

7 Resumen de distancias geográficas entre puntos muestreados

# Seleccionar columnas relevantes
datos_sel <- datos_2020 %>%
  dplyr::select(Longitude, Latitude, Temperatura)

# Crear objeto espacial
datos_sp <- datos_sel
sp::coordinates(datos_sp) <- ~Longitude + Latitude

# Matriz de distancias
distancias <- sp::spDists(sp::coordinates(datos_sp))
dist_summary <- summary(as.vector(distancias))

dist_summary_df <- data.frame(
  Estadístico                    = names(dist_summary),
  `Valor de Distancia (Grados)`  = as.numeric(dist_summary)
)

knitr::kable(
  dist_summary_df,
  caption = "Tabla 1. Resumen de Distancias Geográficas entre Registros Filtrados",
  digits  = 8,
  align   = "lc"
)
Tabla 1. Resumen de Distancias Geográficas entre Registros Filtrados
Estadístico Valor.de.Distancia..Grados.
Min. 0.00000000
1st Qu. 0.00040381
Median 0.00063996
Mean 0.00068140
3rd Qu. 0.00091711
Max. 0.00195913

El cálculo de las distancias euclidianas entre todos los pares de puntos georreferenciados permite evaluar la distribución espacial del muestreo dentro de la finca. La tabla de resumen muestra que la distancia mínima entre puntos es prácticamente cero, lo cual corresponde a la proximidad inmediata entre árboles vecinos. Por otro lado, la distancia máxima registrada es cercana a 0.00196 grados, lo que equivale a un rango espacial coherente con la extensión reducida del área de estudio.

Los valores del primer y tercer cuartil indican que la mayoría de las distancias se encuentran entre 0.00040 y 0.00091 grados, lo cual evidencia una distribución relativamente homogénea de los árboles en el terreno. Esta uniformidad en el espaciamiento es beneficiosa para los métodos geoestadísticos, ya que proporciona suficiente información en las diferentes zonas del área, evitando vacíos que dificultarían la interpolación.

En conjunto, las distancias observadas confirman que el muestreo espacial presenta buena densidad y cobertura dentro del dominio. Esto respalda la estabilidad del semivariograma y contribuye a la fiabilidad del proceso de kriging utilizado en las etapas siguientes.

8 Comparación de modelos teóricos del semivariograma

library(gstat)
library(ggplot2)
library(knitr)

# Semivariograma empírico (gstat)
vgm_emp <- gstat::variogram(Temperatura ~ 1, datos_sp)

# Estimaciones iniciales a partir del empírico
max_dist   <- max(vgm_emp$dist,  na.rm = TRUE)
sill_est   <- max(vgm_emp$gamma, na.rm = TRUE)
nugget_est <- min(vgm_emp$gamma, na.rm = TRUE)
range_est  <- max_dist / 3

vgm_ini_exp <- gstat::vgm(
  psill  = sill_est - nugget_est,
  model  = "Exp",
  range  = range_est,
  nugget = nugget_est
)

vgm_ini_gau <- gstat::vgm(
  psill  = sill_est - nugget_est,
  model  = "Gau",
  range  = range_est,
  nugget = nugget_est
)

vgm_ini_sph <- gstat::vgm(
  psill  = sill_est - nugget_est,
  model  = "Sph",
  range  = range_est,
  nugget = nugget_est
)

# Ajuste de modelos teóricos
fit_exp <- gstat::fit.variogram(vgm_emp, model = vgm_ini_exp)
fit_gau <- gstat::fit.variogram(vgm_emp, model = vgm_ini_gau)
fit_sph <- gstat::fit.variogram(vgm_emp, model = vgm_ini_sph)

# Función para corregir rangos SOLO en modelos distintos de Nugget
fix_range <- function(fit, fallback_range) {
  idx_bad <- fit$model != "Nug" & (fit$range <= 0 | is.na(fit$range))
  fit$range[idx_bad] <- fallback_range
  fit
}

fit_exp <- fix_range(fit_exp, range_est)
fit_gau <- fix_range(fit_gau, range_est)
fit_sph <- fix_range(fit_sph, range_est)

# SSE de cada modelo
loss_exp <- attr(fit_exp, "SSErr")
loss_gau <- attr(fit_gau, "SSErr")
loss_sph <- attr(fit_sph, "SSErr")

# Tabla de parámetros
tabla_modelos <- data.frame(
  Modelo                     = c("Exponencial", "Gaussiano", "Esférico"),
  `Sill Parcial (φ²)`        = c(fit_exp$psill[2], fit_gau$psill[2], fit_sph$psill[2]),
  `Rango (φ)`                = c(fit_exp$range[2], fit_gau$range[2], fit_sph$range[2]),
  `Nugget (τ²)`              = c(fit_exp$psill[1], fit_gau$psill[1], fit_sph$psill[1]),
  `Función de Pérdida (SSE)` = c(loss_exp,       loss_gau,       loss_sph)
)

kable(
  tabla_modelos,
  caption = "Tabla 3. Parámetros finales de los modelos de semivariograma ajustados - Temperatura",
  digits  = 6
)
Tabla 3. Parámetros finales de los modelos de semivariograma ajustados - Temperatura
Modelo Sill.Parcial..φ.. Rango..φ. Nugget..τ.. Función.de.Pérdida..SSE.
Exponencial 2.232047 0.000237 0.554350 9.083002e+90
Gaussiano 0.798355 0.000237 1.028051 2.246121e+12
Esférico 2.169462 0.000275 0.732679 1.228278e+10 
# Comparación visual con ggplot2
emp_df <- as.data.frame(vgm_emp)

pred_exp <- gstat::variogramLine(fit_exp, maxdist = max_dist)
pred_gau <- gstat::variogramLine(fit_gau, maxdist = max_dist)
pred_sph <- gstat::variogramLine(fit_sph, maxdist = max_dist)

pred_exp$model <- "Exponencial"
pred_gau$model <- "Gaussiano"
pred_sph$model <- "Esférico"

pred_all <- rbind(pred_exp, pred_gau, pred_sph)

ggplot() +
  geom_point(data = emp_df,
             aes(x = dist, y = gamma),
             color = "black", size = 1.8, alpha = 0.6) +
  geom_line(data = pred_all,
            aes(x = dist, y = gamma, color = model),
            linewidth = 1.2) +
  scale_color_manual(values = c(
    "Exponencial" = "#2E86C1",
    "Gaussiano"   = "#8E44AD",
    "Esférico"    = "#E74C3C"
  )) +
  labs(
    title = "Comparación de modelos teóricos del semivariograma - Temperatura",
    x     = "Distancia",
    y     = "Semivarianza",
    color = "Modelo"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title     = element_text(face = "bold", hjust = 0.5),
    legend.title   = element_text(face = "bold"),
    legend.position = "bottom"
  )

A partir del semivariograma empírico obtenido con el paquete gstat, se ajustaron tres modelos teóricos comúnmente utilizados en geoestadística: el modelo exponencial, el modelo gaussiano y el modelo esférico. Cada uno fue ajustado mediante mínimos cuadrados con el objetivo de identificar cuál describe de manera más adecuada la estructura espacial de la temperatura.

El gráfico comparativo muestra diferencias claras entre los tres modelos. El modelo esférico (línea roja) presenta un ascenso inicial más rápido y alcanza la meseta de la semivarianza de manera más abrupta, lo cual coincide con el comportamiento observado en los datos empíricos para distancias cortas y medias. El modelo exponencial (línea azul) exhibe un incremento más gradual y tiende a subestimar la semivarianza en las distancias donde el empírico ya se estabiliza. Por su parte, el modelo gaussiano (línea morada) tiene un crecimiento inicial aún más suave y queda por debajo de los valores empíricos durante casi todo el rango, lo que refleja un ajuste menos preciso.

Los valores de la función de pérdida (SSE) muestran que el modelo esférico logra el ajuste más adecuado entre los tres evaluados, seguido del modelo gaussiano. Aunque los valores absolutos de SSE son altos debido a la escala de la semivarianza, lo que importa es la comparación relativa: el modelo esférico presenta el menor error y, por lo tanto, captura mejor la estructura espacial presente en la variable de temperatura.

En conjunto, el análisis gráfico y numérico confirma que el modelo esférico representa de manera más fiel la dependencia espacial observada. Por esta razón, este modelo se considera el más adecuado para la etapa de interpolación espacial mediante kriging o para la validación cruzada, dependiendo de la metodología seleccionada en el estudio.

9 Validación cruzada

9.1 Análisis de la validación cruzada del modelo de kriging

La validación cruzada tipo leave-one-out permite evaluar el desempeño predictivo del modelo de kriging ajustado. A partir de los resultados obtenidos:

Error Medio (ME): –0.0056 El ME es prácticamente cero y de signo negativo, lo que indica que el modelo no presenta sesgos sistemáticos importantes. En otras palabras, el kriging no tiende a sobreestimar ni subestimar la temperatura de forma consistente, lo cual es una propiedad deseable en un modelo espacial.

Raíz del Error Cuadrático Medio (RMSE): 1.093 El RMSE muestra que el error promedio de las predicciones se encuentra alrededor de 1.1 °C. Dado el rango de temperaturas en el dataset (aprox. 22 °C a 30 °C), este error corresponde a un valor relativamente pequeño respecto al total de variación. Esto sugiere que el modelo tiene buena capacidad predictiva dentro del área muestreada.

MSNE: NaN El MSNE (Mean Standardized Normalized Error) aparece como NaN, lo que suele ocurrir cuando alguna desviación estándar de predicción (var1.sd) es cero o extremadamente pequeña. Este indicador evalúa si las varianzas de kriging están correctamente calibradas; al no poder calcularse, no es posible juzgar si la incertidumbre fue estimada de manera adecuada. Sin embargo, el buen comportamiento del ME y el RMSE indica que, aun sin MSNE, el modelo funciona razonablemente bien en términos de precisión.

9.2 Análisis de desempeño de los modelos de Kriging

La tabla muestra tres métricas principales para evaluar la calidad de los modelos de semivariograma utilizados en el kriging: Error Medio (ME), Raíz del Error Cuadrático Medio (RMSE) y Error Medio Normalizado Cuadrático (MSNE).

La evaluación comparativa entre los modelos Exponencial, Gaussiano y Esférico se realizó mediante tres métricas: el Error Medio (ME), la Raíz del Error Cuadrático Medio (RMSE) y el Error Medio Normalizado Cuadrático (MSNE).

En primer lugar, los valores del ME fueron muy cercanos a cero en los tres casos, lo que indica que ningún modelo presenta un sesgo sistemático significativo. Es decir, las predicciones no tienden a sobreestimar ni subestimar de forma consistente los valores observados.

En cuanto al RMSE, que mide la precisión general de las predicciones, el modelo Exponencial presentó el menor error (RMSE ≈ 1.07), seguido muy de cerca por el modelo Esférico (RMSE ≈ 1.09). Por su parte, el modelo Gaussiano mostró el mayor error (RMSE ≈ 1.25), lo que sugiere un desempeño inferior en términos de exactitud.

Respecto al MSNE, los tres modelos arrojaron valores NaN, lo cual se debe a que esta métrica requiere estimaciones válidas de la varianza de predicción que, en este caso, no fueron computables. Sin embargo, esta limitación no afecta la comparación general entre modelos, ya que ME y RMSE aportan información suficiente para evaluar su desempeño.

En conjunto, los resultados indican que el modelo Exponencial es el que mejor representa la estructura espacial de la variable, mostrando el menor error global y una predicción más estable. El modelo Esférico se comporta de manera similar y constituye una alternativa válida, mientras que el modelo Gaussiano presenta el rendimiento más bajo, por lo que no resulta recomendable para la etapa final de interpolación mediante kriging.

10 Conclusiones

El análisis geoestadístico realizado sobre las mediciones de temperatura del 01/10/2020 en la finca de aguacate del Cauca permitió caracterizar de manera detallada su comportamiento espacial y evaluar la capacidad predictiva de distintos modelos de semivariograma y kriging.

En primer lugar, la exploración descriptiva mostró que la temperatura presenta un rango relativamente acotado (aprox. 22–30 °C) y una distribución cercana a la normalidad, sin valores atípicos extremos. La disposición de los 534 árboles fue densa y homogénea dentro del área de estudio, lo que favoreció la estimación de la estructura espacial y la posterior interpolación.

El análisis de tendencia evidenció la existencia de un gradiente espacial moderado: la temperatura varía de forma sistemática con la posición, con ligeras disminuciones hacia el norte y el este del terreno. Aunque el modelo lineal solo explicó una fracción limitada de la variabilidad total, sí confirmó la presencia de patrones espaciales que justifican el uso de técnicas geoestadísticas.

El semivariograma experimental mostró una estructura espacial bien definida, con aumento de la semivarianza a medida que crece la distancia entre puntos y posterior aproximación a un sill. La envolvente de simulación indicó que la mayoría de los valores observados se encuentran fuera de la banda de aleatoriedad, confirmando la existencia de autocorrelación espacial significativa en la temperatura.

Al ajustar distintos modelos teóricos (exponencial, gaussiano y esférico), tanto con geoR como con gstat, se observó que los modelos exponencial y esférico ofrecen los mejores desempeños. En general, el modelo exponencial siguió más de cerca la forma del semivariograma experimental y presentó el menor error en términos de ajuste y validación cruzada (RMSE ≈ 1.07 °C), mientras que el modelo esférico mostró un comportamiento muy similar y se consolidó como una alternativa válida. El modelo gaussiano, en cambio, presentó errores mayores y un ajuste menos consistente, por lo que no se recomienda para la interpolación final.

La interpolación mediante kriging ordinario, basada en el modelo de semivariograma seleccionado, permitió generar un mapa continuo de temperatura sobre la finca. Este mapa reveló la presencia de microzonas ligeramente más cálidas y más frías, coherentes con los gradientes detectados en la exploración inicial, y acompañado de un mapa de desviación estándar que mostró niveles de incertidumbre bajos en la mayor parte del área, especialmente en las zonas con mayor densidad de datos.

Finalmente, la validación cruzada tipo leave-one-out indicó que los modelos de kriging no presentan sesgo sistemático (ME cercano a cero) y que el error promedio de predicción se mantiene alrededor de 1 °C, lo cual es razonable considerando el rango total de temperaturas registradas. Aunque la métrica MSNE no pudo evaluarse adecuadamente, las evidencias combinadas de ME y RMSE confirman que el modelo seleccionado ofrece una capacidad predictiva adecuada para el contexto de estudio.

En conjunto, los resultados muestran que la temperatura en la finca presenta una estructura espacial clara y aprovechable, y que el uso de técnicas geoestadísticas —en particular el kriging basado en un modelo exponencial o esférico— constituye una herramienta útil para describir y predecir patrones térmicos a escala de parcela. Esto abre la puerta a aplicaciones prácticas en la gestión agronómica, como la identificación de zonas con condiciones microclimáticas diferenciadas, la planificación de monitoreos futuros y el diseño de estrategias de manejo más precisas. No obstante, futuras investigaciones podrían incorporar otros periodos de muestreo y variables ambientales adicionales para evaluar la estabilidad temporal de estos patrones y enriquecer la interpretación agronómica de los resultados.