Caso: Interpolación Climatica Cultivo de Aguacate

# library
lapply(c(
  "sf", "readxl", "dplyr", "stringr", "summarytools",
  "gstat", "RColorBrewer", "mapview", "leafsync", "kableExtra"
), library, character.only = TRUE)

Definición del problema objeto de estudio

El objetivo central de este análisis es comprender la variabilidad espacial de la temperatura registrada en las ubicaciones de árboles de aguacate en el Cauca, Colombia, a partir de datos georreferenciados recolectados en una jornada puntual de monitoreo (01/10/2020). Este ejercicio busca identificar patrones espaciales que permitan modelar la dependencia entre observaciones y fundamentar posteriores procesos de interpolación espacial, clave para apoyar la toma de decisiones agrícolas y la planificación de cultivos en el territorio.

# data
df_aguacate <- readxl::read_excel("Data/Datos_Completos_Aguacate.xlsx") %>% 
mutate(Fecha = str_extract(FORMATTED_DATE_TIME, "^\\d{1,2}/\\d{1,2}/\\d{4}")) %>% 
filter(Fecha == "01/10/2020")

Exploración y visualización de los datos

Se trabajó con un conjunto de datos provenientes del monitoreo de distintas variables en 494 árboles de aguacate de la Finca Guadalupe Hass, ubicada en El Tambo, Cauca. Entre las variables registradas a lo largo de varios años, para este análisis se seleccionó la temperatura correspondiente al periodo de muestreo más reciente, fechado el 01/10/2020. El resumen estadístico muestra que la temperatura media registrada durante este periodo fue de 25.8 °C con una desviación estándar de 1.8 °C. Los valores observados oscilaron entre 22.2 °C y 29.7 °C, con una mediana de 25.8 °C y un rango intercuartílico de 2.7 °C.

Resumen de variables registradas por ubicación de árboles de aguacate

Variable	Stats / Values	Freqs (% of Valid)	Graph	Missing
id_arbol [character]	1. 1 2. 10 3. 100 4. 101 5. 102 6. 103 7. 104 8. 105 9. 106 10. 107 [ 524 others ]	1 ( 0.2% ) 1 ( 0.2% ) 1 ( 0.2% ) 1 ( 0.2% ) 1 ( 0.2% ) 1 ( 0.2% ) 1 ( 0.2% ) 1 ( 0.2% ) 1 ( 0.2% ) 1 ( 0.2% ) 524 ( 98.1% )		0 (0.0%)
Latitude [numeric]	Mean (sd) : 2.4 (0) min ≤ med ≤ max: 2.4 ≤ 2.4 ≤ 2.4 IQR (CV) : 0 (0)	447 distinct values		0 (0.0%)
Longitude [numeric]	Mean (sd) : -76.7 (0) min ≤ med ≤ max: -76.7 ≤ -76.7 ≤ -76.7 IQR (CV) : 0 (0)	451 distinct values		0 (0.0%)
FORMATTED_DATE_TIME [character]	1. 01/10/2020  10:11:12 a, m	534 ( 100.0% )		0 (0.0%)
Psychro_Wet_Bulb_Temperature [numeric]	Mean (sd) : 21.6 (1.1) min ≤ med ≤ max: 19 ≤ 21.5 ≤ 24.9 IQR (CV) : 1.6 (0.1)	57 distinct values		0 (0.0%)
Station_Pressure [numeric]	Mean (sd) : 825.5 (1.1) min ≤ med ≤ max: 822.4 ≤ 825.7 ≤ 827.4 IQR (CV) : 1.7 (0)	50 distinct values		0 (0.0%)
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		0 (0.0%)
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		0 (0.0%)
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		0 (0.0%)
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		0 (0.0%)
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		0 (0.0%)
Direction_True [numeric]	Mean (sd) : 222.5 (125.4) min ≤ med ≤ max: 0 ≤ 287.5 ≤ 359 IQR (CV) : 245.8 (0.6)	198 distinct values		0 (0.0%)
Direction_Mag [numeric]	Mean (sd) : 223.3 (125) min ≤ med ≤ max: 0 ≤ 287.5 ≤ 359 IQR (CV) : 246 (0.6)	203 distinct values		0 (0.0%)
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		0 (0.0%)
Heat_Stress_Index [numeric]	Mean (sd) : 27.4 (2.6) min ≤ med ≤ max: 22.8 ≤ 27.2 ≤ 34.5 IQR (CV) : 4 (0.1)	107 distinct values		0 (0.0%)
Altitude [numeric]	Mean (sd) : 1693 (11) min ≤ med ≤ max: 1675 ≤ 1692 ≤ 1724 IQR (CV) : 15.8 (0)	47 distinct values		0 (0.0%)
Dew_Point [numeric]	Mean (sd) : 20.2 (1.1) min ≤ med ≤ max: 17.7 ≤ 20 ≤ 24 IQR (CV) : 1.4 (0.1)	58 distinct values		0 (0.0%)
Density_Altitude [numeric]	Mean (sd) : 2.6 (0.1) min ≤ med ≤ max: 2.4 ≤ 2.6 ≤ 2.7 IQR (CV) : 0.1 (0)	224 distinct values		0 (0.0%)
Wind_Chill [numeric]	Mean (sd) : 25.8 (1.8) min ≤ med ≤ max: 22.2 ≤ 25.8 ≤ 29.6 IQR (CV) : 2.6 (0.1)	74 distinct values		0 (0.0%)
Estado_Fenologico_Predominante [numeric]	Mean (sd) : 717.8 (0.5) min ≤ med ≤ max: 717 ≤ 718 ≤ 719 IQR (CV) : 1 (0)	717 : 147 ( 27.5% ) 718 : 369 ( 69.1% ) 719 : 18 ( 3.4% )		0 (0.0%)
Frutos_Afectados [numeric]	Min : 0 Mean : 0 Max : 1	0 : 531 ( 99.4% ) 1 : 3 ( 0.6% )		0 (0.0%)
Fecha [character]	1. 01/10/2020	534 ( 100.0% )		0 (0.0%)

La visualización espacial de estos datos, realizada mediante la función mapview(), permitió identificar una distribución heterogénea de la temperatura en el área de estudio. El mapa evidencia la presencia de gradientes térmicos y posibles agrupamientos de valores elevados y bajos en distintas zonas del lote, lo que justifica la exploración posterior de autocorrelación y dependencia espacial de la variable. Este análisis visual resulta fundamental para orientar el ajuste de modelos de variograma y la aplicación de técnicas de interpolación espacial.

aguacate_sf <- st_as_sf(df_aguacate, 
                        coords = c("Longitude", "Latitude"),
                        crs = 4326)

cat('<h5>`Valores de temperatura en ubicaciones de los árboles de aguacate`</h5>')

Valores de temperatura en ubicaciones de los árboles de aguacate

maptemp <-mapview(
  aguacate_sf,
  zcol = "Temperature",
  layer.name = "Temperatura [°C]",
  cex = 4,
  alpha = 0,
  col.regions = colorRampPalette(brewer.pal(9, "YlOrRd")),
  legend = TRUE
)
maptemp

Construcción y análisis del variograma muestral

La estructura espacial de la temperatura se exploró mediante la construcción de la nube de variograma y el variograma muestral, utilizando la función variogram() del paquete gstat. La nube de variograma representa la dispersión de las diferencias cuadráticas de temperatura entre todos los pares de puntos según su distancia, mientras que el variograma muestral agrupa estas diferencias en intervalos, mostrando la relación promedio entre semivarianza y distancia.

La interpretación del variograma muestral indica un incremento inicial de la semivarianza con la distancia, hasta alcanzar una meseta en torno a los 2.5 a 3.0 °C² para distancias superiores a 50 metros. El valor de nugget, que representa la variabilidad no explicada a escalas menores al muestreo o debida a error de medición, fue estimado en 1.0 a 1.5 °C². El rango espacial, definido como la distancia a partir de la cual la semivarianza se estabiliza y la autocorrelación espacial deja de ser significativa, fue identificado en aproximadamente 50 metros.

aguacate_sf <- st_transform(aguacate_sf, crs = 32618)

cat('<h5>`Nube de variograma (cloud)`</h5>')

Nube de variograma (cloud)

# Variograma cloud

plot(variogram(Temperature ~ 1, data = aguacate_sf, cloud = TRUE))

cat('<h5>`Variograma muestral`</h5>')

Variograma muestral

# Variograma muestral 

v <- variogram(Temperature ~ 1, data = aguacate_sf)
plot(v)

Ajuste y comparación de modelos de variograma

Para describir la estructura espacial observada en la temperatura, se ajustaron varios modelos teóricos de variograma utilizando la función fit.variogram() del paquete gstat. Se emplearon como parámetros iniciales los valores extraídos del variograma muestral: nugget ≈ 1.25, sill ≈ 2.75 y range ≈ 50. Los modelos evaluados incluyeron el exponencial, esférico, gaussiano, lineal, circular y pentasférico. El ajuste de cada modelo se evaluó mediante el error cuadrático de ajuste (SSErr).

cat('<h4>`Ajuste de modelos de variograma`</h4>')

Ajuste de modelos de variograma

mod <- c("Exp", "Sph", "Gau", "Lin", "Cir", "Pen")
mod_titulos <- c("Exponencial", "Esférico", "Gaussiano", "Linear", "Circular", "Pentasférico")

fits <- lapply(mod, function(m)
  fit.variogram(v, vgm(psill = 2.75, model = m, range = 50, nugget = 1.25))
)

cat('<h5>`Modelo Exponencial`</h5>')

Modelo Exponencial

plot(v, model = fits[[1]])

cat('<h5>`Modelo Esférico`</h5>')

Modelo Esférico

plot(v, model = fits[[2]])

cat('<h5>`Modelo Gaussiano`</h5>')

Modelo Gaussiano

plot(v, model = fits[[3]])

cat('<h5>`Modelo Linear`</h5>')

Modelo Linear

plot(v, model = fits[[4]])

cat('<h5>`Modelo Circular`</h5>')

Modelo Circular

plot(v, model = fits[[5]])

cat('<h5>`Modelo Pentasférico`</h5>')

Modelo Pentasférico

plot(v, model = fits[[6]])

El modelo pentasférico presentó el menor error de ajuste (SSErr = 1.03540), seguido por el esférico (SSErr = 1.05537) y el circular (SSErr = 1.07545), lo que indica que estos modelos representan mejor la estructura espacial de la temperatura en el área de estudio. Por el contrario, el modelo gaussiano presentó un error considerablemente mayor (SSErr = 3.28096), reflejando un ajuste deficiente.

Estos resultados, sumados a la inspección visual de los ajustes sobre el variograma muestral, respaldan la selección del modelo pentasférico como el más adecuado para describir la dependencia espacial de la variable temperatura en el lote evaluado.

par(mfrow = c(1, 1))
#erros
ss_err <- sapply(fits, function(fit) attr(fit, "SSErr"))

# param
params <- do.call(rbind, lapply(fits, function(x) x[2, c("psill", "range", "model")]))
nugget <- sapply(fits, function(x) x[1, "psill"])

tbl <- tibble::tibble(
  Model = mod_titulos,
  Nugget = round(nugget, 3),
  Sill   = round(params[, "psill"], 3),
  Range  = round(params[, "range"], 3),
  SSErr  = round(ss_err, 5)
)

# best model
min_row <- which.min(tbl$SSErr)
tbl$Model[min_row] <- paste0("**", tbl$Model[min_row], "**")

cat('<h5>`Resumen de ajuste de modelos de variograma`</h5>')

Resumen de ajuste de modelos de variograma

tbl %>%
  kable(align = "lrrrr") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

Model	Nugget	Sill	Range	SSErr
Exponencial	0.437	2.636	13.578	1.11254
Esférico	0.753	2.150	30.806	1.05537
Gaussiano	1.015	1.909	12.341	3.28096
Linear	0.922	1.979	24.593	1.36529
Circular	0.792	2.103	27.218	1.07545
Pentasférico	0.725	2.194	37.524	1.03540

Interpolación espacial mediante Kriging

La interpolación de la temperatura en áreas no muestreadas se realizó mediante el método de kriging, utilizando el modelo de variograma pentasférico previamente ajustado. Este procedimiento se implementó con la función gstat::gstat() para definir el objeto de kriging, seguido de predict() para obtener las predicciones sobre una grilla regular abarcando el área de muestreo. Las predicciones resultantes fueron posteriormente recortadas al polígono definido por el convex hull de los puntos observados, con el fin de evitar extrapolaciones fuera del dominio real de los datos.

# Kriging
fit_pen <- fit.variogram(v, model = vgm(psill = 2.75, 
            model = "Pen", range = 50, nugget = 1.25))
k <- gstat(formula = Temperature ~ 1, data = aguacate_sf, model = fit_pen)

# grid
bbox <- st_bbox(aguacate_sf)
grd <- expand.grid(
  x = seq(bbox["xmin"], bbox["xmax"], length.out = 150),
  y = seq(bbox["ymin"], bbox["ymax"], length.out = 150)
)
grd_sf <- st_as_sf(grd, coords = c("x", "y"), crs = st_crs(aguacate_sf))

# prediction
kpred <- predict(k, newdata = grd_sf)

La comparación entre el mapa de temperaturas observadas y el mapa de predicción obtenida por kriging, visualizados de manera sincronizada con leafsync::sync(), evidencia una notable similitud en los patrones espaciales identificados. Ambas visualizaciones muestran la presencia de gradientes térmicos y agrupamientos de valores altos y bajos en regiones coincidentes dentro del lote, lo que indica que el modelo de kriging ha capturado adecuadamente la estructura espacial de la variable temperatura. Esto evidencia que el modelo ajustado logra capturar adecuadamente la dependencia espacial de la variable. Además, la calidad predictiva del modelo se evaluó mediante validación cruzada tipo leave-one-out, obteniendo un error cuadrático medio (RMSE) de 1.07 °C, un error absoluto medio (MAE) de 0.84 °C y un error medio (ME) cercano a cero (-0.006 °C). Estos resultados cuantitativos respaldan la robustez y ausencia de sesgo relevante en las predicciones generadas

# convex hull
hull <- st_convex_hull(st_union(aguacate_sf))
kpred_cropped <- kpred[st_within(st_geometry(kpred), hull, sparse = FALSE), ]

# map pred
map_krig <- mapview(
  kpred_cropped,
  zcol = "var1.pred",
  layer.name = "Predicción kriging [°C]",
  col.regions = colorRampPalette(brewer.pal(9, "YlOrRd")),
  legend = TRUE,
  alpha = 0
)

# maps
cat('<h5>`Comparación: temperaturas observadas y predicción espacial por kriging`</h5>')

Comparación: temperaturas observadas y predicción espacial por kriging

leafsync::sync(maptemp, map_krig)

# leave-one-out (LOOCV) 
cv <- krige.cv(Temperature ~ 1, locations = aguacate_sf, model = fit_pen, nfold = nrow(aguacate_sf))

# tbl
cv_summary <- tibble::tibble(
  RMSE = sqrt(mean(cv$residual^2, na.rm = TRUE)),
  MAE  = mean(abs(cv$residual), na.rm = TRUE),
  ME   = mean(cv$residual, na.rm = TRUE)
)

cat('<h5>`Resultados de validación cruzada (LOOCV) para kriging`</h5>')

Resultados de validación cruzada (LOOCV) para kriging

cv_summary %>%
  kable(align = "rrr") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

RMSE	MAE	ME
1.070817	0.8444086	-0.0058144

Conclusiones

Este análisis demostró la utilidad de las herramientas geoestadísticas para modelar y predecir la variabilidad espacial de la temperatura en plantaciones de aguacate Hass, utilizando datos primarios recolectados en la Finca Guadalupe Hass de El Tambo, Cauca. La metodología aplicada, desde la exploración de datos, el ajuste y selección de modelos de variograma, hasta la interpolación espacial mediante kriging y su validación, permitió obtener resultados sólidos, mostrando patrones espaciales coherentes y estimaciones precisas a escala de lote.

Estos hallazgos adquieren especial relevancia en el contexto de El Tambo, un municipio destacado en la producción de aguacate Hass gracias a sus condiciones de altitud y clima favorables. Además, el desarrollo sostenible y la transformación productiva en la región se fortalecen con la participación de familias campesinas, indígenas y afrocolombianas, así como la implementación de centros de acopio y buenas prácticas agrícolas. El estudio evidencia el valor de integrar datos espaciales y metodologías avanzadas para apoyar la toma de decisiones en territorios rurales en procesos de desarrollo y transición.