options(encoding = "UTF-8")
library(readxl)
library(dplyr)
library(stringr)
library(tidyr)
library(sf)
library(geoR)
library(ggplot2)
library(knitr)
library(robustbase)
library(diptest)

Introducción

El presente trabajo aplica técnicas de interpolación espacial mediante metodología kriginga usando datos meteorológicos registrados por sensores portatiles tomados en árboles de aguacate en el departamento del Cauca, Colombia - municipío de Timbio. Las mediciones corresponden al 1 de octubre de 2020 y provienen de 534 puntos de observación distribuidos en el área de estudio. Las variables analizadas son temperatura del aire,humedad relativa y punto de rocío.

El proyecto comprende tres etapas principales: la descripción de los datos, la estimación del semi- variograma empírico para cada variable y el ajuste de modelos teóricos de variograma con la predicción mediante la metodología kriging. Los resultados se presentan como mapas de superficies, permitiendo evaluar tanto la distribución espacial de cada variable.

Área de estudio

El área de estudio corresponde a una parcela de cultivo de aguacate ubicada en el municipio de Timbío, departamento del Cauca, Colombia, ( Latitud: 2.393° , Longitud: 76.711° ). En area se distribuyen 534 árboles con sensores portatiles, cuyas mediciones meteorologicas, correspondientes al 1 de octubre de 2020.

library(readxl)
library(dplyr)
library(stringr)
library(sf)
library(leaflet)

aguacate <- read_excel("data/aguacate.xlsx")
aguacate_dia <- aguacate %>%
  filter(str_starts(FORMATTED_DATE_TIME, "01/10/2020"))
pts_wgs <- st_as_sf(aguacate_dia,
                    coords = c("Longitude", "Latitude"),
                    crs = 4326)
mapa_df <- pts_wgs %>%
  mutate(lon = st_coordinates(.)[, 1],
         lat = st_coordinates(.)[, 2]) %>%
  st_drop_geometry()
pal <- colorNumeric(palette = "Spectral", domain = mapa_df$Temperature)
leaflet(mapa_df) %>%
  addProviderTiles(providers$Esri.WorldImagery, group = "Satélite") %>%
  addProviderTiles(providers$OpenStreetMap,     group = "Calles") %>%
  addCircleMarkers(
    lng = ~lon, lat = ~lat,
    radius = 5,
    color = ~pal(Temperature),
    stroke = TRUE, weight = 1, fillOpacity = 0.8,
    popup = ~paste0(
      "Árbol: ", id_arbol, "<br>",         
      "Temp: ", round(Temperature, 1), " °C<br>",
      "Humedad: ", round(Relative_Humidity, 1), " %<br>",
      "Estrés calórico: ", round(Heat_Stress_Index, 1), "<br>",
      "Altitud: ", round(Altitude), " m"
    )
  ) %>%
  addLegend("bottomright", pal = pal, values = ~Temperature,
            title = "Temp (°C)", opacity = 0.9) %>%
  addLayersControl(
    baseGroups = c("Satélite", "Calles"),
    options = layersControlOptions(collapsed = FALSE)
  )

Descripción de los datos

library(dplyr)
library(tidyr)
library(knitr)
library(kableExtra)
library(ggplot2)
vars <- c("Temperature", "Relative_Humidity",
          "Dew_Point")

Tabla 1. Estadísticos descriptivos de las variables (Fecha : 1 oct 2020)
Variable	n	Media	DE	CV_pct	Min	Mediana	Max	Rango
Dew_Point	534	20.16	1.13	5.62	17.7	20.00	24.0	6.3
Relative_Humidity	534	71.17	5.34	7.51	59.5	71.25	91.6	32.1
Temperature	534	25.83	1.77	6.86	22.2	25.80	29.7	7.5

La Tabla 1 resume las estadísticas descriptivas de las variables de estudio.

La temperatura (“Temperature”) registra una media de 25.83 °C y una mediana de 25.80 °C. Las temperaturas observadas en los puntos de medición oscilan en un rango relativamente estrecho de 7.5 °C, con un mínimo de 22.2 °C y un máximo de 29.7 °C.

La humedad relativa (“Relative_Humidity”) presenta una media de 71.17% y una mediana de 71.25%, con un amplio rango de 32.1 puntos porcentuales (desde un mínimo de 59.5% hasta un máximo de 91.6%), lo que evidencia una variabilidad considerable en las condiciones de humedad en el terreno de plantación.

Este comportamiento está estrechamente ligado al punto de rocío (“Dew_Point”), el cual registra un promedio de 20.16 °C y un coeficiente de variación bajo (5.62%). Los valores máximos del punto de rocío (24.0 °C) se aproximan a las temperaturas mínimas registradas en el área de estudio, sugiriendo condiciones favorables para alcanzar niveles elevados de saturación del aire y condensación.

p_hist <- aguacate_dia %>%
  dplyr::select(all_of(vars)) %>% 
  pivot_longer(everything(), names_to = "Variable", values_to = "valor") %>%
  ggplot(aes(x = valor)) +        
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "steelblue", alpha = 0.6) +
  geom_density(color = "darkred", linewidth = 0.8) +
  facet_wrap(~ Variable, scales = "free") +
  labs(title = "Ilustracion 1. Distribucion de las variables", 
       x = "Valores registrados", y = "Densidad") +
  theme_minimal()

print(p_hist)

La Ilustración 1 presenta la distribución de las variables de estudio. Se acompañó este diagnóstico visual con el cálculo del medcouple (MC) para evaluarla simetría de las distribuciones.

La temperatura (MC = 0) exhibe una distribución normal con una clara, lo que indica un comportamiento térmico estable alrededor de su media 26 °C.

La humedad relativa presentó una distribución aproximadamente simétrica (MC = -0.0278). Aunque la curva de densidad mostró aparentemente dos picos, la prueba de Hartigan no evidenció desviaciones significativas de la unimodalidad (p = 0.1668)

El punto de rocío presenta una asimetría positiva (MC = 0.143). La mayor densidad de los datos se concentra fuertemente en torno a los 20 °C, pero se extiende hacia valores más altos.

Temperatura

Semivariograma

pts <- st_as_sf(aguacate_dia, coords = c("Longitude", "Latitude"), crs = 4326)
pts <- st_transform(pts, 3115) 

coords_proj <- st_coordinates(pts)

datos_geo <- data.frame(
  X = coords_proj[,1],
  Y = coords_proj[,2],
  Temperature = pts$Temperature
)

geo_temp <- as.geodata(datos_geo, coords.col = 1:2, data.col = 3)

vario_temp <- variog(geo_temp, option = "bin", uvec = seq(0, 130, 15))

## variog: computing omnidirectional variogram

env_temp <- variog.mc.env(geo_temp, obj.variog = vario_temp, nsim = 500)

## variog.env: generating 500 simulations by permutating data values
## variog.env: computing the empirical variogram for the 500 simulations
## variog.env: computing the envelops

plot(vario_temp, envelope = env_temp, pch = 19,
     main = "Ilustración 2.Semivariograma Empírico con Envolvente Monte Carlo",
     xlab = "Distancia (m)", ylab = "Semivarianza")

La Ilustración 2 presenta el semivariograma empírico de la temperatura junto con una envolvente de simulación de Monte Carlo. Se observa que la semivarianza empirico aumenta a medida que se incrementa la distancia entre las observaciones, lo que indica que las observaciones cercanas presentan valores similares frente a las mas lejanas.

Al comparar el semivariograma empírico con la envolvente de Monte Carlo, se aprecia que varios de los valores se sitúan por afuera de las bandas del envolvente. Podemos por lo tanto rechazar la nula autocorrelación espacial y afirmar autocorrelación en los datos de temperatura.

vario_dir_zoom <- variog4(geo_temp, max.dist = 50)

## variog: computing variogram for direction = 0 degrees (0 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing omnidirectional variogram

plot(vario_dir_zoom, 
     xlab = "Distancia (m)", 
     ylab = "Semivarianza")
title(main = "Ilustración 3. Variogramas Direccionales (Hasta 50 m)")

No obstante, la Ilustración 3 presenta los semivariogramas direccionales calculados para los ángulos de 0°, 45°, 90° y 135°. Se observa que la semivarianza presenta comportamientos distintos según el angulo considerado, lo que indica que la autocorrelación espacial no es uniforme en el espacio.Estas diferencias entre las semivarianzas aputan a la construcción de un modelo anisotropico espacial frente a modelo isotrópico.

Selección y validación del modelo.

Tabla 2. Criterio de Información de Akaike (AIC) para modelos isotrópicos
Modelo	AIC
exponential	1715.13
spherical	1716.66
gaussian	NA

Tabla 3. Comparación de AIC entre modelos isotrópicos y anisotrópicos
Modelo	Tipo	AIC	Angulo_rad	Razon
spherical	Anisotrópico	1293.355	0.673	7.036
exponential	Anisotrópico	1295.190	0.657	7.104
exponential	Isotrópico	1715.129	NA	NA

Tabla 4. Resultados de Validación Cruzada (Leave-One-Out)
Modelo	RMSE	ME
Exp. Anisotrópico (Ganador predictivo)	0.6390	0.0066
Esf. Anisotrópico (Ganador ajuste)	0.6439	0.0064

Los modelos isotrópicos ajustados de la variable temperatura se estimaron mediante Máxima Verosimilitud Restringida (REML) y se compararon a través del Criterio de Información de Akaike (AIC). Los resultados muestran que el menor valor correspondió al modelo exponencial isotrópico (AIC = 1715,13), con el modelo esférico isotrópico situándose muy cerca (AIC = 1716,66). El modelo gaussiano, por su parte, no convergió a un ajuste, por lo que no se registró un valor de AIC para este.

Sin embargo, dado que la Ilustración 3 reveló la existencia de anisotropía espacial, se incorporó esta salvedad al análisis. La Tabla 2 presenta los resultados de los modelos anisotrópicos, estimados igualmente por REML, los cuales presentaron una mejora sustancial en el ajuste frente al modelo isotrópico base. El modelo esférico anisotrópico alcanzó el menor AIC (1293,36), seguido por el exponencial anisotrópico (1295,19). Ambos valores se ubican muy por debajo del valor del modelo isotrópico (1715,13), una diferencia que confirma que considerar la anisotropía produce mejores estimaciones.

Por último, la Tabla 4 recoge los resultados de la validación cruzada Leave-One-Out. El modelo exponencial anisotrópico produjo el menor Error Cuadrático Medio de predicción (RMSE = 0,6390), y en ambos casos el Error Medio (ME) resultó prácticamente igual a cero. Dado que el propósito principal del estudio es la interpolación espacial mediante Kriging, el modelo exponencial anisotrópico se perfila como la opción más adecuada gracias a su mayor precisión predictiva respaldada por el RMSE.

Resultados

paso <- 2
grilla_x <- seq(min(datos_geo$X) - 5, max(datos_geo$X) + 5, by = paso)
grilla_y <- seq(min(datos_geo$Y) - 5, max(datos_geo$Y) + 5, by = paso)
grilla_pred <- expand.grid(x = grilla_x, y = grilla_y)


krig_iso <- krige.conv(geo_temp, locations = grilla_pred, 
                       krige = krige.control(obj.model = fit_iso_ganador))

## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood

krig_aniso <- krige.conv(geo_temp, locations = grilla_pred, 
                         krige = krige.control(obj.model = fit_aniso_exp))

## krige.conv: model with constant mean
## krige.conv: anisotropy correction performed
## krige.conv: Kriging performed using global neighbourhood

paleta <- hcl.colors(50, "Inferno")
limites_temp <- range(c(krig_iso$predict, krig_aniso$predict, geo_temp$data))
cortes <- seq(limites_temp[1], limites_temp[2], length.out = 51)
colores_puntos <- paleta[findInterval(geo_temp$data, cortes)]

par(mar = c(4, 4, 3, 6)) 

image(krig_iso, col = paleta, zlim = limites_temp,
      xlab = "Este (m)", ylab = "Norte (m)")
points(geo_temp$coords, pch = 21, bg = colores_puntos, col = "white", cex = 1.2)

title(main = "Mapa 1. Distribución espacial de temperatura (Modelo Isotrópico)")

legend.krige(x.leg = c(max(grilla_x) + 2, max(grilla_x) + 5), 
             y.leg = c(min(grilla_y), max(grilla_y)), 
             values = limites_temp, vertical = TRUE, col = paleta)

par(mar = c(4, 4, 3, 6)) 

image(krig_aniso, col = paleta, zlim = limites_temp,
      xlab = "Este (m)", ylab = "Norte (m)")
points(geo_temp$coords, pch = 21, bg = colores_puntos, col = "white", cex = 1.2)

title(main = "Mapa 2. Distribución espacial de temperatura (Modelo Anisotrópico)")

legend.krige(x.leg = c(max(grilla_x) + 2, max(grilla_x) + 5), 
             y.leg = c(min(grilla_y), max(grilla_y)), 
             values = limites_temp, vertical = TRUE, col = paleta)

Finalmente el Mapa 1 (modelo isotrópico) y el Mapa 2 (modelo anisotrópico) presenta ambos modelos. Bajo el modelo con isotropía, la superficie presenta un patrón fragmentado en islas, con zonas de influencia locales en torno a cada arbol. Por su parte, el modelo anisotrópico genera bandas continuas de temperatura con una dirreción cercana a los 38°. Ambos mapas conservan el carácter de interpolador exacto del kriging, dado el efecto pepita nulo estimado en el modelo.

Humedad relativa (Relative_Humidity)

Semivariograma

coords_proj <- st_coordinates(pts)
datos_geo_rh <- data.frame(
  X = coords_proj[, 1],
  Y = coords_proj[, 2],
  Relative_Humidity = pts$Relative_Humidity
)

datos_geo_rh <- datos_geo_rh[!is.na(datos_geo_rh$Relative_Humidity), ]

geo_rh <- as.geodata(datos_geo_rh, coords.col = 1:2, data.col = 3)

vario_rh <- variog(geo_rh, option = "bin", uvec = seq(0, 130, 15))

## variog: computing omnidirectional variogram

env_rh   <- variog.mc.env(geo_rh, obj.variog = vario_rh, 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

plot(vario_rh, envelope = env_rh, pch = 19,
     main = "Ilustracion 4.Semivariograma empirico de Humedad Relativa",  
     xlab = "Distancia (m)", ylab = "Semivarianza")

vario_dir_rh <- variog4(geo_rh, max.dist = 50)

## variog: computing variogram for direction = 0 degrees (0 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing omnidirectional variogram

plot(vario_dir_rh, 
     main = "Ilustración 5. Variogramas direccionales de Humedad Relativa (Hasta 50 m)",
     xlab = "Distancia (m)", ylab = "Semivarianza")

La Ilustración 4 presenta el semivariograma empírico de la humedad relativa junto con su envolvente de Monte Carlo. Se observa que la semivarianza crece conforme se incrementa la distancia entre los puntos de observación, indicando que las observaciones cercanas presentan valores de humedad más similares. Este comportamiento evidencia la existencia de autocorrelación espacial positiva, al igual que con la temperatura se evaluo un analisis de variogramas dirrecionales, la Ilustración 5 presenta estos semivariogramas direccionales calculados para los ángulos de 0°, 45°, 90° y 135°. Se observa que la semivarianza presenta comportamientos distintos según el angulo considerado, lo que indica que la autocorrelación espacial no es uniforme en el espacio.Estas diferencias entre las semivarianzas aputan a la construcción de un modelo anisotropico espacial frente a modelo isotrópico.

Selección y validación del modelo.

var_inicial <- var(geo_rh$data, na.rm = TRUE)
ini_cov <- c(var_inicial * 0.7, 30)
ini_nugget <- var_inicial * 0.3
modelos <- c("exponential", "spherical", "gaussian")
tabla_rh_aic <- data.frame(Modelo = character(), Tipo = character(), AIC = numeric(), stringsAsFactors = FALSE)

invisible(capture.output({

  for(m in modelos) {
    num_aic <- tryCatch({
      fit <- likfit(geo_rh, ini.cov.pars = ini_cov, cov.model = m, nugget = ini_nugget, lik.method = "REML")
      fit$AIC
    }, error = function(e) return(NA))
    tabla_rh_aic <- rbind(tabla_rh_aic, data.frame(Modelo = m, Tipo = "Isotrópico", AIC = num_aic))
  }
  
  fit_iso_rh <- likfit(geo_rh, ini.cov.pars = ini_cov, cov.model = "exponential", nugget = ini_nugget, lik.method = "REML")

  for(m in modelos) {
    num_aic_a <- tryCatch({
      fit_a <- likfit(geo_rh, ini.cov.pars = ini_cov, cov.model = m, nugget = ini_nugget, 
                      lik.method = "REML", fix.psiA = FALSE, psiA = pi/4, fix.psiR = FALSE, psiR = 2)
      fit_a$AIC
    }, error = function(e) return(NA))
    tabla_rh_aic <- rbind(tabla_rh_aic, data.frame(Modelo = m, Tipo = "Anisotrópico", AIC = num_aic_a))
  }
}))

tabla_rh_final <- tabla_rh_aic[order(tabla_rh_aic$AIC), ]
rownames(tabla_rh_final) <- NULL

kable(tabla_rh_final, 
      caption = "Tabla 5. Comparación de criterios de información (AIC) para Humedad Relativa", 
      digits = 2)

Tabla 5. Comparación de criterios de información (AIC) para Humedad Relativa
Modelo	Tipo	AIC
exponential	Anisotrópico	2360.82
spherical	Anisotrópico	2361.58
exponential	Isotrópico	2829.23
spherical	Isotrópico	2840.92
gaussian	Isotrópico	NA
gaussian	Anisotrópico	NA

invisible(capture.output({
  fit_aniso_exp_rh <- likfit(geo_rh, ini.cov.pars = ini_cov, cov.model = "exponential", nugget = ini_nugget, 
                             lik.method = "REML", fix.psiA = FALSE, psiA = pi/4, fix.psiR = FALSE, psiR = 2)

  fit_aniso_sph_rh <- likfit(geo_rh, ini.cov.pars = ini_cov, cov.model = "spherical", nugget = ini_nugget, 
                             lik.method = "REML", fix.psiA = FALSE, psiA = pi/4, fix.psiR = FALSE, psiR = 2)

  cv_exp_rh <- xvalid(geo_rh, model = fit_aniso_exp_rh)
  cv_sph_rh <- xvalid(geo_rh, model = fit_aniso_sph_rh)
}))

tabla_cv_rh <- data.frame(
  Modelo = c("Exponencial Anisotrópico", "Esférico Anisotrópico"),
  RMSE = c(sqrt(mean((cv_exp_rh$data - cv_exp_rh$predicted)^2)), 
           sqrt(mean((cv_sph_rh$data - cv_sph_rh$predicted)^2))), # <-- Corregido aquí
  ME   = c(mean(cv_exp_rh$data - cv_exp_rh$predicted), 
           mean(cv_sph_rh$data - cv_sph_rh$predicted))
)

tabla_cv_rh_final <- tabla_cv_rh[order(tabla_cv_rh$RMSE), ]
rownames(tabla_cv_rh_final) <- NULL

kable(tabla_cv_rh_final, 
      caption = "Tabla 6. Métricas de validación cruzada para Humedad Relativa", 
      digits = 4)

Tabla 6. Métricas de validación cruzada para Humedad Relativa
Modelo	RMSE	ME
Esférico Anisotrópico	1.6735	0.0399
Exponencial Anisotrópico	1.6795	0.0372

Los modelos isotrópicos ajustados a la variable humedad relativa se estimaron mediante Máxima Verosimilitud Restringida (REML) y se compararon a través del Criterio de Información de Akaike (AIC). Los resultados muestran que el menor valor correspondió al modelo exponencial isotrópico (AIC = 2829,24) y el modelo esférico isotrópico situándose muy cerca (AIC = 2840,92). El modelo gaussiano, por su parte, no convergió.

Al igual que en el caso de la temperatura, los semivariogramas direccionales Ilustración 5 revelaron la existencia de anisotropía espacial,. La Tabla 5 presenta los resultados de los modelos anisotrópicos, estimados igualmente por REML.En esta variable, el menor AIC correspondió al modelo exponencial anisotrópico (AIC = 2360,82), seguido muy de cerca por el esférico anisotrópico (2361,58). Ambos valores se ubican muy por debajo del mejor modelo isotrópico.

Por último, la Tabla 6 recoge los resultados de la validación cruzada Leave-One-Out. El modelo esférico anisotrópico obtuvo el menor Error Cuadrático Medio de predicción (RMSE = 1,6736), seguido muy de cerca por el exponencial anisotrópico (RMSE = 1,6795). En ambos modelos el Error Medio (ME) resultó prácticamente igual. Dado que las diferencias entre ambos modelos son mínimas tanto en AIC y con el fin de mantener la consistencia con el modelo seleccionado para la temperatura, se optó por el modelo exponencial anisotrópico como tecnica de interpolación para la variable humeda.

Resultados

paso <- 2
grilla_x <- seq(min(datos_geo_rh$X) - 5, max(datos_geo_rh$X) + 5, by = paso)
grilla_y <- seq(min(datos_geo_rh$Y) - 5, max(datos_geo_rh$Y) + 5, by = paso)
grilla_pred <- expand.grid(x = grilla_x, y = grilla_y)

krig_aniso_rh <- krige.conv(geo_rh, locations = grilla_pred, 
                            krige = krige.control(obj.model = fit_aniso_exp_rh))

## krige.conv: model with constant mean
## krige.conv: anisotropy correction performed
## krige.conv: Kriging performed using global neighbourhood

paleta_rh <- hcl.colors(50, "Viridis")
limites_rh <- range(c(krig_aniso_rh$predict, geo_rh$data))
cortes_rh <- seq(limites_rh[1], limites_rh[2], length.out = 51)
colores_puntos_rh <- paleta_rh[findInterval(geo_rh$data, cortes_rh)]

par(mar = c(4, 4, 3, 6)) 
image(krig_aniso_rh, col = paleta_rh, zlim = limites_rh,
      main = "Mapa 3. Distribucion espacial de humedad relativa (Modelo Anisotropico)", 
      xlab = "Este (m)", ylab = "Norte (m)")
points(geo_rh$coords, pch = 21, bg = colores_puntos_rh, col = "white", cex = 1.2)
legend.krige(x.leg = c(max(grilla_x) + 2, max(grilla_x) + 5), 
             y.leg = c(min(grilla_y), max(grilla_y)), 
             values = limites_rh, vertical = TRUE, col = paleta_rh)

El Mapa 3 presenta la distribución de la humedad relativa en el area de cultivo de la plantación. Como se puede observar, la humedad presenta un patrón de líneas diagonales. Este patrón espacial es un reflejo directo de la forma en la que están sembrados los árboles en el cultivo que se puede observar en el apartado área de estudio.

En el mapa, las franjas de colores más claros indican las zonas con mayor humedad. estas franjas coinciden con las zonas deonde se encuentran los sensores (los arboles). La sombra proyectada hace que el ambiente sea más fresco, reteniendo mejor la humedad.Por el contrario, las franjas de colores oscuros corresponden a los pasillos, donde el sol golpea directamente la tierra, reduciendo la humedad relativa.

Punto de rocío (Dew_Point)

Semivariograma

coords_proj <- st_coordinates(pts)

datos_geo_dew <- data.frame(
  X = coords_proj[,1],
  Y = coords_proj[,2],
  Dew_Point = pts$Dew_Point
)

geo_dew <- as.geodata(datos_geo_dew, coords.col = 1:2, data.col = 3)

vario_dew <- variog(geo_dew, option = "bin", uvec = seq(0, 70, 10))

## variog: computing omnidirectional variogram

env_dew <- variog.mc.env(geo_dew, obj.variog = vario_dew, nsim = 500)

## variog.env: generating 500 simulations by permutating data values
## variog.env: computing the empirical variogram for the 500 simulations
## variog.env: computing the envelops

plot(vario_dew, envelope = env_dew, pch = 19,
     main = "Ilustracion 6. Semivariograma empirico de Punto de Rocio (Omnidireccional)",
     xlab = "Distancia (m)", ylab = "Semivarianza")

vario_dir_dew <- variog4(geo_dew, max.dist = 50)

## variog: computing variogram for direction = 0 degrees (0 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing omnidirectional variogram

plot(vario_dir_dew, 
     main = "Ilustración 7. Variogramas direccionales de Punto de Rocío (Hasta 50 m)",
     xlab = "Distancia (m)", ylab = "Semivarianza")

La Ilustración 6 presenta el semivariograma omnidireccional del punto de rocío, en el que los valores se ubican mayoritariamente dentro de la envolvente de Monte Carlo, lo que en una primera lectura sugeriría ausencia de autocorrelación espacial. No obstante, el semivariograma direccional (Ilustración 7) muestra que en la dirección de 45° (NE-SO) la semivarianza se mantiene baja a cortas distancias, evidenciando una estructura de continuidad espacial que no aparece en las direcciones perpendiculares. Por tanto, el comportamiento plano del semivariograma omnidireccional se interpreta como el resultado de promediar una dirección estructurada (45°) con direcciones sin estructura, lo que diluye la señal. El punto de rocío presenta así una dependencia espacial marcadamente anisotrópica, coincidente con la dirección detectada en la temperatura y la humedad relativa.

Selección y validación del modelo.

Tabla 7. Comparación de criterios de información (AIC) para el Punto de Rocío
Modelo	Tipo	AIC	Angulo_rad	Razon
gaussian	Anisotrópico	1227.301	0.643	4.934
exponential	Anisotrópico	1245.198	0.645	5.991
spherical	Anisotrópico	1268.195	0.633	5.659
gaussian	Isotrópico	1542.057	NA	NA
exponential	Isotrópico	1544.522	NA	NA
spherical	Isotrópico	1554.193	NA	NA

invisible(capture.output({
  fit_aniso_exp_dew <- likfit(geo_dew, ini.cov.pars = ini_cov_dew, cov.model = "exponential", 
                              nugget = ini_nugget_dew, lik.method = "REML", 
                              fix.psiA = FALSE, psiA = pi/4, fix.psiR = FALSE, psiR = 2)

  fit_aniso_sph_dew <- likfit(geo_dew, ini.cov.pars = ini_cov_dew, cov.model = "spherical", 
                              nugget = ini_nugget_dew, lik.method = "REML", 
                              fix.psiA = FALSE, psiA = pi/4, fix.psiR = FALSE, psiR = 2)

  fit_aniso_gau_dew <- likfit(geo_dew, ini.cov.pars = ini_cov_dew, cov.model = "gaussian", 
                              nugget = ini_nugget_dew, lik.method = "REML", 
                              fix.psiA = FALSE, psiA = pi/4, fix.psiR = FALSE, psiR = 2)

  cv_exp_dew <- xvalid(geo_dew, model = fit_aniso_exp_dew)
  cv_sph_dew <- xvalid(geo_dew, model = fit_aniso_sph_dew)
  cv_gau_dew <- xvalid(geo_dew, model = fit_aniso_gau_dew) # Validación cruzada Gaussiano
}))

tabla_cv_dew <- data.frame(
  Modelo = c("Exponencial Anisotrópico", "Esférico Anisotrópico", "Gaussiano Anisotrópico"),
  RMSE = c(sqrt(mean((cv_exp_dew$data - cv_exp_dew$predicted)^2)), 
           sqrt(mean((cv_sph_dew$data - cv_sph_dew$predicted)^2)),
           sqrt(mean((cv_gau_dew$data - cv_gau_dew$predicted)^2))),
  ME   = c(mean(cv_exp_dew$data - cv_exp_dew$predicted), 
           mean(cv_sph_dew$data - cv_sph_dew$predicted),
           mean(cv_gau_dew$data - cv_gau_dew$predicted))
)

tabla_cv_dew_final <- tabla_cv_dew[order(tabla_cv_dew$RMSE), ]
rownames(tabla_cv_dew_final) <- NULL

kable(tabla_cv_dew_final, 
      caption = "Tabla 8. Métricas de validación cruzada para el Punto de Rocío", 
      digits = 4)

Tabla 8. Métricas de validación cruzada para el Punto de Rocío
Modelo	RMSE	ME
Gaussiano Anisotrópico	0.6087	0.0141
Esférico Anisotrópico	0.6191	0.0116
Exponencial Anisotrópico	0.6240	0.0129

Los modelos isotrópicos ajustados a la variable punto de rocío se estimaron mediante Máxima Verosimilitud Restringida (REML) y se compararon a través del Criterio de Información de Akaike (AIC). Los resultados muestran que el menor valor correspondió al modelo gaussiano isotrópico (AIC = 1542,06) y el modelo exponencial isotrópico situándose muy cerca (AIC = 1544,52). El modelo esférico, por su parte, registró un valor ligeramente superior (AIC = 1554,19), logrando todos la convergencia.

Al igual que en el caso de la temperatura y la humedad, los semivariogramas direccionales (Ilustración X) revelaron la existencia de anisotropía espacial. La Tabla 7 presenta los resultados de los modelos anisotrópicos, estimados igualmente por REML. En esta variable, el menor AIC correspondió al modelo gaussiano anisotrópico (AIC = 1227,30), seguido por el exponencial anisotrópico (AIC = 1245,20).

Por último, la Tabla 8 recoge los resultados de la validación cruzada Leave-One-Out. El modelo gaussiano anisotrópico obtuvo el menor (RMSE = 0,6087), seguido por el esférico anisotrópico (RMSE = 0,6191) y el exponencial anisotrópico (RMSE = 0,6240). En los tres modelos el Error Medio (ME) resultó prácticamente igual. Dado que el modelo gaussiano anisotrópico presentó los mejores resultados, se optó por este modelo como técnica de interpolación para la variable punto de rocío.

Resultados

paso <- 2
grilla_x <- seq(min(datos_geo_dew$X) - 5, max(datos_geo_dew$X) + 5, by = paso)
grilla_y <- seq(min(datos_geo_dew$Y) - 5, max(datos_geo_dew$Y) + 5, by = paso)
grilla_pred <- expand.grid(x = grilla_x, y = grilla_y)

krig_aniso_dew <- krige.conv(geo_dew, locations = grilla_pred, 
                             krige = krige.control(obj.model = fit_aniso_gau_dew))

## krige.conv: model with constant mean
## krige.conv: anisotropy correction performed
## krige.conv: Kriging performed using global neighbourhood

paleta_dew <- hcl.colors(50, "YlGnBu", rev = TRUE) 
limites_dew <- range(c(krig_aniso_dew$predict, geo_dew$data))
cortes_dew <- seq(limites_dew[1], limites_dew[2], length.out = 51)
colores_puntos_dew <- paleta_dew[findInterval(geo_dew$data, cortes_dew)]

par(mar = c(4, 4, 3, 6)) 

image(krig_aniso_dew, col = paleta_dew, zlim = limites_dew,
      main = "Mapa 5. Distribución espacial del Punto de Rocío (Modelo Anisotrópico)", 
      xlab = "Este (m)", ylab = "Norte (m)")

points(geo_dew$coords, pch = 21, bg = colores_puntos_dew, col = "white", cex = 1.2)

legend.krige(x.leg = c(max(grilla_x) + 2, max(grilla_x) + 5), 
             y.leg = c(min(grilla_y), max(grilla_y)), 
             values = limites_dew, vertical = TRUE, col = paleta_dew)

El Mapa 5 presenta la distribución espacial del punto de rocío. Tal como se evidenció en las variables anteriores, el mapa dibuja un patrón indiscutible de franjas . Lo que evidenciaria que la organización cultivo es el factor que controla las variables climaticas del lote de aguacates.

En el mapa las franjas de color azul oscuro representan las zonas con el punto de rocío más alto las cuales corresponden a las hileras de los árboles de aguacate. Por el contrario, las franjas de tonos claros marcan los pasillos.

Conclusiones

Los Mapas 1 a 5 presentan la interpolación espacial de las variables seleccionadas mediante la metodología de kriging. A lo largo del análisis se evidenció la importancia de evaluar los semivariogramas omnidireccionales como direccionales para caracterizar el comportamiento de la autocorrelación espacial de cada variable, pues fue esta evaluación la que permitió detectar y modelar la anisotropia.

Cabe aclarar que las mediciones corresponden a una fecha y hora específicas (1 de octubre de 2020, aproximadamente las 10:11 a. m.), instante en el que el sol se encontraba en un azimut de 102° (este-sureste) y una altura de 63°. Es muy probable que esta dispocisión , sumada a la organización del cultivo en hileras más la inclinación del terreno, expliquen el patrón anisotrópico. No obstante, dado que se desconoce el objetivo específico del proyecto del que provienen los datos, esta interpretación se plantea como una hipótesis y no como una conclusión definitiva.

En consecuencia, la interpretación de los resultados resulta limitada por el desconocimiento de la materia de estudio. Un análisis más completo exigiría un abordaje interdisciplinario que integre conocimientos agronómicos, climatológicos. Para trabajos futuros, sería importante incorporar la dimensión temporal mediante análisis de series de tiempo (En distintas horas , días y estaciones) así como modelos que incluyan covariables adicionales capaces de explicar la variabilidad espacial de forma más realista.

Actividad 3 - Interpolación mediante metodología Kriging

Introducción

Área de estudio

Descripción de los datos

Temperatura

Semivariograma

Selección y validación del modelo.

Resultados

Humedad relativa (Relative_Humidity)

Semivariograma

Selección y validación del modelo.

Resultados

Punto de rocío (Dew_Point)

Semivariograma

Selección y validación del modelo.

Resultados

Conclusiones