Introducción

El cultivo de aguacate (Persea americana) ha adquirido gran importancia en el departamento del Cauca debido a su potencial económico y las condiciones agroclimáticas favorables. La producción y calidad del fruto están fuertemente influenciadas por variables ambientales como la temperatura, humedad relativa, presión atmosférica y velocidad del viento. Por esta razón, analizar la distribución espacial de la temperatura durante las primeras etapas del ciclo productivo permite identificar patrones microclimáticos y orientar decisiones de manejo agrícola basadas en datos.

En este estudio se aplica una metodología geoestadística para caracterizar la variabilidad espacial de la temperatura, evaluar la existencia de autocorrelación espacial y generar un mapa de predicción mediante el método de Kriging Ordinario, utilizando mediciones registradas en una finca productora de aguacate ubicada en el Cauca.

Objetivos

Objetivo general: Analizar la dependencia espacial de la variable temperatura en árboles de aguacate del departamento del Cauca, aplicando técnicas geoestadísticas para generar mapas de predicción.

Objetivos específicos: - Seleccionar y limpiar los datos de la primera fecha de medición. - Diagnosticar la estructura espacial de la temperatura. - Ajustar un modelo teórico de semivariograma (esférico, exponencial, gaussiano). - Interpolar los valores de temperatura mediante Kriging Ordinario. - Evaluar la validez del modelo mediante validación cruzada.

Marco Teórico

La geoestadística permite analizar fenómenos que varían en el espacio, asumiendo que las observaciones más cercanas son más similares que las distantes (principio de autocorrelación espacial de Tobler).
El semivariograma es la herramienta fundamental para cuantificar esta dependencia, mostrando cómo cambia la semivarianza con la distancia entre puntos de muestreo.
El modelo teórico (esférico, exponencial o gaussiano) describe esta relación y permite estimar valores en ubicaciones no muestreadas mediante interpolación Kriging, un método óptimo insesgado que minimiza la varianza del error de predicción.
En el contexto agrícola, estas técnicas ayudan a caracterizar microclimas y apoyar la agricultura de precisión mediante mapas de variabilidad espacial.

Carga y filtrado de datos

A continuación, se presenta la ubicación general de los puntos muestreados en el estudio.
El conjunto de datos corresponde a mediciones georreferenciadas de árboles de aguacate ubicados en el departamento del Cauca, específicamente en una finca productora cercana al municipio de Popayán.
Cada punto representa un árbol con registros simultáneos de variables meteorológicas, siendo la temperatura la variable de interés principal para este análisis.

# =====================================================
# 📊 3. Carga y filtrado de datos
# =====================================================
library(readxl)
library(dplyr)
library(sf)
library(sp)
library(gstat)
library(ggplot2)
library(tidyr)

# Cargar datos
datos <- read_excel("Datos_Completos_Aguacate.xlsx")

# Ver estructura
str(datos)

## tibble [20,271 × 21] (S3: tbl_df/tbl/data.frame)
##  $ id_arbol                      : chr [1:20271] "1" "2" "3" "4" ...
##  $ Latitude                      : num [1:20271] 2.38 2.38 2.38 2.38 2.38 ...
##  $ Longitude                     : num [1:20271] -76.6 -76.6 -76.6 -76.6 -76.6 ...
##  $ FORMATTED_DATE_TIME           : chr [1:20271] "21/08/2019  9:22:57 a, m," "21/08/2019  9:27:13 a, m," "21/08/2019  9:36:36 a, m," "21/08/2019  9:38:02 a, m," ...
##  $ Psychro_Wet_Bulb_Temperature  : num [1:20271] 14.8 11.6 12.9 14.1 14.3 14.2 14.4 12.8 15 14 ...
##  $ Station_Pressure              : num [1:20271] 805 805 806 806 805 ...
##  $ Relative_Humidity             : num [1:20271] 33.6 36.8 31.5 33.2 34.3 33.8 34.9 34.2 33.6 34 ...
##  $ Crosswind                     : num [1:20271] 0.2 3.6 0.4 0.6 0.4 0.5 0.6 2.9 0.4 0.7 ...
##  $ Temperature                   : num [1:20271] 25.7 20.8 23.7 25 25 25 24.9 22.9 26.2 24.6 ...
##  $ Barometric_Pressure           : num [1:20271] 805 805 806 806 805 ...
##  $ Headwind                      : num [1:20271] 0.7 3.5 0.7 0.7 0.4 0.1 0 2.2 0.2 0.7 ...
##  $ Direction_True                : num [1:20271] 166 314 332 139 129 83 93 128 118 137 ...
##  $ Direction_Mag                 : num [1:20271] 165 313 331 139 128 83 92 127 118 137 ...
##  $ Wind_Speed                    : num [1:20271] 0.8 5.1 0.8 0.9 0.6 0.5 0.6 3.7 0.5 1 ...
##  $ Heat_Stress_Index             : num [1:20271] 24.1 19.5 22 23.2 23.3 23.3 23.2 21.5 24.6 22.9 ...
##  $ Altitude                      : num [1:20271] 1896 1895 1889 1890 1894 ...
##  $ Dew_Point                     : num [1:20271] 8.6 5.5 5.8 7.7 8.1 7.9 8.3 6.3 8.9 7.7 ...
##  $ Density_Altitude              : num [1:20271] 2.74 2.57 2.66 2.71 2.71 ...
##  $ Wind_Chill                    : num [1:20271] 25.7 20.8 23.7 24.9 24.9 24.9 24.8 22.9 26.2 24.5 ...
##  $ Estado_Fenologico_Predominante: num [1:20271] 715 715 715 715 715 715 715 715 715 715 ...
##  $ Frutos_Afectados              : num [1:20271] 0 3 0 0 1 0 0 0 0 21 ...

# Filtrar primer periodo (ajusta con la fecha de tu dataset)
primer_periodo <- datos %>%
  filter(FORMATTED_DATE_TIME == "01/10/2020  10:11:12 a, m,") # ajustar según tus datos

# Seleccionar columnas relevantes
datos_sel <- primer_periodo %>%
  select(id_arbol, Latitude, Longitude, 
         Temperatura = Psychro_Wet_Bulb_Temperature, 
         Station_Pressure, Relative_Humidity, Crosswind)

head(datos_sel)

library(leaflet)

leaflet() %>%
  addProviderTiles("CartoDB.Positron") %>%
  addCircleMarkers(
    lng = mean(datos$Longitude, na.rm = TRUE),
    lat = mean(datos$Latitude, na.rm = TRUE),
    radius = 6, color = "darkgreen", fillColor = "limegreen",
    fillOpacity = 0.7, label = "Área de estudio - Finca de aguacate"
  ) %>%
  addScaleBar(position = "bottomleft")

La distribución espacial de los puntos evidencia una cobertura homogénea del área, condición favorable para el análisis geoestadístico posterior.

Exploración espacial y descriptiva

Los datos del primer periodo corresponden a 534 registros georreferenciados, representando la posición de árboles de aguacate muestreados dentro del área de estudio.

La inspección espacial muestra una distribución homogénea de los puntos, lo cual es favorable para el modelado geoestadístico.

El histograma de temperatura evidencia un comportamiento ligeramente simétrico y centrado alrededor de los 21–22 °C, con valores mínimos de 19 °C y máximos cercanos a 25 °C. Esta distribución indica que la variable cumple de manera razonable el supuesto de normalidad requerido por el Kriging, sin necesidad de transformaciones adicionales.

Resumen de distancias geográficas entre puntos de muestreo

# Crear objeto espacial (si no lo tienes)
library(sp)
datos_sp <- datos_sel
coordinates(datos_sp) <- ~Longitude + Latitude

# Calcular distancias euclidianas entre puntos
distancias <- sp::spDists(coordinates(datos_sp))

# Resumen estadístico
dist_summary <- summary(as.vector(distancias))
dist_summary_df <- data.frame(
  Estadístico = names(dist_summary),
  `Valor de Distancia (Grados)` = as.numeric(dist_summary)
)

# Mostrar tabla bonita
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

Antes del cálculo del semivariograma, se evaluó la dispersión espacial de los puntos de muestreo. La tabla anterior muestra el resumen estadístico de las distancias euclidianas entre pares de puntos georreferenciados. Estas distancias, de magnitud homogénea y sin valores extremos, confirman una cobertura espacial adecuada para la estimación de la dependencia espacial mediante el modelo de semivariograma.

Análisis del semivariograma

Análisis Estadístico

# =====================================================
# 5. Análisis del Semivariograma
# =====================================================

# Calcular semivariograma empírico
vgm_emp <- variogram(Temperatura ~ 1, datos_sp)

# Ajustar modelos teóricos
fit_exp <- fit.variogram(vgm_emp, model = vgm("Exp"))
fit_gau <- fit.variogram(vgm_emp, model = vgm("Gau"))
fit_sph <- fit.variogram(vgm_emp, model = vgm("Sph"))

# Calcular errores SSE
loss_exp <- attr(fit_exp, "SSErr")
loss_gau <- attr(fit_gau, "SSErr")
loss_sph <- attr(fit_sph, "SSErr")

# Tabla comparativa
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)
)

knitr::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	0.789368	0.000108	0.538822	1047288643
Gaussiano	0.576587	0.000026	0.719313	107883347357
Esférico	0.634339	0.000261	0.655524	785921931

Desde el punto de vista estadístico:

El SSE (Sum of Squared Errors) mide la pérdida del ajuste; cuanto menor es el SSE, mejor es el modelo. 👉 En este caso, el modelo Esférico tiene el SSE más bajo (7.86×10⁹), indicando el ajuste más preciso a los valores empíricos.

El rango (φ = 0.000261) del modelo esférico es mayor que el de los otros modelos, lo que indica una dependencia espacial más extendida, coherente con la continuidad gradual esperada en la variable temperatura.

El nugget (τ = 0.65) sugiere una pequeña componente de variabilidad no explicada (ruido o error de medición), pero dentro de valores razonables para datos experimentales de campo.

Análisis Visual

# =====================================================
# 5.1 Comparación visual de modelos teóricos (final)
# =====================================================

library(gstat)
library(ggplot2)

# Calcular distancia máxima
max_dist <- max(vgm_emp$dist, na.rm = TRUE)

# Convertir el semivariograma empírico a data.frame
emp_df <- as.data.frame(vgm_emp)

# Calcular valores teóricos de cada modelo hasta la distancia máxima
pred_exp <- variogramLine(fit_exp, maxdist = max_dist)
pred_gau <- variogramLine(fit_gau, maxdist = max_dist)
pred_sph <- variogramLine(fit_sph, maxdist = max_dist)

# Agregar nombres de modelo
pred_exp$model <- "Exponencial"
pred_gau$model <- "Gaussiano"
pred_sph$model <- "Esférico"

# Combinar en un solo data frame
pred_all <- rbind(pred_exp, pred_gau, pred_sph)

# Crear gráfico comparativo con ggplot2
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.position = "bottom",
    legend.title = element_text(face = "bold")
  )

El gráfico comparativo muestra las tres funciones teóricas (Exponencial, Gaussiana y Esférica) superpuestas sobre el semivariograma empírico obtenido a partir de la variable Temperatura. Visualmente, se observa que:

El modelo Gaussiano (curva morada) presenta un ascenso muy abrupto y alcanza el “sill” (meseta) casi de forma inmediata, lo que sugiere una autocorrelación espacial excesiva a distancias muy cortas. Esto indica que sobredimensiona la dependencia local y no representa bien la variación real observada entre los puntos.

El modelo Exponencial (azul) muestra un crecimiento más progresivo que el Gaussiano, pero subestima la semivarianza en las primeras distancias, y su ajuste a los puntos empíricos es menos preciso en la región intermedia (0.0001–0.0003 grados).

El modelo Esférico (rojo) logra el mejor equilibrio visual y estadístico: ajusta con precisión la tendencia creciente inicial del semivariograma y alcanza una meseta suave en concordancia con los valores empíricos, representando adecuadamente la variabilidad espacial observada.

Conclusión: Tanto por criterios visuales (mejor ajuste a los puntos empíricos) como por criterios numéricos (menor SSE y rango coherente con la estructura espacial esperada), el modelo Esférico se identifica como el modelo más adecuado para describir la dependencia espacial de la Temperatura en el cultivo de aguacate.

Este modelo será utilizado para la interpolación mediante Kriging Ordinario, garantizando una estimación más robusta y confiable de la variable en el área de estudio. En conjunto, estos resultados indican una estructura espacial definida que justifica la interpolación mediante Kriging.

Interpretación kriging Ordinario

# =====================================================
# 🌐 6. Interpolación Kriging Ordinario
# =====================================================
# Selección final (elige el mejor modelo según ajuste)
vgm_model <- fit_sph   # ← Usa fit_exp o fit_gau si es el mejor
# Crear grilla de interpolación
bbox <- st_bbox(st_as_sf(datos_sp))
x.range <- seq(bbox["xmin"], bbox["xmax"], length=100)
y.range <- seq(bbox["ymin"], bbox["ymax"], length=100)
grd <- expand.grid(Longitude = x.range, Latitude = y.range)
coordinates(grd) <- ~Longitude + Latitude
gridded(grd) <- TRUE

# Kriging
kriging_result <- krige(Temperatura ~ 1, datos_sp, grd, model=vgm_model)

## [using ordinary kriging]

# Mapa de predicción
spplot(kriging_result["var1.pred"], main="Mapa de Predicción de Temperatura")

# --- Mapa de error (Desviación estándar del Kriging)
plot(kriging_result, type = "variance", main = "Mapa de Error (Desviación estándar del Kriging)")

Análisis

Mapa de predicción - La interpolación mediante Kriging Ordinario permitió estimar la temperatura en toda el área de estudio a partir de los 534 puntos muestreados. El mapa de predicción evidencia zonas con gradientes térmicos suaves, con temperaturas ligeramente más altas en la zona central y valores más bajos hacia los límites del área. Esto puede estar asociado con diferencias topográficas o de exposición al sol.

El histograma de errores del Kriging presenta una distribución centrada en cero con valores comprendidos entre −2.45 °C y +2.63 °C, lo cual indica que el modelo es prácticamente insesgado y presenta errores moderados dentro de un rango aceptable para aplicaciones agrícolas.

Mapa de error

El mapa de error del Kriging muestra la incertidumbre asociada a la estimación de temperatura. Las zonas con errores bajos (en tonos fríos) representan regiones donde la densidad de puntos de muestreo es mayor y, por tanto, las predicciones son más confiables. Por el contrario, las áreas con errores altos (en tonos cálidos) indican regiones con menor soporte de datos y mayor incertidumbre en la interpolación.

Validación cruzada

# =====================================================
# 🧪 7. Validación cruzada
# =====================================================

cv <- krige.cv(Temperatura ~ 1, datos_sp, model = vgm_model)
summary(cv$residual)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -2.451188 -0.609707 -0.010900 -0.001691  0.555265  2.632204

hist(cv$residual, main="Distribución de errores", col="lightblue")

Análisis validación cruzada del modelo

Para evaluar la calidad del modelo de interpolación, se aplicó una validación cruzada (Leave-One-Out), donde cada punto fue predicho utilizando los demás datos.
Los errores de predicción se distribuyeron de manera aproximadamente normal y centrada en cero, confirmando la insesgadez del modelo.
La desviación estándar de los errores fue inferior a 1.5 °C, lo que indica una buena precisión del Kriging Ordinario aplicado.

Conclusiones

La temperatura presenta una correlación espacial significativa, evidenciada por la forma ascendente del semivariograma y el ajuste satisfactorio del modelo esférico.
El Kriging Ordinario resultó un método adecuado para generar un mapa continuo de la variable, permitiendo identificar zonas térmicamente homogéneas dentro del cultivo.
Los errores de predicción son bajos y no muestran sesgo sistemático, lo que demuestra la validez del modelo predictivo.
Este análisis constituye una base sólida para incorporar otras variables (humedad, altitud o radiación solar) en futuros estudios de agricultura de precisión y zonificación agroclimática.

Análisis Geoestadístico de Temperatura en Cultivos de Aguacate – Cauca

María Hernández

10/11/2025