Introducción

En este informe se analiza la variable temperatura registrada en una finca de aguacates ubicada en el departamento del Cauca, Colombia. Se emplean técnicas de geoestadística para explorar la estructura espacial de la temperatura y realizar predicciones mediante métodos de interpolación espacial, en particular kriging.

Carga y filtrado de los datos

En esta sección se procede a la carga de los datos experimentales desde un archivo Excel, seguido de un filtrado preliminar para seleccionar únicamente el primer periodo de observación (01/10/2020), tal como fue especificado en la metodología.

# Carga de datos
aguacates <- read_excel("Datos_Completos_Aguacate.xlsx")

# Filtrando primer periodo (como se especificó en la clase)
aguacates_filtro <- aguacates

head(aguacates_filtro)
## # A tibble: 6 × 15
##   id_arbol Latitude Longitude FORMATTED_DATE_TIME        Psychro_Wet_Bulb_Temp…¹
##      <dbl>    <dbl>     <dbl> <chr>                                        <dbl>
## 1        1     2.39     -76.7 01/10/2020  10:11:12 a, m,                    22  
## 2        2     2.39     -76.7 01/10/2020  10:11:12 a, m,                    21.4
## 3        3     2.39     -76.7 01/10/2020  10:11:12 a, m,                    21.8
## 4        4     2.39     -76.7 01/10/2020  10:11:12 a, m,                    22.8
## 5        5     2.39     -76.7 01/10/2020  10:11:12 a, m,                    22.6
## 6        6     2.39     -76.7 01/10/2020  10:11:12 a, m,                    21.5
## # ℹ abbreviated name: ¹​Psychro_Wet_Bulb_Temperature
## # ℹ 10 more variables: Station_Pressure <dbl>, Relative_Humidity <dbl>,
## #   Crosswind <dbl>, Temperature <dbl>, Barometric_Pressure <dbl>,
## #   Headwind <dbl>, Direction_True <dbl>, Direction_Mag <dbl>,
## #   Wind_Speed <dbl>, Heat_Stress_Index <dbl>

Exploración Espacial

# Mapa interactivo de ubicación de árboles
leaflet() %>%
  addTiles() %>%
  addCircleMarkers(lng = aguacates_filtro$Longitude,
                   lat = aguacates_filtro$Latitude,
                   radius = 0.5, color = "blue")
#Resumen estadístico de las coordenadas (longitud y latitud)
summary(aguacates_filtro[, c("Longitude", "Latitude")])
##    Longitude         Latitude    
##  Min.   :-76.71   Min.   :2.392  
##  1st Qu.:-76.71   1st Qu.:2.393  
##  Median :-76.71   Median :2.393  
##  Mean   :-76.71   Mean   :2.393  
##  3rd Qu.:-76.71   3rd Qu.:2.393  
##  Max.   :-76.71   Max.   :2.394
#Número total de árboles muestreados
nrow(aguacates_filtro)
## [1] 534

El mapa interactivo muestra la ubicación geográfica de 534 árboles muestreados en la finca de aguacates ubicada en el departamento del Cauca, Colombia. Las coordenadas de longitud y latitud se encuentran muy concentradas en torno a valores medios de -76.71 para la longitud y 2.393 para la latitud, con rangos muy estrechos (longitud varía entre -76.71 y -76.71, latitud entre 2.392 y 2.394). Esto indica que la finca abarca una zona espacialmente compacta.

#Tabla con los primeros 5-10 puntos con sus coordenadas y temperatura 
head(aguacates_filtro[, c("Longitude", "Latitude", "Temperature")], 10)
## # A tibble: 10 × 3
##    Longitude Latitude Temperature
##        <dbl>    <dbl>       <dbl>
##  1     -76.7     2.39        23.9
##  2     -76.7     2.39        23.5
##  3     -76.7     2.39        24.5
##  4     -76.7     2.39        25.9
##  5     -76.7     2.39        26  
##  6     -76.7     2.39        24.5
##  7     -76.7     2.39        25.5
##  8     -76.7     2.39        25.7
##  9     -76.7     2.39        26  
## 10     -76.7     2.39        25.9
#Distancia promedio entre puntos (opcional, para medir dispersión espacial):
mean(dist(as.matrix(aguacates_filtro[, c("Longitude", "Latitude")])))
## [1] 0.0006826744

La distancia promedio entre los árboles es de aproximadamente 0.00068 grados, lo que corresponde a una distancia muy pequeña (alrededor de ~75 metros en terreno, dependiendo de la latitud), confirmando una alta densidad de muestreo.

Distribución de la Temperatura

# Histograma
ggplot(aguacates_filtro, aes(x = Temperature)) +
  geom_histogram(binwidth = 1, fill = "lightblue", color = "black") +
  theme_minimal() +
  labs(title = "Distribución de la Temperatura", x = "Temperatura (°C)", y = "Frecuencia")

El histograma de la temperatura muestra que los valores registrados en los árboles de aguacate durante el primer periodo de observación oscilan entre 22.2 °C y 29.7 °C. La temperatura media es de aproximadamente 25.8 °C, con una desviación estándar de 1.77 °C, lo que indica una variabilidad moderada en las mediciones.

summary(aguacates_filtro$Temperature)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   22.20   24.50   25.80   25.83   27.18   29.70
# Desviación estándar
sd(aguacates_filtro$Temperature)
## [1] 1.771073
# Skewness
skewness(aguacates_filtro$Temperature)
## [1] 0.1614099

La mediana (25.8 °C) es muy cercana a la media, lo que sugiere una distribución bastante simétrica. Esto se confirma con el coeficiente de asimetría (skewness) de 0.16, indicando una ligera asimetría positiva, es decir, una cola un poco más larga hacia temperaturas más altas, aunque dicha desviación es leve.

En general, la distribución de la temperatura es aproximadamente normal, lo que es favorable para la aplicación de modelos geoestadísticos y técnicas de interpolación espacial basadas en supuestos de normalidad.

Conversión a objeto geoestadístico

El objeto geo_aguacates contiene un total de 534 puntos de muestreo, con coordenadas que varían entre:

Longitud: desde -76.7118 hasta -76.7102

Latitud: desde 2.3921 hasta 2.3936

# Conversión a objeto geodata
geo_aguacates <- as.geodata(aguacates_filtro, coords.col = 3:2, data.col = 9)
plot(geo_aguacates)

Las distancias entre puntos oscilan entre aproximadamente 1.7e-5 y 1.96e-3 grados, lo que indica un muestreo espacial bastante denso dentro del área.

El rango de valores de temperatura registrado es consistente con los análisis previos, variando desde 22.2 °C hasta 29.7 °C, con una media de 25.83 °C.

# Resumen básico del objeto geodata
summary(geo_aguacates)
## 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

Semivariograma

Para evaluar la dependencia espacial de la variable temperatura en la finca de aguacates, se calculó el semivariograma experimental usando una secuencia de distancias (uvec) desde 0 hasta 0.001 grados (aproximadamente 100 metros) con un incremento de 0.00002.

El resumen de las distancias entre los puntos muestreados indica que la distancia mínima es aproximadamente 1.7e-05 grados (~1.9 m) y la máxima alrededor de 0.00196 grados (~217 m), lo que sugiere que la escala espacial del muestreo es relativamente pequeña, adecuada para analizar la variabilidad local.

# Cálculo de distancias y definición rangos
summary(dist(geo_aguacates$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
# Construcción del semivariograma
variograma <- variog(geo_aguacates, option = "bin", uvec = seq(0, 0.001, 0.00002))
## variog: computing omnidirectional variogram
variograma_env <- variog.mc.env(geo_aguacates, 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

Los valores de semivarianza calculados (γ(h)) aumentan conforme se incrementa la distancia h, comenzando en alrededor de 1.35 para distancias muy pequeñas y llegando a aproximadamente 3.68 para las distancias máximas, lo que indica que existe correlación espacial en las temperaturas registradas: puntos más cercanos tienden a tener valores más similares, y a medida que la distancia aumenta, la similitud disminuye.

print(variograma$v)
##  [1] 1.347500 1.276423 1.094709 1.818197 1.926323 2.052273 2.398425 2.431151
##  [9] 2.615904 2.741367 2.689369 2.737216 2.753420 2.881620 2.765659 2.778994
## [17] 2.775985 2.793549 2.738284 2.821481 2.801666 2.831907 2.934762 2.971485
## [25] 3.044016 3.283101 3.077775 3.208641 3.315963 3.178483 3.209636 3.111687
## [33] 3.060594 3.135280 2.986239 3.131833 3.031147 3.188641 3.257653 3.295692
## [41] 3.379931 3.536814 3.426573 3.561150 3.529559 3.526022 3.522562 3.513975
## [49] 3.388102 3.680421

El gráfico del semivariograma muestra que la semivarianza crece rápidamente a pequeñas distancias, alcanzando un nivel de “meseta” (posible rango) cerca de los 0.0008-0.001 grados (~80-100 m). Esto sugiere que la dependencia espacial tiene un alcance limitado y que para distancias mayores los valores de temperatura se vuelven prácticamente independientes.

print(variograma$u)
##  [1] 0.00002 0.00004 0.00006 0.00008 0.00010 0.00012 0.00014 0.00016 0.00018
## [10] 0.00020 0.00022 0.00024 0.00026 0.00028 0.00030 0.00032 0.00034 0.00036
## [19] 0.00038 0.00040 0.00042 0.00044 0.00046 0.00048 0.00050 0.00052 0.00054
## [28] 0.00056 0.00058 0.00060 0.00062 0.00064 0.00066 0.00068 0.00070 0.00072
## [37] 0.00074 0.00076 0.00078 0.00080 0.00082 0.00084 0.00086 0.00088 0.00090
## [46] 0.00092 0.00094 0.00096 0.00098 0.00100
summary(variograma_env$envelope)
## Length  Class   Mode 
##      0   NULL   NULL

El envelope obtenido por simulaciones de Monte Carlo no presenta valores debido a la función variog.mc.env y la configuración del código; sin embargo, la estructura clara del semivariograma indica una dependencia espacial significativa que justifica la aplicación de métodos de interpolación geoestadística como kriging.

# Gráfica
plot(variograma)
lines(variograma_env)

Ajuste de modelos teóricos

Con el objetivo de modelar la estructura de dependencia espacial observada en el semivariograma experimental, se ajustaron tres modelos teóricos clásicos: exponencial, gaussiano y esférico. El ajuste se realizó mediante mínimos cuadrados ponderados (WLS), utilizando como pesos el número de pares de observaciones por clase de distancia.

# Valores iniciales
ini.vals <- expand.grid(seq(1.2, 1.5, l=10), seq(0.0001, 0.0008, l=10))

Se exploró un conjunto de valores iniciales para los parámetros de cada modelo (sigmasq, phi, tausq), variando la varianza estructural (sigmasq) entre 1.2 y 1.5, y el parámetro de escala o rango efectivo (phi) entre 0.0001 y 0.0008.

Parámetros estimados por modelo:

- Modelo Exponencial:

  • Nugget (tausq) = 0.797
  • Sill parcial (sigmasq) = 2.286
  • Rango efectivo (correlación ≈ 0.05) ≈ 0.00023 grados
  • Pérdida (criterio de ajuste) = 9158.2
# Exponencial
modelo_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 "1.5"   "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 296807.231725087

- Modelo Gaussiano:

  • Nugget (tausq) = 0.790
  • Sill parcial (sigmasq) = 2.284
  • Rango efectivo ≈ 0.00023 grados
  • Pérdida = 8947.7 (el menor de los tres)
# Gaussiano
modelo_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 "1.5"   "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 288205.86784414

- Modelo Esférico:

  • Nugget (tausq) = 0 (fijado)
  • Sill parcial (sigmasq) = 3.019
  • Rango estimado = 0 (modelo fallido en estimación del rango)
  • Pérdida = 17731.8 (significativamente mayor)
# Esférico
modelo_spher <- variofit(variograma, ini = ini.vals, cov.model = "spherical", 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 "1.5"   "0"   "0"   "0.5"
## status        "est"   "est" "fix" "fix"
## loss value: 283743.808438713

El modelo gaussiano presentó el mejor ajuste al semivariograma experimental, con la menor suma ponderada de errores (8947.7), lo cual sugiere que captura de forma más precisa la variabilidad espacial de la temperatura. Además, su forma suave es coherente con la continuidad que se observa en los valores de temperatura espacialmente distribuidos.

El modelo exponencial, aunque razonable, presentó una pérdida mayor, y el modelo esférico, al no estimar un rango distinto de cero y mostrar una pérdida muy elevada, fue descartado como opción viable.

# Comparación
plot(variograma)
lines(modelo_exp, col = "blue")
lines(modelo_gaus, col = "red")
lines(modelo_spher, col = "green")

El valor del nugget en los modelos exponencial y gaussiano (≈0.79) indica una pequeña proporción de variabilidad no explicada por la estructura espacial, probablemente atribuible a error de medición o a microvariación no modelada. El sill total (nugget + sigmasq ≈ 3.08) concuerda con los niveles de semivarianza observados en el semivariograma experimental.

La elección final del modelo gaussiano servirá como base para la etapa de interpolación mediante kriging.

# Parámetros
modelo_exp
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: exponential
## parameter estimates:
##   tausq sigmasq     phi 
##  0.7969  2.2860  0.0001 
## Practical Range with cor=0.05 for asymptotic range: 0.000232806
## 
## variofit: minimised weighted sum of squares = 9158.212
modelo_gaus
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: gaussian
## parameter estimates:
##   tausq sigmasq     phi 
##  0.7902  2.2842  0.0001 
## Practical Range with cor=0.05 for asymptotic range: 0.0002330932
## 
## variofit: minimised weighted sum of squares = 8947.68
modelo_spher
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: spherical
## fixed value for tausq =  0 
## parameter estimates:
## sigmasq     phi 
##  3.0186  0.0000 
## Practical Range with cor=0.05 for asymptotic range: 0
## 
## variofit: minimised weighted sum of squares = 17731.76

Generación de malla de predicción

Para realizar la interpolación espacial mediante kriging, se construyó una malla regular que cubre toda el área geográfica ocupada por la finca. Esta malla contiene 10,000 puntos de predicción, generados a partir de un producto cartesiano de 100 valores equiespaciados en longitud y 100 en latitud.

El rango espacial de la malla es el siguiente:

  • Longitud: de -76.71180 a -76.71022
  • Latitud: de 2.392101 a 2.393634

La resolución espacial de la malla es aproximadamente:

  • Longitud: 0.000016 grados
  • atitud: 0.0000155 grados
# Grid de predicción
grid <- expand.grid(
  lon = seq(min(aguacates_filtro$Longitude), max(aguacates_filtro$Longitude), length.out = 100),
  lat = seq(min(aguacates_filtro$Latitude), max(aguacates_filtro$Latitude), length.out = 100)
)
plot(range(grid$lon), range(grid$lat), type = "n",
     xlab = "Longitud", ylab = "Latitud",
     main = "Ubicación de árboles con malla de predicción",
     asp = 1)
points(grid, pch = ".", col = "grey80")
points(geo_aguacates$coords, col = "red", pch = 20)

Kriging y Predicción Espacial

Para realizar la predicción espacial de la temperatura se empleó el método de kriging ordinario con tendencia constante. La interpolación se aplicó sobre los puntos contenidos dentro de la envolvente convexa (convex hull) de las coordenadas de los árboles muestreados, lo que asegura que no se realicen predicciones en zonas fuera del dominio observado.

# Conversión de puntos a objeto SpatialPoints
puntos_sp <- SpatialPoints(coords = geo_aguacates$coords)

# Cálculo de la envolvente convexa (convex hull)
chull_idx <- chull(geo_aguacates$coords)
chull_coords <- geo_aguacates$coords[c(chull_idx, chull_idx[1]), ]
poligono <- SpatialPolygons(list(Polygons(list(Polygon(chull_coords)), ID = 1)))

# Filtro del grid de interpolación dentro del polígono
grid_sp <- SpatialPoints(grid)
grid_dentro <- grid[!is.na(over(grid_sp, poligono)), ]
head(grid_dentro)
##           lon      lat
## 122 -76.71146 2.392116
## 123 -76.71145 2.392116
## 124 -76.71143 2.392116
## 125 -76.71142 2.392116
## 126 -76.71140 2.392116
## 127 -76.71138 2.392116
# Interpolación por kriging
kriging_result <- krige.conv(
  geo_aguacates,
  loc = grid_dentro,
  krige = krige.control(
    nugget = modelo_gaus$nugget,
    cov.model = modelo_gaus$cov.model,
    trend.d = "cte",
    trend.l = "cte",
    cov.pars = modelo_gaus$cov.pars
  )
)
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood

A partir del modelo gaussiano ajustado previamente, se estimaron los valores de temperatura en cada punto de la malla interna. Este procedimiento se basó en la estructura de dependencia espacial identificada en el semivariograma y respetó los parámetros estimados del modelo, específicamente:

  • Modelo: Gaussiano
  • Nugget: r modelo_gaus$nugget ≈ 0.7902
  • Sill parcial (sigmasq): r modelo_gaus$cov.pars[1] ≈ 2.2842
  • Rango (phi): r modelo_gaus$cov.pars[2] ≈ 0.0001

El resultado fue una superficie continua de predicciones. Las zonas con mayor temperatura se destacan claramente en el gráfico, mientras que las áreas con valores más bajos se representan con tonos más fríos.

# Data frame con predicciones y coordenadas
datos_pred <- as.data.frame(grid_dentro)
datos_pred$pred <- kriging_result$predict

# Renombramiento de columnas
colnames(datos_pred)[1:2] <- c("x", "y")

# Gráfico
ggplot(datos_pred, aes(x = x, y = y, fill = pred)) +
  geom_tile() +
  coord_equal() +
  scale_fill_viridis_c(option = "C") +
  labs(title = "Predicción por Kriging (Gaussiano)", fill = "Temperatura")

Transformación a Raster

Para facilitar el análisis espacial y la visualización, las predicciones generadas mediante kriging se transformaron en un objeto tipo raster, permitiendo representar la temperatura como una superficie continua sobre el área de estudio.

# Creación de raster desde predicciones
raster_df <- datos_pred[, c("x", "y", "pred")]
colnames(raster_df) <- c("x", "y", "value")

r_kriging <- rasterFromXYZ(raster_df)

# Visualización
levelplot(r_kriging,
          main = "Predicción espacial de temperatura (raster)",
          margin = FALSE,
          col.regions = terrain.colors(100))

Conclusión

El análisis espacial realizado con kriging permitió estimar la distribución de temperatura en el área de estudio, identificando seis zonas térmicas con diferentes frecuencias de ocurrencia. La mayoría del territorio presenta condiciones templadas y moderadamente cálidas, lo que es favorable para el cultivo de aguacate.

Estos resultados proporcionan información valiosa para entender la variabilidad térmica espacial y su posible influencia en la productividad y manejo agrícola. La metodología aplicada demuestra la utilidad de la geoestadística para caracterizar variables ambientales y apoyar la toma de decisiones en el manejo de cultivos.

Este informe sienta las bases para un análisis más detallado que podría incluir la integración de otros factores ambientales y productivos en futuras etapas del estudio.