1 Introducción

El análisis geoestadístico aplicado a datos climáticos constituye una herramienta fundamental para comprender la variabilidad espacial de las condiciones ambientales en sistemas agrícolas. En el presente estudio se trabaja con información recolectada en una finca de aguacates ubicada en el departamento del Cauca, donde se registraron variables de altura, latitud y longitud asociadas a cada árbol. El objetivo es evaluar la dependencia espacial de la altura mediante el uso de semivariogramas, seleccionar el modelo teórico que mejor represente la estructura de autocorrelación y realizar una interpolación espacial con la metodología de kriging.

2 Importar datos y librerias

# Librerias a utilizar
library(readxl)
library(geoR)
library(raster)
library(ggplot2)
library(leaflet)
require(rasterVis)

#Cargue de la base de datos
aguacate= read_excel("Datos_Completos_Aguacate.xlsx")

# Filtro de columnas a trabajar
aguacate_filtrado <- subset(aguacate, select= c(Latitude, Longitude, Altitude))
aguacate_filtrado

2.1 Ubicación de la finca caso de estudio

leaflet() %>% addTiles() %>% addCircleMarkers(lng = aguacate_filtrado$Longitude,lat = aguacate_filtrado$Latitude,radius = 0.2,color = "black")

3 Geoestadística

3.1 Creación de la variable regionalizada

geo_aguacate=as.geodata(aguacate_filtrado, coords.col = 1:2, data.col = 3)
plot(geo_aguacate)

La imagen representa el análisis exploratorio de los datos de altura registrados en la finca de aguacates, mostrando tanto su distribución espacial como su comportamiento estadístico. En los gráficos de dispersión se observa que las mediciones se encuentran bien distribuidas en el espacio, siguiendo un patrón regular que sugiere una toma sistemática de datos dentro del terreno. La altura no presenta una tendencia lineal evidente con las coordenadas geográficas. El histograma con la curva de densidad revela una distribución asimétrica hacia la derecha, ccon mayor concentración de valores entre 1680 y 1700 metros sobre el nivel del mar (MSNM), rango que podría considerarse óptimo para el cultivo.

3.2 Estadísticas

estadistica=summary(dist(geo_aguacate$coords))
estadistica
##      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

3.3 Semivariograma

# Semivariograma
variograma <- variog(geo_aguacate, option="bin", uvec=seq(0, 0.0009178, 0.00009178))
## variog: computing omnidirectional variogram
plot(variograma)

variograma_mc=variog.mc.env(geo_aguacate,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 la semivarianza aumenta rápidamente al inicio y luego se estabiliza alrededor de un valor cercano a 3, lo que indica que existe autocorrelación espacial fuerte a distancias cortas: los puntos cercanos presentan alturas similares. A medida que la distancia crece, las diferencias dejan de aumentar significativamente, lo que sugiere que más allá de cierto rango las mediciones son prácticamente independientes.

Las líneas superior e inferior corresponden al intervalo de confianza generado por simulaciones Monte Carlo; el hecho de que los puntos del semivariograma se mantengan dentro de ese rango confirma que la estructura espacial observada es consistente y no producto del azar.

4 Ajuste del modelo teorico

# Ajuste del modelo teórico
ini.vals <- expand.grid(
  sill = seq(100, 140, length.out = 10),   # valores alrededor del sill observado
  phi  = seq(0.0005, 0.0007, length.out = 10) # valores alrededor del rango observado
)
model_gaus <- variofit(variograma, ini=ini.vals, cov.model="gaussian", wei="npair", min="optim")
## 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 "140"   "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 44622020.9226709
model_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 "140"   "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 10902702.8852714
model_sph  <- variofit(variograma, ini=ini.vals, cov.model="spheric", 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 "113.33" "0"   "0"   "0.5"
## status        "est"    "est" "est" "fix"
## loss value: 10002366.6387751
# Primero graficar el semivariograma
plot(variograma)

# Luego añadir las líneas de los modelos
lines(model_gaus, col="red")
lines(model_exp, col="blue")
lines(model_sph, col="purple")

El gráfico muestra el semivariograma experimental junto con los tres modelos teóricos ajustados: gaussiano (línea roja), exponencial (línea azul) y esférico (línea morada). Se observa que la semivarianza aumenta rápidamente con la distancia, reflejando una autocorrelación espacial fuerte entre puntos cercanos, y luego se estabiliza alrededor de un valor cercano a 120, que corresponde al sill o varianza total del fenómeno.

El modelo gaussiano (rojo) describe una transición suave y continua, ajustándose mejor a los puntos experimentales en las distancias cortas y medias, lo que sugiere una variación gradual de la altura en el espacio. El modelo exponencial (azul) muestra un incremento más lineal, mientras que el esférico (morado) alcanza el sill de forma más abrupta. En conjunto, el comportamiento del semivariograma confirma la existencia de una estructura espacial bien definida y una correlación significativa hasta aproximadamente 0.0006 unidades de distancia, validando el uso del modelo gaussiano para la interpolación por kriging en el área de estudio.

5 Interpolación

min(aguacate_filtrado[,1])
## [1] 2.392101
max(aguacate_filtrado[,1])
## [1] 2.393634
min(aguacate_filtrado[,2])
## [1] -76.7118
max(aguacate_filtrado[,2])
## [1] -76.71022
geodatos_grid2 = expand.grid(
  lon = seq(min(geo_aguacate$coords[,2]), max(geo_aguacate$coords[,2]), l = 100),
  lat = seq(min(geo_aguacate$coords[,1]), max(geo_aguacate$coords[,1]), l = 100)
)
plot(geodatos_grid2)
points(aguacate_filtrado[,2:1],col="red")

La figura representa la malla de interpolación espacial generada a partir de las coordenadas geográficas de los puntos de muestreo, donde cada celda del grid define una ubicación potencial para la estimación mediante kriging. Los puntos rojos corresponden a las observaciones empíricas de la variable analizada, distribuidas sobre un plano de coordenadas geográficas (latitud y longitud). La disposición regular del entramado evidencia la homogeneidad del espacio de predicción, mientras que la concentración y alineación de los puntos sugieren una estructura espacial coherente con el patrón de muestreo en campo.

geodatos_grid = expand.grid(
  lon = seq(min(geo_aguacate$coords[,1]), max(geo_aguacate$coords[,1]), l = 100),
  lat = seq(min(geo_aguacate$coords[,2]), max(geo_aguacate$coords[,2]), l = 100)
)
geodatos_ko=krige.conv(geo_aguacate, loc=geodatos_grid,
      krige= krige.control(nugget=0.2,trend.d="cte", 
      trend.l="cte",cov.pars=c(sigmasq=120, phi=0.0006 )))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
par(mfrow=c(1,2))
image(geodatos_ko, main="kriging Predict", xlab="East", ylab="North")
contour(geodatos_ko, main = "Kriging Contours", xlab = "East", ylab = "North",
        drawlabels = TRUE)

La figura presenta los resultados del proceso de interpolación espacial mediante el método de kriging ordinario, aplicado sobre la variable regionalizada correspondiente a la altura en la finca de estudio. En el panel izquierdo, titulado Kriging Predict, se observa la superficie continua de predicción, donde los tonos cromáticos reflejan la variación espacial de los valores estimados. Las zonas de color más oscuro indican áreas con valores superiores, mientras que las tonalidades claras representan regiones con menor intensidad del fenómeno. Este patrón evidencia una gradiente espacial bien definida, coherente con la estructura de autocorrelación identificada en el semivariograma.

El panel derecho, Kriging Contours, muestra las líneas de contorno derivadas de la misma interpolación, las cuales permiten visualizar con mayor precisión las transiciones entre intervalos de valores. La disposición de las curvas revela una continuidad espacial suave y una tendencia diagonal que coincide con la orientación del muestreo original.

par(mfrow=c(1,2))
image(geodatos_ko, main="kriging StDv Predicted",val=sqrt(geodatos_ko$krige.var), xlab="East", ylab="North")
contour(geodatos_ko,main="kriging StDv Predict",val=sqrt(geodatos_ko$krige.var), drawlabels=TRUE)

La figura ilustra la distribución espacial de la desviación estándar de las predicciones obtenidas mediante kriging ordinario, la cual refleja el grado de incertidumbre asociado a cada estimación. En el panel izquierdo (Kriging StDv Predicted), se observa una superficie continua donde los tonos oscuros indican zonas de mayor error de predicción, mientras que las áreas claras representan regiones con menor incertidumbre. Este patrón evidencia que la precisión del modelo es más alta en las zonas con mayor densidad de puntos de muestreo y disminuye hacia los extremos del área de interpolación, donde la información empírica es escasa.

El panel derecho (Kriging StDv Predict) complementa la interpretación mediante líneas de contorno, que delimitan intervalos de variabilidad del error. La disposición de las curvas muestra una estructura espacial coherente con la distribución de los datos originales, confirmando que la incertidumbre se incrementa gradualmente conforme aumenta la distancia entre puntos observados

6 Transformación a raster

pred=cbind(geodatos_grid2,geodatos_ko$predict)
temp_predict=rasterFromXYZ(cbind(geodatos_grid2,geodatos_ko$predict))
plot(temp_predict)

levelplot(temp_predict,par.settings =BuRdTheme)

temp_error=rasterFromXYZ(cbind(geodatos_grid2,sqrt(geodatos_ko$krige.var)))
levelplot(temp_error,par.settings =BuRdTheme)

La transformación de las predicciones de kriging a formato raster constituye un paso metodológico esencial para la representación cartográfica y el análisis espacial avanzado. Al convertir los resultados en objetos raster, se obtiene una estructura matricial que facilita la visualización continua de la variable regionalizada y permite integrar la información en sistemas de información geográfica (SIG). Esta operación garantiza que cada celda del raster contenga un valor estimado de la altura, lo cual posibilita la generación de mapas temáticos de alta resolución y la comparación con otras capas ambientales o productivas.

7 Conclusiones

El análisis geoestadístico realizado permitió caracterizar de manera precisa la estructura espacial de la altura de los árboles de aguacate en la finca objeto de estudio. El semivariograma experimental evidenció una autocorrelación espacial significativa a distancias cortas, confirmando que las observaciones próximas presentan valores similares y que la variabilidad se estabiliza a partir de un rango aproximado de 0.0006 unidades geográficas. Este comportamiento valida la existencia de una dependencia espacial real y justifica el uso de modelos geoestadísticos para la estimación de valores no observados.

El ajuste de los modelos teóricos mostró que el modelo gaussiano ofrece el mejor desempeño, al reproducir con mayor fidelidad la tendencia del semivariograma experimental y reflejar una transición suave entre las distancias cortas y medias. Dicho modelo fue empleado en la interpolación mediante kriging ordinario, generando una superficie continua que representa la distribución espacial de la altura con alta coherencia respecto al patrón de muestreo.

Los mapas de predicción y de desviación estándar revelaron una gradiente espacial bien definida, con zonas de mayor altura concentradas en el sector sur de la finca y una disminución progresiva hacia el norte. Asimismo, el análisis de incertidumbre mostró que la precisión del modelo es mayor en las áreas con mayor densidad de datos, mientras que la variabilidad aumenta en los extremos del área de interpolación, donde la información empírica es limitada.