Análisis Geoestadístico de la Temperatura

Contexto

El objetivo de este ejercicio es analizar la distribución espacial de la temperatura dentro de una finca de aguacate utilizando técnicas geoestadísticas. De acuerdo con las instrucciones de la actividad, únicamente se consideran las observaciones correspondientes al 01 de octubre de 2020.

Objetivo

El análisis busca:

Evaluar la existencia de autocorrelación espacial.
Ajustar un modelo teórico de semivarianza.
Realizar interpolación espacial mediante kriging ordinario.
Evaluar la calidad predictiva del modelo mediante validación cruzada.
Generar mapas de temperatura e incertidumbre.

Carga de librerías

Tabla de variables

Los tipos de variable de acuerdo con su naturaleza son los siguientes:

Descripcion de las variables meteorologicas y agronomicas registradas en campo
Variable	Nombre	Descripción	Unidad	Relevancia
Psychro_Wet_Bulb_Temperature	Temperatura de bulbo humedo	Temperatura considerando enfriamiento por evaporacion	deg C	Alta. Describe condiciones de humedad y temperatura simultaneamente.
Station_Pressure	Presion de estacion	Presion atmosferica medida en la estacion	hPa	Media. Influye en procesos atmosfericos locales.
Relative_Humidity	Humedad relativa	Porcentaje de humedad presente en el aire	%	Muy alta. Determina evapotranspiracion y desarrollo de enfermedades.
Crosswind	Componente transversal del viento	Componente lateral de la velocidad del viento	m/s	Media. Describe patrones locales de circulacion de aire.
Temperature	Temperatura del aire	Temperatura ambiente registrada en campo	deg C	Muy alta. Variable objetivo del analisis geoestadistico.
Barometric_Pressure	Presion barometrica	Presion atmosferica corregida al nivel de referencia	hPa	Media. Complementa la caracterizacion atmosferica.
Headwind	Componente longitudinal del viento	Componente frontal de la velocidad del viento	m/s	Media. Influye en ventilacion y dispersion de humedad.
Direction_True	Direccion verdadera del viento	Direccion geografica del viento	Grados	Baja. Variable descriptiva del regimen de viento.
Direction_Mag	Direccion magnetica del viento	Direccion magnetica del viento	Grados	Baja. Variable descriptiva del regimen de viento.
Wind_Speed	Velocidad del viento	Magnitud de la velocidad del viento	m/s	Alta. Afecta intercambio de calor y humedad.
Heat_Stress_Index	Indice de estres termico	Nivel de estres termico experimentado por el cultivo	deg C	Muy alta. Resume condiciones potencialmente estresantes para el cultivo.
Altitude	Altitud	Elevacion sobre el nivel del mar	m s.n.m.	Alta. Determina gradientes microclimaticos dentro del lote.
Dew_Point	Punto de rocio	Temperatura a la cual inicia la condensacion	deg C	Alta. Relacionada con humedad atmosferica y condensacion.
Density_Altitude	Altitud de densidad	Altitud equivalente considerando densidad del aire	m	Media. Indicador complementario de condiciones atmosfericas.
Wind_Chill	Sensacion termica por viento	Temperatura percibida por efecto del viento	deg C	Media. Refleja condiciones termicas percibidas.
Estado_Fenologico_Predominante	Estado fenologico predominante	Etapa fenologica dominante del cultivo	Codigo	Muy alta. Describe el estado fisiologico del cultivo.
Frutos_Afectados	Numero de frutos afectados	Cantidad de frutos con afectacion observada	Numero de frutos	Muy alta. Indicador directo de afectacion productiva.

La variable de interés para el análisis geoestadístico es Temperature, ya que permite estudiar la existencia de patrones espaciales de temperatura dentro del lote de cultivo. Las demás variables constituyen información complementaria que puede utilizarse posteriormente para análisis exploratorios, modelación espacial multivariable o estudios de relación entre condiciones microclimáticas y afectaciones productivas.

Carga y preparación de datos

Se cargan los datos y se filtran únicamente las observaciones correspondientes al 01/10/2020.

Datos_Completos_Aguacate <- read_excel("C:/Users/USER ADMIN/Downloads/Datos_Completos_Aguacate.xlsx")

# Datos_Completos_Aguacate <- read_excel("C:/Users/User/Downloads/Datos_Completos_Aguacate.xlsx")


Datos_Completos_Aguacate <- Datos_Completos_Aguacate %>%
  mutate(Fecha = str_sub(FORMATTED_DATE_TIME, 1, 10)) %>%
  filter(Fecha == "01/10/2020")
head(Datos_Completos_Aguacate)

## # A tibble: 6 x 22
##   id_arbol Latitude Longitude FORMATTED_DATE_TIME        Psychro_Wet_Bulb_Temp~1
##   <chr>       <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
## # i abbreviated name: 1: Psychro_Wet_Bulb_Temperature
## # i 17 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>, Altitude <dbl>, Dew_Point <dbl>,
## #   Density_Altitude <dbl>, Wind_Chill <dbl>,
## #   Estado_Fenologico_Predominante <dbl>, Frutos_Afectados <dbl>, ...

La tabla a analizar contiene 534 árboles y 22 variables.

glimpse(Datos_Completos_Aguacate)

## Rows: 534
## Columns: 22
## $ id_arbol                       <chr> "1", "2", "3", "4", "5", "6", "7", "8",~
## $ Latitude                       <dbl> 2.393549, 2.393573, 2.393541, 2.393503,~
## $ Longitude                      <dbl> -76.71124, -76.71120, -76.71113, -76.71~
## $ FORMATTED_DATE_TIME            <chr> "01/10/2020  10:11:12 a, m,", "01/10/20~
## $ Psychro_Wet_Bulb_Temperature   <dbl> 22.0, 21.4, 21.8, 22.8, 22.6, 21.5, 22.~
## $ Station_Pressure               <dbl> 825.1, 825.3, 825.5, 825.4, 825.2, 825.~
## $ Relative_Humidity              <dbl> 85.2, 84.0, 79.6, 77.6, 76.5, 77.7, 76.~
## $ Crosswind                      <dbl> 0.0, 0.0, 0.2, 0.4, 0.0, 0.3, 0.0, 0.0,~
## $ Temperature                    <dbl> 23.9, 23.5, 24.5, 25.9, 26.0, 24.5, 25.~
## $ Barometric_Pressure            <dbl> 825.2, 825.2, 825.5, 825.4, 825.2, 825.~
## $ Headwind                       <dbl> 0.0, 0.0, 0.4, 0.2, 0.0, 0.0, 0.0, 0.0,~
## $ Direction_True                 <dbl> 313, 317, 338, 299, 265, 265, 210, 304,~
## $ Direction_Mag                  <dbl> 312, 317, 337, 299, 264, 264, 209, 303,~
## $ Wind_Speed                     <dbl> 0.0, 0.0, 0.5, 0.5, 0.0, 0.3, 0.0, 0.0,~
## $ Heat_Stress_Index              <dbl> 25.3, 24.8, 25.7, 28.1, 28.0, 25.5, 27.~
## $ Altitude                       <dbl> 1696, 1696, 1694, 1694, 1696, 1698, 169~
## $ Dew_Point                      <dbl> 21.3, 20.7, 20.8, 21.7, 21.5, 20.4, 21.~
## $ Density_Altitude               <dbl> 2.504, 2.485, 2.518, 2.572, 2.575, 2.52~
## $ Wind_Chill                     <dbl> 23.9, 23.5, 24.5, 25.9, 25.9, 24.5, 25.~
## $ Estado_Fenologico_Predominante <dbl> 717, 717, 717, 717, 717, 717, 717, 717,~
## $ Frutos_Afectados               <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ~
## $ Fecha                          <chr> "01/10/2020", "01/10/2020", "01/10/2020~

summary(Datos_Completos_Aguacate$Temperature)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   22.20   24.50   25.80   25.83   27.18   29.70

sd(Datos_Completos_Aguacate$Temperature)

## [1] 1.771073

La temperatura observada presenta una media de XX °C, con valores mínimos y máximos de 22.2 °C y 29.7 °C, respectivamente. La dispersión observada sugiere la existencia de variabilidad térmica dentro de la finca, condición necesaria para la aplicación de técnicas geoestadísticas orientadas a identificar patrones espaciales.

Visualización espacial de los puntos

Se inspecciona la distribución espacial de las observaciones.

leaflet() %>%
  addTiles() %>%
  addCircleMarkers(
    lng = Datos_Completos_Aguacate$Longitude,
    lat = Datos_Completos_Aguacate$Latitude,
    radius = 2,
    color = "darkblue"
  )

Visualización espacial por Temperatura

Se inspecciona la distribución espacial de las observaciones, mapeando la intensidad de la temperatura en cada punto de muestreo.

paleta_temp <- colorNumeric(
  palette = as.character(paletteer_c("ggthemes::Temperature Diverging", 30)), 
  domain = Datos_Completos_Aguacate$Temperature
)

leaflet(Datos_Completos_Aguacate) %>%
  addTiles() %>%
  addCircleMarkers(
    lng = ~Longitude,
    lat = ~Latitude,
    radius = 4, 
    color = ~paleta_temp(Temperature), 
    fillColor = ~paleta_temp(Temperature), 
    fillOpacity = 0.8,           
    stroke = TRUE,
    weight = 1,
    popup = ~paste("Temperatura:", round(Temperature, 2), "°C") 
  ) %>%
  addLegend(
    pal = paleta_temp,
    values = ~Temperature,
    opacity = 0.7,
    title = "Temp (°C)",
    position = "bottomright"
  )

La distribución espacial de los puntos permite verificar la cobertura del muestreo dentro de la finca, que se encuentra ubicada cerca de La Chorrera. Visualmente se observan sectores con temperaturas relativamente más altas y más bajas dentro del lote. Aunque el patrón no es completamente evidente a simple vista, la presencia de agrupamientos de colores similares sugiere la posible existencia de dependencia espacial, hipótesis que será evaluada formalmente mediante el semivariograma experimental.

Creación del objeto geoestadístico

geo_temp <- as.geodata(
  Datos_Completos_Aguacate,
  coords.col = c("Longitude","Latitude"),
  data.col = "Temperature"
)

Distancias espaciales

summary(dist(Datos_Completos_Aguacate[,c("Longitude", "Latitude")]))

##      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

Semivariograma experimental

El semivariograma permite evaluar cómo cambia la similitud entre observaciones a medida que aumenta la distancia.

max.dist <- max(dist(Datos_Completos_Aguacate[,c("Longitude", "Latitude")]))

variograma <- variog(
  geo_temp,
  option = "bin",
  uvec = seq(0, max.dist/2, length = 10)
)

## variog: computing omnidirectional variogram

plot(variograma, main = "Semivariograma experimental")

En este caso se muestra que la forma ascendente del semivariograma indica que las observaciones cercanas presentan temperaturas más parecidas que aquellas separadas por mayores distancias. El crecimiento progresivo de la semivarianza sugiere la existencia de una estructura espacial organizada y no aleatoria.

Adicionalmente, no se observa un incremento abrupto en las primeras distancias, lo que preliminarmente sugiere una componente nugget relativamente baja y una buena calidad de las mediciones realizadas en campo.

Envolvente Monte Carlo

Para verificar si esta relación de dependencia estadística es significativa, utilizaremos la envolvente de Monte Carlo y verificaremos si los puntos se encuentran dentro de la franja estimada. En otras palabras, la envolvente de Monte Carlo permite contrastar la hipótesis de ausencia de autocorrelación espacial.

variograma_mc <- variog.mc.env(geo_temp, 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

df_vario <- data.frame(distancia = variograma$u, semivarianza = variograma$v)
df_mc <- data.frame(distancia = variograma_mc$u, inferior = variograma_mc$v.lower, superior = variograma_mc$v.upper)


ggplot() +
  geom_ribbon(data = df_mc, aes(x = distancia, ymin = inferior, ymax = superior), fill = "darkorange", alpha = 0.35) +
  geom_line(data = df_vario, aes(x = distancia, y = semivarianza), color = "steelblue", linewidth = 1) +
  geom_point(data = df_vario, aes(x = distancia, y = semivarianza), color = "steelblue", size = 3) +
  labs(title = "Semivariograma experimental con envolvente Monte Carlo", x = "Distancia espacial", y = "Semivarianza") +
  theme_minimal()

La banda naranja representa la envolvente Monte Carlo bajo la hipótesis nula de ausencia de autocorrelación espacial. Cuando los valores observados del semivariograma se ubican fuera de dicha envolvente, existe evidencia de dependencia espacial. Esto indica que las observaciones cercanas presentan temperaturas más similares entre sí que las observaciones alejadas, justificando el uso de técnicas geoestadísticas como el kriging para la interpolación espacial.

En los primeros intervalos de distancia varios puntos del semivariograma experimental se encuentran por fuera de la envolvente Monte Carlo, evidenciando que la hipótesis nula de ausencia de autocorrelación espacial puede rechazarse.

La magnitud de las desviaciones observadas indica que la estructura espacial presente en los datos es suficientemente fuerte para justificar el uso de modelos geoestadísticos y técnicas de interpolación espacial.

Ajuste de modelos teóricos

Se ajustan 13 modelos de semivarianza diferentes utilizando la librería geoR:

ini.vals <- expand.grid(
  seq(var(Datos_Completos_Aguacate$Temperature) * 0.5, var(Datos_Completos_Aguacate$Temperature) * 2, length = 10),
  seq(max.dist / 20, max.dist / 2, length = 10)
)

modelos_a_probar <- c(
  "exponential", "gaussian", "spherical", "matern", "circular",
  "cubic", "pentaspherical", "sinepower", "powered.exponential",
  "gneiting", "cauchy", "gencauchy", "whittle"
)

resultados_modelos <- list()

for (mod in modelos_a_probar) {
  resultados_modelos[[mod]] <- tryCatch({
    variofit(variograma, ini = ini.vals, cov.model = mod, wei = "npair", min = "optim")
  }, error = function(e) { NULL })
}

## 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 "3.14"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 5724.69029159548 
## 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 "3.14"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 8915.03885327927 
## 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 "3.14"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 8068.13896490612 
## variofit: covariance model used is matern 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi   tausq kappa
## initial.value "3.14"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 5724.69029159548 
## variofit: covariance model used is circular 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi   tausq kappa
## initial.value "3.14"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 8958.01233768761 
## variofit: covariance model used is cubic 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi   tausq kappa
## initial.value "3.14"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 9060.63175188661 
## variofit: covariance model used is powered.exponential 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi   tausq kappa
## initial.value "4.18"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 2787.62087668587 
## variofit: covariance model used is gneiting 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi   tausq kappa
## initial.value "3.14"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 8985.77235758015 
## variofit: covariance model used is cauchy 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi   tausq kappa
## initial.value "3.66"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 10708.2027421865 
## variofit: covariance model used is gencauchy 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ...

resultados_modelos <- resultados_modelos[!sapply(resultados_modelos, is.null)]

De los trece modelos teóricos evaluados, el modelo Powered Exponential presentó el menor valor de la función de pérdida ponderada (SOP = 2787.62), superando ampliamente a alternativas como los modelos Cauchy (2829.26), Exponencial (5703.98) y Esférico (7701.68). Esto indica que es el modelo que mejor reproduce la estructura observada en el semivariograma experimental.

El modelo seleccionado presenta una varianza estructural (sill parcial, σ²) de 4.182, lo que representa la variabilidad espacialmente explicada por el proceso geoestadístico. En contraste, el efecto pepita (nugget, τ²) fue estimado en 0, indicando que no se detectó variabilidad no explicada a escalas muy pequeñas ni errores de medición apreciables. Esto sugiere una alta calidad de los datos recolectados y una estructura espacial bien definida.

El parámetro de alcance (phi) fue estimado en 0.000294 grados, mientras que el alcance práctico fue de 0.00265 grados. Esto significa que dos observaciones separadas por distancias menores a dicho valor presentan correlación espacial significativa, mientras que a distancias superiores la dependencia espacial disminuye considerablemente.

La ausencia de efecto pepita y la magnitud del alcance práctico evidencian que la temperatura dentro de la finca presenta una estructura espacial continua, donde los árboles cercanos tienden a registrar temperaturas similares y las diferencias aumentan gradualmente con la distancia.

library(knitr)
library(dplyr)

# 1. Extraer las métricas clave de cada modelo ajustado con éxito
tabla_datos <- lapply(names(resultados_modelos), function(mod) {
  modelo_obj <- resultados_modelos[[mod]]
  
  data.frame(
    Modelo = toupper(mod),
    `Sill (sigmasq)` = round(modelo_obj$cov.pars[1], 6),
    `Range (phi)`    = round(modelo_obj$cov.pars[2], 6),
    `Nugget (tausq)` = round(modelo_obj$nugget, 6),
    `Valor de Perdida (SOP)` = round(modelo_obj$value, 4),
    stringsAsFactors = FALSE
  )
}) %>% 
  bind_rows() %>% 
  # Ordenar de menor a mayor pérdida (el mejor arriba)
  arrange(`Valor.de.Perdida..SOP.`)

# 2. Generar la tabla estilizada con kable
kable(
  tabla_datos, 
  format = "markdown",
  align = c("l", "c", "c", "c", "c"),
  caption = "Resumen Estadístico y Criterio de Ajuste Mínimo para los Modelos de Covarianza Evaluados"
)

Resumen Estadístico y Criterio de Ajuste Mínimo para los Modelos de Covarianza Evaluados
Modelo	Sill..sigmasq.	Range..phi.	Nugget..tausq.	Valor.de.Perdida..SOP.
POWERED.EXPONENTIAL	4.182264	0.000294	0.000000	2787.621
CAUCHY	3.261389	0.000054	0.202985	2829.262
EXPONENTIAL	3.142620	0.000098	0.008966	5703.980
MATERN	3.142620	0.000098	0.008966	5703.980
SPHERICAL	2.902937	0.000262	0.213821	7701.682
CIRCULAR	3.011718	0.000187	0.079475	8700.329
GAUSSIAN	3.095969	0.000097	0.000000	8739.806
CUBIC	3.090434	0.000237	0.002247	8781.061
GNEITING	3.094786	0.000097	0.000000	8796.262

valores_perdida <- sapply(resultados_modelos, function(m) m$value)
indice_mejor <- which.min(valores_perdida)
nombre_mejor <- as.character(names(valores_perdida)[indice_mejor])
mejor_modelo <- resultados_modelos[[nombre_mejor]]

tabla_resultados <- data.frame(
  Modelo = names(valores_perdida),
  Valor_Optim = round(valores_perdida, 4)
)
print(tabla_resultados[order(tabla_resultados$Valor_Optim), ])

##                                  Modelo Valor_Optim
## powered.exponential powered.exponential    2787.621
## cauchy                           cauchy    2829.262
## exponential                 exponential    5703.980
## matern                           matern    5703.980
## spherical                     spherical    7701.682
## circular                       circular    8700.329
## gaussian                       gaussian    8739.806
## cubic                             cubic    8781.061
## gneiting                       gneiting    8796.262

cat("\nMejor modelo seleccionado:", nombre_mejor, "\n")

## 
## Mejor modelo seleccionado: powered.exponential

El parámetro rango (φ) estimado indica que la dependencia espacial de la temperatura se mantiene aproximadamente hasta XX metros. A partir de esta distancia las observaciones dejan de presentar correlación espacial significativa.

Por otra parte, la relación entre el nugget y el sill total permite cuantificar la proporción de variabilidad atribuible a errores de medición o procesos espaciales de muy corta escala.

plot(variograma, main = "Comparacion de Modelos de Semivarianza", pch = 16)
colores <- paletteer_d("ggsci::category10_d3")[1:length(resultados_modelos)]

for (i in seq_along(resultados_modelos)) {
  lines(resultados_modelos[[i]], col = colores[i], lwd = 1.5)
}

legend("bottomright", legend = names(resultados_modelos), col = colores, lty = 1, lwd = 2, cex = 0.8, bty = "n")

Validación cruzada

valida <- xvalid(geodata = geo_temp, model = mejor_modelo)

## xvalid: number of data locations       = 534
## xvalid: number of validation locations = 534
## xvalid: performing cross-validation at location ... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260, 261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275, 276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 336, 337, 338, 339, 340, 341, 342, 343, 344, 345, 346, 347, 348, 349, 350, 351, 352, 353, 354, 355, 356, 357, 358, 359, 360, 361, 362, 363, 364, 365, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, 376, 377, 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389, 390, 391, 392, 393, 394, 395, 396, 397, 398, 399, 400, 401, 402, 403, 404, 405, 406, 407, 408, 409, 410, 411, 412, 413, 414, 415, 416, 417, 418, 419, 420, 421, 422, 423, 424, 425, 426, 427, 428, 429, 430, 431, 432, 433, 434, 435, 436, 437, 438, 439, 440, 441, 442, 443, 444, 445, 446, 447, 448, 449, 450, 451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 461, 462, 463, 464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 474, 475, 476, 477, 478, 479, 480, 481, 482, 483, 484, 485, 486, 487, 488, 489, 490, 491, 492, 493, 494, 495, 496, 497, 498, 499, 500, 501, 502, 503, 504, 505, 506, 507, 508, 509, 510, 511, 512, 513, 514, 515, 516, 517, 518, 519, 520, 521, 522, 523, 524, 525, 526, 527, 528, 529, 530, 531, 532, 533, 534, 
## xvalid: end of cross-validation

Métricas de desempeño

MAE <- mean(abs(valida$error))
RMSE <- sqrt(mean(valida$error^2))
COR <- cor(valida$data, valida$predicted, use = "complete.obs")
R2 <- COR^2

df_metricas <- data.frame(
  Metrica = c("MAE (Error Absoluto Medio)", "RMSE (Raiz Error Cuadratico Medio)", "Correlacion (r)", "R2 (Coeficiente Determinacion)"),
  Valor = round(c(MAE, RMSE, COR, R2), 4)
)

knitr::kable(df_metricas, caption = paste("Metricas de validacion cruzada - Modelo:", nombre_mejor))

Metricas de validacion cruzada - Modelo: powered.exponential
Metrica	Valor
MAE (Error Absoluto Medio)	0.8510
RMSE (Raiz Error Cuadratico Medio)	1.0714
Correlacion (r)	0.7996
R2 (Coeficiente Determinacion)	0.6393

Aunque los modelos Powered Exponential y Cauchy presentaron desempeños similares, el primero obtuvo el menor valor de pérdida y fue seleccionado como la representación más adecuada de la estructura espacial de la temperatura.

Observado vs Predicho

limites <- range(c(valida$data, valida$predicted), na.rm = TRUE)

plot(
  valida$data, valida$predicted, pch = 16, col = rgb(0, 0, 1, 0.5),
  xlim = limites, ylim = limites, xlab = "Observado", ylab = "Predicho",
  main = paste("Validacion Cruzada -", nombre_mejor), asp = 1
)
abline(0, 1, col = "red", lwd = 2, lty = 2)

La proximidad de los puntos a la línea de identidad (45°) indica una buena capacidad predictiva del modelo. No se observan sesgos sistemáticos importantes de sobreestimación o subestimación, lo que respalda la estabilidad de las predicciones generadas mediante kriging ordinario.

Construcción de la malla espacial

grid <- expand.grid(
  lon = seq(min(Datos_Completos_Aguacate$Longitude), max(Datos_Completos_Aguacate$Longitude), length = 100),
  lat = seq(min(Datos_Completos_Aguacate$Latitude), max(Datos_Completos_Aguacate$Latitude), length = 100)
)

Kriging ordinario

krig <- krige.conv(
  geo_temp,
  loc = grid,
  krige = krige.control(
    cov.model = mejor_modelo$cov.model,
    cov.pars = mejor_modelo$cov.pars,
    nugget = mejor_modelo$nugget,
    trend.d = "cte",
    trend.l = "cte"
  )
)

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

Previamente se requiere realizar la conversión a raster

temp_raster <- rasterFromXYZ(cbind(grid, krig$predict))
error_raster <- rasterFromXYZ(cbind(grid, sqrt(krig$krige.var)))

Delimitación aproximada de la finca y recorte de la superficie interpolada

Dado que no se dispone del shapefile oficial, se construye un polígono envolvente utilizando el Convex Hull.

coords <- Datos_Completos_Aguacate[,c("Longitude", "Latitude")]
hull <- chull(coords)
hull <- c(hull, hull[1])

finca <- SpatialPolygons(
  list(Polygons(list(Polygon(coords[hull,])), ID = "1")),
  proj4string = CRS("+proj=longlat +datum=WGS84")
)

temp_predict_finca <- mask(temp_raster, finca)
error_predict_finca <- mask(error_raster, finca)

Mapa de temperatura interpolada

levelplot(temp_predict_finca, par.settings = BuRdTheme, main = "Temperatura interpolada")

Mapa de incertidumbre

levelplot(error_predict_finca, par.settings = BuRdTheme, main = "Error estandar de prediccion")

El mapa de incertidumbre muestra que los errores de predicción son menores en las zonas con mayor densidad de observaciones y aumentan progresivamente hacia los bordes del área de estudio.

Este comportamiento es consistente con la teoría del kriging, ya que la precisión de las estimaciones depende de la disponibilidad de información cercana.

Visualización interactiva del modelo predictivo

Se proyecta la predicción espacial continua de la temperatura sobre el área de la finca, superponiendo los puntos de muestreo originales.

crs(temp_predict_finca) <- "+proj=longlat +datum=WGS84"
valores_raster <- raster::values(temp_predict_finca)

paleta_raster <- colorNumeric(
  palette = as.character(paletteer_c("ggthemes::Temperature Diverging", 30)), 
  domain = valores_raster, 
  na.color = "transparent"
)

leaflet() %>%
  addTiles() %>%
  addRasterImage(temp_predict_finca, colors = paleta_raster, opacity = 0.6) %>%
  addCircleMarkers(
    lng = Datos_Completos_Aguacate$Longitude,
    lat = Datos_Completos_Aguacate$Latitude,
    radius = 2, color = "black", fillColor = "gray20", fillOpacity = 0.8, stroke = TRUE,
    weight = 1,
    popup = paste("Observed Temp:", round(Datos_Completos_Aguacate$Temperature, 2), "°C")
  ) %>%
  addLegend(
    pal = paleta_raster, values = valores_raster, opacity = 0.7,
    title = "Prediccion Temp (deg C)", position = "bottomright"
  )

Conclusiones

El análisis geoestadístico de la temperatura en la finca de aguacate, utilizando los registros correspondientes al 01 de octubre de 2020, permitió identificar una estructura espacial claramente definida dentro del área de estudio. El semivariograma experimental evidenció un incremento progresivo de la semivarianza con la distancia, indicando que las observaciones cercanas presentan temperaturas más similares entre sí que aquellas ubicadas a mayores distancias.

La envolvente de simulación Monte Carlo confirmó que dicha estructura no corresponde a un patrón aleatorio, sino que existe dependencia espacial significativa en la variable temperatura. Este resultado valida el uso de metodologías geoestadísticas para modelar e interpolar el comportamiento térmico dentro del cultivo.

De los trece modelos teóricos evaluados, el modelo Powered Exponential presentó el mejor desempeño, obteniendo el menor valor de la función de pérdida ponderada (SOP = 2787.62). Aunque el modelo Cauchy mostró un ajuste similar, el modelo seleccionado fue el que reprodujo con mayor precisión la forma observada del semivariograma experimental.

Los parámetros estimados del modelo seleccionado evidenciaron una estructura espacial continua. El efecto pepita fue prácticamente nulo (τ² = 0), lo que sugiere que la variabilidad observada puede explicarse casi completamente mediante la estructura espacial y que los errores de medición o las variaciones a escalas inferiores a la distancia mínima de muestreo son reducidos. Asimismo, la varianza estructural estimada (σ² = 4.18) indicó una importante componente espacial en la variabilidad de la temperatura.

El alcance práctico estimado fue de aproximadamente 0.00265 grados, equivalente a cerca de 294 metros, lo que indica que las observaciones mantienen correlación espacial hasta distancias cercanas a dicho valor. En consecuencia, árboles ubicados dentro de este rango tienden a experimentar condiciones térmicas similares, mientras que la dependencia espacial disminuye progresivamente para distancias mayores.

La validación cruzada permitió evaluar la capacidad predictiva del modelo seleccionado. Los valores obtenidos para el MAE (0.8510), RMSE (1.0714), correlación y coeficiente de determinación (0.6393) evidenciaron la capacidad del modelo para reproducir adecuadamente las temperaturas observadas y respaldan la utilización del kriging ordinario como herramienta de interpolación espacial en este contexto.

La superficie generada mediante kriging ordinario permitió identificar gradientes térmicos dentro de la finca, revelando la existencia de microclimas locales que podrían estar asociados a diferencias de altitud, exposición solar, circulación de aire u otras características propias del terreno. Este resultado demuestra que la temperatura no se distribuye de forma homogénea dentro del cultivo.

Por su parte, el mapa de incertidumbre mostró el comportamiento esperado de la varianza de predicción, con menores errores en las zonas con mayor densidad de observaciones y mayores niveles de incertidumbre hacia los bordes del área de estudio. Este patrón aporta confianza sobre la calidad de las estimaciones realizadas en la mayor parte de la finca.

En conjunto, los resultados obtenidos demuestran que la geoestadística constituye una herramienta eficaz para caracterizar la variabilidad espacial de variables microclimáticas en sistemas agrícolas. La identificación de zonas con diferentes condiciones térmicas puede contribuir al diseño de estrategias de manejo localizado, optimización del riego, monitoreo fitosanitario y toma de decisiones orientadas a mejorar la productividad y sostenibilidad del cultivo de aguacate.

Anexo

Análisis Automatizado Multivariable (Variables Climáticas)

A continuación, se automatiza el flujo geoestadístico completo para mapear y evaluar el resto de variables microclimáticas críticas registradas el 01/10/2020. Los resultados se organizan en pestañas para facilitar su comparación.

library(dplyr)
library(geoR)
library(raster)
library(ggplot2)
library(paletteer)
library(stringr)

# 1. Definir la lista de variables climáticas a procesar
variables_climaticas <- c(
  "Temperature",
  "Psychro_Wet_Bulb_Temperature",
  "Relative_Humidity",
  "Heat_Stress_Index",
  "Dew_Point",
  "Wind_Chill"
)

# 2. Malla espacial base (ajustada a 80x80 para balancear velocidad y resolución)
grid_base <- expand.grid(
  lon = seq(min(Datos_Completos_Aguacate$Longitude), max(Datos_Completos_Aguacate$Longitude), length = 80),
  lat = seq(min(Datos_Completos_Aguacate$Latitude), max(Datos_Completos_Aguacate$Latitude), length = 80)
)

# Distancia máxima para los límites del variograma
max.dist <- max(dist(Datos_Completos_Aguacate[,c("Longitude", "Latitude")]))

# 3. Lista de modelos teóricos a evaluar de forma iterativa
modelos_pool <- c(
  "exponential", "gaussian", "spherical", "matern", "circular",
  "cubic", "pentaspherical", "sinepower", "powered.exponential",
  "gneiting", "cauchy", "gencauchy", "whittle"
)

# 4. FUNCIÓN MAESTRA: Ejecuta todo el flujo para una variable específica
procesar_variable_espacial <- function(nombre_var) {

  cat("\n### Variable: ", nombre_var, " {.tabset}\n")

  # a. Crear objeto geoestadístico dinámico
  geo_dinamico <- as.geodata(
    Datos_Completos_Aguacate,
    coords.col = c("Longitude", "Latitude"),
    data.col = nombre_var
  )

  # b. Calcular variograma empírico de forma silenciosa
  v_empirico <- suppressMessages(
    variog(geo_dinamico, option = "bin", uvec = seq(0, max.dist/2, length = 10), messages = FALSE)
  )

  # c. Definir valores iniciales dinámicos basados en la varianza real de la variable actual
  var_actual <- var(Datos_Completos_Aguacate[[nombre_var]], na.rm = TRUE)
  ini.vals_dinamico <- expand.grid(
    sigmasq = seq(var_actual * 0.3, var_actual * 1.8, length = 5),
    phi = seq(max.dist / 50, max.dist / 2, length = 5)
  )

  # d. Ajuste iterativo de los 13 modelos teóricos
  ajustes <- list()
  for (mod in modelos_pool) {
    ajustes[[mod]] <- tryCatch({
      suppressMessages(
        variofit(v_empirico, ini = ini.vals_dinamico, cov.model = mod, wei = "npair", min = "optim", messages = FALSE)
      )
    }, error = function(e) { NULL })
  }
  ajustes <- ajustes[!sapply(ajustes, is.null)]

  if(length(ajustes) == 0) {
    cat("\nNo se pudo ajustar ningun modelo para la variable:", nombre_var, "\n")
    return(NULL)
  }

  # e. Selección del mejor modelo según el menor valor de pérdida ópitma
  vals_perdida <- sapply(ajustes, function(m) m$value)
  lbl_mejor <- names(vals_perdida)[which.min(vals_perdida)]
  mod_optimo <- ajustes[[lbl_mejor]]

  # f. Validación Cruzada protegida y silenciada
  val_cruzada <- tryCatch({
    suppressMessages(xvalid(geodata = geo_dinamico, model = mod_optimo, messages = FALSE))
  }, error = function(e) { NULL })

  # g. Kriging Ordinario continuo (CORREGIDO: Sin el argumento inválido messages = FALSE)
  krig_dinamico <- suppressMessages(
    krige.conv(
      geo_dinamico, loc = grid_base,
      krige = krige.control(
        cov.model = mod_optimo$cov.model,
        cov.pars = mod_optimo$cov.pars,
        nugget = mod_optimo$nugget,
        trend.d = "cte", trend.l = "cte"
      )
    )
  )

  # h. Crear Raster y aplicar la máscara geográfica de la finca
  r_pred <- rasterFromXYZ(cbind(grid_base, krig_dinamico$predict))
  r_pred_finca <- mask(r_pred, finca)
  df_mapa <- as.data.frame(r_pred_finca, xy = TRUE)
  colnames(df_mapa) <- c("Longitude", "Latitude", "Prediccion")
  df_mapa <- filter(df_mapa, !is.na(Prediccion))

  # ---- SUBPESTAÑA 1: MAPA DE PREDICCIÓN ----
  cat("\n#### Mapa de Prediccion\n")

  # Ajuste de paleta de color según la naturaleza física de la variable
  paleta_seleccionada <- if(stringr::str_detect(nombre_var, "Humidity|Dew")) "viridis::viridis" else "ggthemes::Temperature Diverging"

  p1 <- ggplot() +
    geom_raster(data = df_mapa, aes(x = Longitude, y = Latitude, fill = Prediccion)) +
    scale_fill_paletteer_c(paleta_seleccionada) +
    geom_point(data = Datos_Completos_Aguacate, aes(x = Longitude, y = Latitude), color = "black", size = 0.6, alpha = 0.4) +
    labs(title = paste("Prediccion Espacial:", nombre_var), fill = "Unidades") +
    theme_minimal() +
    coord_fixed()

  print(p1)
  cat("\n")

  # ---- SUBPESTAÑA 2: DIAGNÓSTICO ----
  cat("\n#### Diagnostico del Modelo\n")
  cat("**Estructura Matematica Optima Seleccionada:** ", toupper(as.character(lbl_mejor)), "\n\n")

  if(!is.null(val_cruzada)) {
    mae_v <- mean(abs(val_cruzada$error), na.rm = TRUE)
    rmse_v <- sqrt(mean(val_cruzada$error^2, na.rm = TRUE))
    cor_v <- cor(val_cruzada$data, val_cruzada$predicted, use = "complete.obs")

    df_res_met <- data.frame(
      Metrica = c("MAE (Error Absoluto Medio)", "RMSE (Raiz Error Cuadratico Medio)", "Correlacion (r)"),
      Valor = round(c(mae_v, rmse_v, cor_v), 4)
    )
    print(knitr::kable(df_res_met, format = "markdown"))
  } else {
    cat("La validacion cruzada no logro ejecutarse para esta variable.\n")
  }
  cat("\n")
}

# 5. Iterar de forma limpia sobre el vector de variables climáticas
for (variable in variables_climaticas) {
  procesar_variable_espacial(variable)
}