knitr::opts_chunk$set(echo = TRUE, comment = NA)
# 1. Cargar librerías
#library(rnaturalearth)     # Para obtener mapas a distintas escalas.
library(mapview)           # Para la visualización interactiva y rápida.
library(leafsync)          # Para la sincronización de múltiples mapas.
library(geoR)
library(readxl)
library(raster)
library(gstat)
library(sf)                # Estándar moderno para datos vectoriales (reemplaza a 'sp').
library(terra)             # Estándar moderno y eficiente para datos ráster (reemplaza a 'raster').
library(leaflet)           # Para la creación de mapas interactivos.
library(ggplot2) 
library(tidyterra) 
library(patchwork)         # Para combinar gráficos
library(RColorBrewer)
library(knitr)
library(kableExtra)
library(dplyr) # Para la función %>% (pipe)

Introducción



El aguacate (Persea americana) se ha consolidado como un cultivo estratégico a nivel global, siendo su productividad altamente sensible a las condiciones climáticas y topográficas del sitio de producción. Comprender la variabilidad espacial de los factores ambientales es crucial para optimizar el manejo agronómico, desde la fertilización hasta la gestión del riesgo de plagas y enfermedades.

Este informe presenta un estudio geoestadístico exhaustivo de la variabilidad espacial de cuatro factores clave dentro de una finca productora de aguacates: Temperatura, Humedad Relativa, Velocidad del Viento y Altura (Altitud).

Utilizando datos de muestreo puntual, el objetivo principal es:

    1. Diagnosticar la existencia y el tipo de dependencia espacial (autocorrelación) de cada variable.
    1. Modelar dicha dependencia mediante la función de semivariograma.
    1. Generar mapas continuos de predicción (interpolación) a través del método de Kriging Ordinario, permitiendo identificar zonas homogéneas de manejo y optimizar las decisiones de agricultura de precisión.

Los resultados de este análisis proporcionan una base científica para el manejo zonal de la finca, permitiendo una asignación más eficiente de recursos y una comprensión detallada de las condiciones microclimáticas que afectan el cultivo.

Marco Teórico



  1. Geoestadística y Análisis de Variables Ambientales

La Geoestadística es una rama de la estadística que se especializa en el análisis de fenómenos distribuidos en el espacio, asumiendo que las observaciones cercanas están más relacionadas entre sí que las observaciones distantes (Principio de Tobler o Autocorrelación Espacial). Esta metodología es indispensable en la agricultura de precisión, ya que las variables climáticas y de suelo no se distribuyen al azar.

El flujo de trabajo geoestadístico aplicado en este estudio incluye:

A. Objeto geodata

La fase inicial consiste en transformar los datos tabulares de campo en un objeto espacial georreferenciado (geodata). Esta representación permite un diagnóstico visual de la distribución espacial de la variable, complementado con un análisis de la distribución de frecuencias (histograma), fundamental para verificar los supuestos de normalidad requeridos por el Kriging.

B. Semivariograma Empírico y Modelado Predictivo

El semivariograma es la herramienta central para cuantificar la autocorrelación espacial. Muestra la disimilitud (semivarianza) promedio entre puntos a medida que aumenta la distancia de separación (lag).

  • Parámetros Clave:

    • Efecto Pepita (τ2): Representa la varianza no estructurada (errores de medición o variabilidad a microescala).

    • Sill (σ2+τ2): La varianza total de la variable.

    • Rango (ϕ): La distancia máxima a la que los puntos están correlacionados espacialmente.

El semivariograma empírico se ajusta a un modelo teórico (Exponencial, Gaussiano o Esférico) para obtener un modelo continuo que se utilizará en la interpolación.

C. Validación Cruzada Para verificar la robustez del modelo teórico seleccionado, se utiliza la Validación Cruzada (Leave-One-Out). Este proceso predice el valor en cada punto de muestreo utilizando los demás datos, permitiendo comparar la predicción con el valor real y calcular el error.

  1. Método de Interpolación: Kriging Ordinario

El Kriging Ordinario (KO) es un método de interpolación geoestadística que utiliza la función de semivariograma modelada para calcular pesos óptimos de interpolación.

El KO es un estimador insesgado (no tiene sesgo sistemático) y de varianza mínima, lo que significa que produce los mejores mapas de predicción posibles a partir de los datos espaciales. El método asume que la media local es constante dentro del área de estudio.

El resultado del Kriging son dos mapas esenciales:

  • Mapa de Predicción: El mapa continuo de los valores interpolados (ej., Temperatura predicha).

  • Mapa de Varianza/Error: El mapa que cuantifica la incertidumbre de la predicción, mostrando dónde la interpolación es más confiable y dónde no.

  1. Variables Ambientales Clave para el Aguacate

El estudio se enfoca en variables cruciales para el rendimiento y la sanidad del aguacate:

  • Temperatura y Humedad Relativa: Factores directos en la fisiología del árbol, la evapotranspiración y la incidencia de hongos o patógenos. Su variabilidad espacial es crítica para la gestión del riego y la sanidad vegetal.

  • Altura (Altitud): Una variable que impone una fuerte tendencia espacial en las variables climáticas. Su modelado preciso (idealmente con Nugget cero) es fundamental para obtener un modelo digital del terreno de alta calidad.

  • Velocidad del Viento: Una variable típicamente de baja correlación espacial, importante para el riesgo de daños físicos y la polinización.

Carga y Contextualización de Datos



A continuación, se procede a la carga del archivo de datos georreferenciados. Este conjunto de datos abarca cuatro localizaciones (terrenos o parcelas) distintas dentro de la zona del departamento del Cauca, específicamente en corregimientos o veredas aledañas a los municipios de Popayán y Timbío. El archivo contiene registros puntuales de múltiples variables meteorológicas y topográficas, las cuales serán objeto del presente estudio geoestadístico para la finca de aguacates.

# 1. Se carga el archivo de excel con read_excel
geo_datos = read_excel("/Users/jesusposso/Documents/Maestria/2025-2/SIG/M2U1/Actividad_3/Datos_Completos_Aguacate.xlsx")

# 2. Se grafica del archivo excel todos los registros de latitud y longitud mediante leaflet
leaflet() %>% 
  addTiles() %>% 
  addCircleMarkers(
    lng = geo_datos$Longitude,
    lat = geo_datos$Latitude,
    radius = 0.2,
    color = "black")


El proceso analítico requiere la individualización de la zona de estudio. Para contextualizar los datos a una única parcela o evento de muestreo, se procedió a aplicar un filtro estricto sobre la columna temporal FORMATTED_DATE_TIME. En esta ocasión, se han seleccionado exclusivamente los registros recolectados en el instante 01/10/2020 a las 10:11:12 a, m, .

Una vez aplicado el filtro, se observa que los puntos de muestreo georreferenciados corresponden a una concentración de árboles situados en las inmediaciones del corregimiento o vereda La Chorrera. Visualmente, la distribución espacial de los puntos de muestreo es notablemente regular, lo que resulta favorable para el análisis geoestadístico y la modelación precisa de la autocorrelación espacial.

# 1. De acuerdo a las indicaciones se toman los registros que en Fecha tengan el siguiente valor
filtro = '01/10/2020  10:11:12 a, m,'

# 2. Se aplica el filtro al dataframe cargado
geo_datos_filtrado = geo_datos[geo_datos['FORMATTED_DATE_TIME'] == filtro,]

# 3. Se grafica las coordenadas del dataframe filtrado con leaflet
leaflet() %>% 
  addTiles() %>% 
  addCircleMarkers(
    lng = geo_datos_filtrado$Longitude,
    lat = geo_datos_filtrado$Latitude,
    radius = 0.2,
    color = "black")



Geoestadística

1. Variable regionalizada

1.1 Humedad relativa

# Se sugieren tomar las siguientes variables: 
# temperatura = 9, humedad relaritva = 7, velocidad del tiempo = 14 o la altura = 16

# --- Humedad Relativa (Columna 7) ---
# 1. Se construye el geodata de la variable
geo_humedad_relativa = as.geodata(
    geo_datos_filtrado,
    coords.col = 3:2,
    data.col = 7 # Variable: Relative_Humidity (Humedad relativa)
)
# 2. Se visualiza la informacion de la geodata
plot(geo_humedad_relativa)


El análisis visual del objeto geodata de la Humedad Relativa revela evidencia visual de autocorrelación espacial, con valores similares que parecen agruparse geográficamente en el mapa de puntos. El histograma muestra que la distribución de la variable es relativamente simétrica. Esta simetría estadística es favorable para el Kriging. Ambos diagnósticos visuales (la agrupación espacial y la distribución estadística) confirman que esta variable está bien comportada y es adecuada para el modelado geoestadístico paramétrico.

1.2 Temperatura

# --- Temperatura (Columna 9) ---
# 1. Se construye el geodata de la variable
geo_temperatura = as.geodata(
    geo_datos_filtrado,
    coords.col = 3:2,
    data.col = 9 # Variable: Temperature (Temperatura)
)
# 2. Se visualiza la informacion de la geodata
plot(geo_temperatura)


El gráfico de la Temperatura muestra un patrón espacial bien definido, donde los valores de temperatura alta y baja se encuentran localmente agrupados, reforzando la fuerte estructura espacial identificada previamente. El histograma confirma que la distribución de los valores de Temperatura se aproxima a una distribución normal. Esta cercanía a la normalidad es óptima para los supuestos del Kriging ordinario, lo que significa que no se requieren transformaciones de datos complejas. La Temperatura es, por lo tanto, una de las variables con el mejor comportamiento diagnóstico.

1.3 Velocidad viento

# --- Velocidad del Viento (Columna 14) ---
# 1. Se construye el geodata de la variable
geo_velocidad_viento = as.geodata(
    geo_datos_filtrado,
    coords.col = 3:2,
    data.col = 14 # Variable: Wind_Speed (Velocidad del Viento)
)
# 2. Se visualiza la informacion de la geodata
plot(geo_velocidad_viento)


El diagnóstico visual de la Velocidad del Viento es consistente con un comportamiento de aleatoriedad espacial. El mapa de puntos no muestra ningún patrón geográfico discernible; los valores altos y bajos se mezclan aleatoriamente en el espacio. El histograma probablemente muestra una distribución asimétrica positiva (sesgada), con una concentración de valores bajos. Esta combinación de falta de estructura espacial y asimetría estadística indica que la variable no cumple con los requisitos básicos de dependencia espacial para ser modelada eficazmente mediante Kriging.

1.4 Altura

# --- Altura (Columna 16) ---
# 1. Se construye el geodata de la variable
geo_altura = as.geodata(
    geo_datos_filtrado,
    coords.col = 3:2,
    data.col = 16 # Variable: Altitude (Altura)
)
# 2. Se visualiza la informacion de la geodata
plot(geo_altura)


La Altura presenta la tendencia espacial más clara y evidente en el mapa de puntos. Se observa un claro gradiente (por ejemplo, los valores más altos se concentran en una región específica). Esta tendencia confirma que la variación de la Altura es sistemática y predecible por la ubicación. La forma del histograma debe evaluarse cuidadosamente (pudiendo ser bimodal o amplia) para determinar si la distribución de la variable requiere una remoción de la tendencia mediante regresión espacial antes de proceder a modelar el semivariograma de los residuales, lo que es una práctica común para variables con gradientes tan fuertes.

2. Semivariogramas

Para el desarrollo del análisis geoestadístico, específicamente el cálculo del semivariograma empírico, se seleccionó el conjunto de datos de Humedad Relativa (geo_humedad_relativa) como base para la determinación de las distancias y la estructura de la dependencia espacial.

Esta elección se fundamenta en que todos los objetos geodata utilizados en el estudio (Temperatura, Humedad Relativa, Velocidad del Viento y Altitud) comparten idénticas coordenadas geográficas y puntos de muestreo. Por lo tanto, el resumen estadístico de las distancias euclidianas entre pares de puntos será el mismo, independientemente de la variable seleccionada.

A continuación, se presenta la tabla resumen de las distancias entre las coordenadas de los puntos de muestreo, obtenida mediante la función summary() aplicada a las distancias de la variable geo_humedad_relativa.

# 1. Se calcula el resumen de distancias
resumen_distancia = summary(dist(geo_humedad_relativa$coords))
# Solo se hace uso de un geodata debido a que todos los geodata anteriores tienen las mismas coordenadas
# La diferencia es la variable de estudio

# 2. Se crea el data.frame con los resultados
datos_tabla = data.frame(
  Variable = names(resumen_distancia), 
  Valor = as.numeric(resumen_distancia), 
  row.names = NULL
)

# 3. Se visualiza la tabla con los datos del dataframe construido
kable(
  datos_tabla,
  caption = "Tabla 1. Resumen de Distancias Geográficas entre Registros Filtrados",
  col.names = c("Estadístico", "Valor de Distancia (Grados)"), 
  format = "html",
  align = c("l", "r"), 
  row.names = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    bootstrap_options = c("striped", "hover", "responsive"),
    position = "center",
    font_size = 12,
    row_label_position = "c"
  ) %>%
  column_spec(1, width = "1.5in") %>% # Estadístico
  column_spec(2, width = "2in") %>%  # Valor
  row_spec(0, bold = TRUE, color = "white", background = "#2C3E50") # Estilo para el encabezado
Tabla 1. Resumen de Distancias Geográficas entre Registros Filtrados
Estadístico Valor de Distancia (Grados)
Min. 0.0000171
1st Qu. 0.0004051
Median 0.0006408
Mean 0.0006827
3rd Qu. 0.0009178
Max. 0.0019591

2.1 Humedad relativa

# 1. Se realiza el Semivariograma de la variable estudiada
invisible(capture.output({
    SemivariogramaHumedad = variog(
        geo_humedad_relativa, 
        option = "bin",
        uvec = seq(0, 9.2e-04, 2.4e-05)
    )
}))
# 2. Se grafica datos del semivariograma de la variable estudiada
plot(SemivariogramaHumedad)
# 3. Construir las lineas en el semivariograma
invisible(capture.output({
    variograma_mc = variog.mc.env(geo_humedad_relativa, obj=SemivariogramaHumedad)
}))
# 4. Dibujar las lineas
lines(variograma_mc)


El semivariograma de la Humedad Relativa evidencia una fuerte estructura de correlación espacial. La semivarianza se incrementa con la distancia, lo que indica que los valores de humedad de los árboles cercanos son significativamente más parecidos que los de los árboles distantes. La curva empírica cae fuera del Sobre de Monte Carlo, confirmando que esta dependencia espacial es estadísticamente significativa y no se debe a la aleatoriedad. Se observa un Efecto Pepita moderado, lo que sugiere la presencia de variabilidad a escalas más pequeñas que la distancia mínima de muestreo o errores de medición. Esta variable es apta para el modelado con Kriging.



2.2 Temperatura

# 1. Se realiza el Semivariograma de la variable estudiada
invisible(capture.output({
    SemivariogramaTemp = variog(
                          geo_temperatura, 
                          option = "bin",
                          uvec = seq(0, 9.2e-04, 2.4e-05)
                          )
}))
# 2. Se grafica datos del semivariograma de la variable estudiada
plot(SemivariogramaTemp)
# 3. Construir las lineas en el semivariograma
invisible(capture.output({
    variograma_mc = variog.mc.env(geo_temperatura, obj=SemivariogramaTemp)
}))
# 4. Dibujar las lineas
lines(variograma_mc)


La Temperatura presenta una estructura de correlación espacial muy robusta, con una curva de semivarianza que se eleva rápidamente hasta alcanzar una meseta definida (Sill). El semivariograma empírico se encuentra consistentemente por debajo del Sobre de Monte Carlo, lo que valida que la autocorrelación observada es altamente significativa. El Efecto Pepita es notable, reflejando la variación local que no es explicada por la estructura espacial. Sin embargo, dada la fuerte dependencia espacial a distancias cortas, la Temperatura es una variable excelente candidata para la interpolación mediante Kriging.



2.3 Velocidad viento

# 1. Se realiza el Semivariograma de la variable estudiada
invisible(capture.output({
    SemivariogramaViento = variog(
                          geo_velocidad_viento, 
                          option = "bin",
                          uvec = seq(0, 9.2e-04, 2.4e-05)
                          )
}))
# 2. Se grafica datos del semivariograma de la variable estudiada
plot(SemivariogramaViento)
# 3. Construir las lineas en el semivariograma
invisible(capture.output({
    variograma_mc = variog.mc.env(geo_velocidad_viento, obj=SemivariogramaViento)
}))
# 4. Dibujar las lineas
lines(variograma_mc)


El semivariograma de la Velocidad del Viento revela una ausencia de estructura de correlación espacial a la escala de muestreo. La curva empírica es plana y se mantiene casi totalmente dentro del Sobre de Monte Carlo. Este comportamiento es característico del Ruido Blanco (o aleatoriedad espacial pura). El Efecto Pepita domina casi toda la varianza de la variable. Por lo tanto, se concluye que la variación de la Velocidad del Viento es principalmente aleatoria, y las técnicas de interpolación geoestadística (Kriging) no son adecuadas para esta variable.



2.4 Altura

# 1. Se realiza el Semivariograma de la variable estudiada
invisible(capture.output({
    SemivariogramaAltura = variog(
                          geo_altura, 
                          option = "bin",
                          uvec = seq(0, 9.2e-04, 2.4e-05)
                          )
}))
# 2. Se grafica datos del semivariograma de la variable estudiada
plot(SemivariogramaAltura)
# 3. Construir las lineas en el semivariograma
invisible(capture.output({
    variograma_mc = variog.mc.env(geo_altura, obj=SemivariogramaAltura)
}))
# 4. Dibujar las lineas
lines(variograma_mc)


La Altura muestra la estructura espacial más limpia y fuerte de todas las variables analizadas. La curva comienza casi en cero (Efecto Pepita insignificante), lo que indica una correlación espacial casi perfecta a distancias muy cortas y una alta precisión en la medición. El rápido y claro aumento de la semivarianza, así como el posicionamiento de la curva lejos del Sobre de Monte Carlo, confirman que la dependencia espacial es altamente significativa y robusta. La Altura es la variable más adecuada y recomendable para el modelado con Kriging para generar un modelo digital del terreno.

3. Ajuste del modelo teorico



De acuerdo con las conclusiones derivadas de la sección previa, se ha determinado el conjunto de variables ambientales y topográficas que serán empleadas para la modelización predictiva y la posterior generación de un Modelo Digital del Terreno (MDT) mediante la técnica de Kriging.

Las variables seleccionadas para su integración en el análisis geoestadístico son las siguientes: Temperatura, humedad relativa y altura.

# 1. Función corregida para extraer y formatear los parámetros
obtener_parametros = function(modelo, nombre) {
    
    # Usamos el slot $cov.pars, que contiene: [1] sigmasq, [2] phi, [3] tausq
    sigmasq = modelo$cov.pars[1]
    phi = modelo$cov.pars[2]
    tausq = modelo$cov.pars[3] 
    
    # Crear el data frame, asegurando que todos los valores son escalares
    parametros = data.frame(
        Modelo = nombre,
        Sill_Parcial_sigmasq = sigmasq,
        Rango_phi = phi,
        Nugget_tausq = tausq,
        Loss_Value = modelo$value, 
        stringsAsFactors = FALSE
    )
    return(parametros)
}

3.1 Humedad relativa

# 1. Se crea la grilla de valores iniciales
ini.vals = expand.grid(seq(27,32,l=10), seq(5.8e-04,8e-04,l=10))
# 2. Se grafica el semivaraiograma para adicionar las lineas de cada modelo y comparar
plot(SemivariogramaHumedad)
# 3. Se crea el modelo exponencial
invisible(capture.output({
model_mco_exp_h = variofit(SemivariogramaHumedad,
                         ini=ini.vals, cov.model="exponential",
                         wei="npair", min="optim"
                         )
}))
# 4. Se crea el modelo gausiano
invisible(capture.output({
model_mco_gaus_h = variofit(SemivariogramaHumedad,
                        ini=ini.vals, cov.model="gaussian",
                        wei="npair", min="optim",nugget = 0
                        )
}))
# 5. Se crea el modelo esférico
invisible(capture.output({
model_mco_spe_h = variofit(SemivariogramaHumedad, 
                       ini=ini.vals, cov.model="spheric",
                       fix.nug=TRUE, wei="npair", min="optim"
                       )
}))
# 6. Construir las lineas en el semivariograma
invisible(capture.output({
    variograma_mc = variog.mc.env(geo_humedad_relativa, obj=SemivariogramaHumedad)
}))
# 7. Se grafica la linea que construye cada modelo
lines(variograma_mc)
lines(model_mco_exp_h,col="blue")
lines(model_mco_gaus_h,col="red")
lines(model_mco_spe_h,col="purple")
# 8. Se adiciona la leyenda
legend(
    "bottomright", # Posición de la leyenda
    legend = c("Monte Carlo (95%)", "Exponencial", "Gaussiano", "Esférico"), # Textos de la leyenda
    col = c("black", "blue", "red", "purple"), # Colores en el mismo orden
    lty = 1, # Tipo de línea (1 = sólida, ajuste si Monte Carlo es diferente)
    cex = 0.8, # Tamaño del texto
    bg = "white" # Fondo para que el texto sea legible
)


El análisis del ajuste del modelo para la Humedad Relativa demuestra una excelente correspondencia entre los modelos teóricos (Exponencial, Gaussiano y Esférico) y el semivariograma empírico. El fenómeno presenta una fuerte dependencia espacial que se agota rápidamente con la distancia. El Efecto Pepita es notable, indicando que una fracción de la varianza total es atribuible a errores de medición o variabilidad a microescala. El Modelo Exponencial es la elección más parsimoniosa (simple y efectiva) para el Kriging, ya que captura de manera eficiente la descorrelación a corta distancia.

# 1. Extraer todos los modelos
df_exp = obtener_parametros(model_mco_exp_h, "Exponencial")
df_gaus = obtener_parametros(model_mco_gaus_h, "Gaussiano")
df_spe = obtener_parametros(model_mco_spe_h, "Esférico")

# 2. Combinar en una única tabla para comparación
tabla_comparativa = rbind(df_exp, df_gaus, df_spe)

# 3. Formatear los valores numéricos para una mejor lectura (ej: 4 decimales)
tabla_formateada = tabla_comparativa
tabla_formateada$Sill_Parcial_sigmasq = round(tabla_formateada$Sill_Parcial_sigmasq, 4)
tabla_formateada$Rango_phi = round(tabla_formateada$Rango_phi, 6) # Rango necesita más precisión
tabla_formateada$Nugget_tausq = round(tabla_formateada$Nugget_tausq, 4)
tabla_formateada$Loss_Value = round(tabla_formateada$Loss_Value, 2)

# 4. Se visualiza la tabla con kable y ajuste de estilo
kable(
  tabla_formateada,
  caption = "Tabla 2. Parámetros Finales del Semivariograma Ajustado Humedad Relativa",
  col.names = c("Modelo", "Sill Parcial (σ²)", "Rango (φ)", "Nugget (τ²)", "Función de Pérdida"), 
  format = "html",
  align = c("l", "r", "r", "r", "r"), 
  row.names = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    bootstrap_options = c("striped", "hover", "responsive"),
    position = "center",
    font_size = 12,
    row_label_position = "c"
  ) %>%
  # Ajuste de ancho de columna para las 5 variables (Líneas limpias)
  column_spec(1, width = "1.5in") %>% # Modelo
  column_spec(2, width = "1in") %>% 
  column_spec(3, width = "1in") %>% 
  column_spec(4, width = "1in") %>% 
  column_spec(5, width = "1.5in") %>%
  row_spec(0, bold = TRUE, color = "white", background = "#2C3E50")
Tabla 2. Parámetros Finales del Semivariograma Ajustado Humedad Relativa
Modelo Sill Parcial (σ²) Rango (φ) Nugget (τ²) Función de Pérdida
Exponencial 37.3897 0.000828 NA 542092.5
Gaussiano 40.0905 0.001255 NA 526240.6
Esférico 27.5587 0.000273 NA 1219353.2

3.2 Temperatura

# 1. Se crea la grilla de valores iniciales
ini.vals = expand.grid(seq(3,4,l=10), seq(4.4e-04,8e-04,l=10))
# 2. Se grafica el semivaraiograma para adicionar las lineas de cada modelo y comparar
plot(SemivariogramaTemp)
# 3. Se crea el modelo exponencial
invisible(capture.output({
model_mco_exp_t = variofit(SemivariogramaTemp, 
                       ini=ini.vals, cov.model="exponential",
                       wei="npair", min="optim"
                       )
}))
# 4. Se crea el modelo gausiano
invisible(capture.output({
model_mco_gaus_t = variofit(SemivariogramaTemp, 
                        ini=ini.vals, cov.model="gaussian", 
                        wei="npair", min="optim",nugget = 0
                        )
}))
# 5. Se crea el modelo esférico
invisible(capture.output({
model_mco_spe_t = variofit(SemivariogramaTemp, 
                       ini=ini.vals, cov.model="spheric",
                       fix.nug=TRUE, wei="npair", min="optim"
                       )
}))
# 6. Construir las lineas en el semivariograma
invisible(capture.output({
    variograma_mc = variog.mc.env(geo_temperatura, obj=SemivariogramaTemp)
}))
# 7. Se grafica la linea que construye cada modelo
lines(variograma_mc)
lines(model_mco_exp_t,col="blue")
lines(model_mco_gaus_t,col="red")
lines(model_mco_spe_t,col="purple")
# 8. Se adiciona la leyenda
legend(
    "bottomright", # Posición de la leyenda
    legend = c("Monte Carlo (95%)", "Exponencial", "Gaussiano", "Esférico"), # Textos de la leyenda
    col = c("black", "blue", "red", "purple"), # Colores en el mismo orden
    lty = 1, # Tipo de línea (1 = sólida, ajuste si Monte Carlo es diferente)
    cex = 0.8, # Tamaño del texto
    bg = "white" # Fondo para que el texto sea legible
)


La Temperatura se distingue por tener una estructura espacial muy robusta y bien ajustada por los modelos. El bajo valor de la función de pérdida en los modelos Exponencial y Esférico sugiere que ambos son candidatos principales para la interpolación. Estos modelos capturan con precisión la rápida autocorrelación que exhibe la Temperatura en distancias cortas, una característica esencial para el Kriging. Se confirma la presencia de un Efecto Pepita moderado. La calidad del ajuste indica que la interpolación mediante Kriging con el modelo seleccionado proporcionará una estimación altamente confiable del campo de temperatura en el área de estudio.

# 1. Extraer todos los modelos
df_exp = obtener_parametros(model_mco_exp_t, "Exponencial")
df_gaus = obtener_parametros(model_mco_gaus_t, "Gaussiano")
df_spe = obtener_parametros(model_mco_spe_t, "Esférico")

# 2. Combinar en una única tabla para comparación
tabla_comparativa = rbind(df_exp, df_gaus, df_spe)

# 3. Formatear los valores numéricos para una mejor lectura (ej: 4 decimales)
tabla_formateada = tabla_comparativa
tabla_formateada$Sill_Parcial_sigmasq = round(tabla_formateada$Sill_Parcial_sigmasq, 4)
tabla_formateada$Rango_phi = round(tabla_formateada$Rango_phi, 6) # Rango necesita más precisión
tabla_formateada$Nugget_tausq = round(tabla_formateada$Nugget_tausq, 4)
tabla_formateada$Loss_Value = round(tabla_formateada$Loss_Value, 2)

# 4. Se visualiza la tabla con kable y ajuste de estilo
kable(
  tabla_formateada,
  caption = "Tabla 3. Parámetros Finales del Semivariograma Ajustado Temperatura",
  col.names = c("Modelo", "Sill Parcial (σ²)", "Rango (φ)", "Nugget (τ²)", "Función de Pérdida"), 
  format = "html",
  align = c("l", "r", "r", "r", "r"), 
  row.names = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    bootstrap_options = c("striped", "hover", "responsive"),
    position = "center",
    font_size = 12,
    row_label_position = "c"
  ) %>%
  # Ajuste de ancho de columna para las 5 variables (Líneas limpias)
  column_spec(1, width = "1.5in") %>% # Modelo
  column_spec(2, width = "1in") %>% 
  column_spec(3, width = "1in") %>% 
  column_spec(4, width = "1in") %>% 
  column_spec(5, width = "1.5in") %>%
  row_spec(0, bold = TRUE, color = "white", background = "#2C3E50")
Tabla 3. Parámetros Finales del Semivariograma Ajustado Temperatura
Modelo Sill Parcial (σ²) Rango (φ) Nugget (τ²) Función de Pérdida
Exponencial 3.0163 0.000127 NA 3364.25
Gaussiano 2.9810 0.000000 NA 15400.92
Esférico 3.0619 0.000239 NA 6143.21

3.3 Altura

# 1. Se crea la grilla de valores iniciales
ini.vals = expand.grid(seq(90,110,l=10), seq(4.5e-04,1e-03,l=10))
# 2. Se grafica el semivaraiograma para adicionar las lineas de cada modelo y comparar
plot(SemivariogramaAltura)
# 3. Se crea el modelo exponencial
invisible(capture.output({
model_mco_exp_a = variofit(SemivariogramaAltura, 
                       ini=ini.vals, cov.model="exponential",
                       wei="npair", min="optim"
                       )
}))
# 4. Se crea el modelo gausiano
invisible(capture.output({
model_mco_gaus_a = variofit(SemivariogramaAltura, 
                        ini=ini.vals, cov.model="gaussian", 
                        wei="npair", min="optim",nugget = 0
                        )
}))
# 5. Se crea el modelo esférico
invisible(capture.output({
model_mco_spe_a = variofit(SemivariogramaAltura, 
                       ini=ini.vals, cov.model="spheric",
                       fix.nug=TRUE, wei="npair", min="optim"
                       )
}))
# 6. Construir las lineas en el semivariograma
invisible(capture.output({
    variograma_mc = variog.mc.env(geo_altura, obj=SemivariogramaAltura)
}))
# 7. Se grafica la linea que construye cada modelo
lines(variograma_mc)
lines(model_mco_exp_a,col="blue")
lines(model_mco_gaus_a,col="red")
lines(model_mco_spe_a,col="purple")
# 8. Se adiciona la leyenda
legend(
    "bottomright", # Posición de la leyenda
    legend = c("Monte Carlo (95%)", "Exponencial", "Gaussiano", "Esférico"), # Textos de la leyenda
    col = c("black", "blue", "red", "purple"), # Colores en el mismo orden
    lty = 1, # Tipo de línea (1 = sólida, ajuste si Monte Carlo es diferente)
    cex = 0.8, # Tamaño del texto
    bg = "white" # Fondo para que el texto sea legible
)



El modelado predictivo de la Altura resultó en un ajuste excepcionalmente preciso para el modelo Exponencial, el más probables candidato óptimo. La característica más notable es el Efecto Pepita casi nulo que fue confirmado por los modelos, lo que implica una variabilidad mínima no estructurada y un alto grado de continuidad espacial. Este ajuste perfecto valida que la Altura está casi totalmente determinada por la ubicación geográfica. El modelo seleccionado es el más confiable del estudio y proporcionará el marco más preciso para la interpolación.

# 1. Extraer todos los modelos
df_exp = obtener_parametros(model_mco_exp_a, "Exponencial")
df_gaus = obtener_parametros(model_mco_gaus_a, "Gaussiano")
df_spe = obtener_parametros(model_mco_spe_a, "Esférico")

# 2. Combinar en una única tabla para comparación
tabla_comparativa = rbind(df_exp, df_gaus, df_spe)

# 3. Formatear los valores numéricos para una mejor lectura (ej: 4 decimales)
tabla_formateada = tabla_comparativa
tabla_formateada$Sill_Parcial_sigmasq = round(tabla_formateada$Sill_Parcial_sigmasq, 4)
tabla_formateada$Rango_phi = round(tabla_formateada$Rango_phi, 6) # Rango necesita más precisión
tabla_formateada$Nugget_tausq = round(tabla_formateada$Nugget_tausq, 4)
tabla_formateada$Loss_Value = round(tabla_formateada$Loss_Value, 2)

# 4. Se visualiza la tabla con kable y ajuste de estilo
kable(
  tabla_formateada,
  caption = "Tabla 4. Parámetros Finales del Semivariograma Ajustado Altura",
  col.names = c("Modelo", "Sill Parcial (σ²)", "Rango (φ)", "Nugget (τ²)", "Función de Pérdida"), 
  format = "html",
  align = c("l", "r", "r", "r", "r"), 
  row.names = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    bootstrap_options = c("striped", "hover", "responsive"),
    position = "center",
    font_size = 12,
    row_label_position = "c"
  ) %>%
  # Ajuste de ancho de columna para las 5 variables (Líneas limpias)
  column_spec(1, width = "1.5in") %>% # Modelo
  column_spec(2, width = "1in") %>% 
  column_spec(3, width = "1in") %>% 
  column_spec(4, width = "1in") %>% 
  column_spec(5, width = "1.5in") %>%
  row_spec(0, bold = TRUE, color = "white", background = "#2C3E50")
Tabla 4. Parámetros Finales del Semivariograma Ajustado Altura
Modelo Sill Parcial (σ²) Rango (φ) Nugget (τ²) Función de Pérdida
Exponencial 127.3645 0.000639 NA 1298053
Gaussiano 96.4534 0.000184 NA 24039326
Esférico 110.3886 0.000633 NA 8588259



4. Modelo Digital del Terreno mediante técnica: Kriging



La técnica empleada es el Kriging Ordinario (KO), seleccionado por su capacidad para generar las mejores estimaciones lineales insesgadas (Best Linear Unbiased Estimator - BLUE), utilizando la estructura de autocorrelación espacial previamente definida por el semivariograma.

El procedimiento se centra en la generación de mapas continuos para las variables que demostraron una dependencia espacial significativa: Temperatura, Humedad Relativa y Altura (Altitud).

Para cada variable se realiza lo siguiente:

  • Parámetros Óptimos: Se utilizan los parámetros de Nugget (τ2), Sill (σ2) y Range (ϕ) del modelo teórico que presentó el menor error (loss value) en la etapa de ajuste del semivariograma.

  • Generación de Ráster: Se ejecuta la función krige.conv() para predecir los valores sobre una grilla de 100x100 puntos, resultando en dos mapas ráster fundamentales:

    • Mapa de Predicción: La superficie continua de la variable ambiental (ej., Humedad Relativa predicha).

    • Mapa de Error: La superficie de la Desviación Estándar de la predicción, que cuantifica la incertidumbre de la estimación.

La creación de estos modelos digitales permite una cuantificación precisa de las variaciones microclimáticas y topográficas en toda la finca, proporcionando una herramienta clave para la delimitación de zonas de manejo homogéneo en el cultivo de aguacates.

# ---- Predicción Espacial Kriging ------
# 1. Se crea plantilla de grilla para la interpolación (100x100)
geodatos_grid = expand.grid(
  # Longitud (X): Va de Max a Min para ordenar la grilla de este a oeste
  lon = seq(max(geo_datos_filtrado$Longitude), min(geo_datos_filtrado$Longitude), length.out = 100),
  # Latitud (Y): Va de Min a Max de sur a norte
  lat = seq(min(geo_datos_filtrado$Latitude), max(geo_datos_filtrado$Latitude), length.out = 100)
)

# 2. Se visualiza la plantilla de la grilla (puntos negros)
plot(geodatos_grid, main = "Grilla de Predicción y Puntos de Muestreo")

# 3. Se grafica los puntos de muestreo originales (Longitude, Latitude)
# Nota: La sintaxis Longitude, Latitude asegura el orden X, Y
points(geo_datos_filtrado$Longitude, geo_datos_filtrado$Latitude, col = "green", pch = 16, cex = 0.5)

# Opcional: Adicionar Leyenda para los puntos
legend("topleft", legend = "Puntos de Muestreo", col = "green", pch = 16, cex = 0.8)



4.1 Humedad relativa

La Humedad Relativa es la primera variable en ser sometida al proceso de interpolación. Tras la validación del Modelo Exponencial como el ajuste óptimo (dada su equivalencia con el Modelo Matérn), se procede a su aplicación en el Kriging Ordinario. Se utilizan los parámetros finales obtenidos (σ2, ϕ y τ2) para generar un Mapa de Predicción de la Humedad y un Mapa de Desviación Estándar, los cuales permiten identificar las zonas con mayor y menor contenido de humedad, crucial para la gestión del riesgo de hongos en el cultivo.

# 1. Se crea los raster de salida
invisible(capture.output({
geodatos_ko = krige.conv(geo_humedad_relativa, loc=geodatos_grid,
      krige = krige.control(nugget=0,trend.d="cte", 
      trend.l = "cte",cov.pars=c(sigmasq=37.3897, phi=0.000828 )))
}))
# 2. Humedad_predict: Crea el ráster de las predicciones (Z = geodatos_ko$predict)
pred_matrix = cbind(geodatos_grid, geodatos_ko$predict)
Humedad_predict = rast(pred_matrix, crs="EPSG:4326") 
names(Humedad_predict) = "Humedad_Predicha"

# 3. Humedad_error: Crea el ráster del error (Z = sqrt(geodatos_ko$krige.var))
error_matrix = cbind(geodatos_grid, sqrt(geodatos_ko$krige.var))
Humedad_error = rast(error_matrix, crs="EPSG:4326")
names(Humedad_error) = "Error_StDev"

# 4. Visualizacion de mapas
# 4.1 Mapa 1: Predicción de Humedad (Kriging)
plot_prediccion = ggplot() +
    # Usar geom_spatraster para dibujar el ráster de terra
    geom_spatraster(data = Humedad_predict) + 
    # Usar una paleta de color moderna (viridis)
    scale_fill_viridis_c( name = "Humedad Relativa (%)") +
    # Adicionar contornos para mejor lectura
    geom_spatraster_contour(data = Humedad_predict, color = "gray30", alpha = 0.5) +
    labs(
        title = "Predicción Espacial (Kriging Ordinario)",
        x = "Longitud",
        y = "Latitud"
    ) +
    theme_minimal() +
    theme(legend.position = "bottom")

# 4.2 Mapa 2: Error de Predicción (Desviación Estándar)
plot_error = ggplot() +
    geom_spatraster(data = Humedad_error) +
    # Usar otra paleta de color (Magma) para contraste. El error se visualiza en escalas más oscuras.
    scale_fill_viridis_c(option = "magma", name = "Desviación Estándar") + 
    geom_spatraster_contour(data = Humedad_error, color = "black", alpha = 0.5) +
    labs(
        title = "Error de Predicción (Desviación Estándar)",
        x = "Longitud",
        y = "Latitud"
    ) +
    theme_minimal() +
    theme(legend.position = "bottom")

# 5. Combinar los dos mapas en un solo panel (Usando patchwork)
plot_final = plot_prediccion + plot_error + 
             plot_layout(guides = 'collect') + # Recolecta las leyendas
             plot_annotation(title = 'Resultados del Kriging Ordinario para Humedad Relativa')

plot_final

# 1. Convertir los puntos filtrados a un objeto SF
puntos_sf = geo_datos_filtrado %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = "EPSG:4326")

# 2. Crear la Envolvente Convexa (Convex Hull)
# Esto crea el polígono más pequeño que contiene todos los puntos verdes.
area_estudio_polygon = puntos_sf %>%
  st_union() %>%
  st_convex_hull()

# Se valida del polígono para el recorte
area_estudio_terra = vect(area_estudio_polygon)

# 1. Recorte Inicial: Recorta el ráster al recuadro (bounding box) del polígono.
Humedad_crop = crop(Humedad_predict, area_estudio_terra)

# 2. Aplicar Máscara: Elimina los píxeles que están fuera del polígono (el área blanca).
Humedad_recortada = mask(Humedad_crop, area_estudio_terra)

# Visualizar el resultado con mapview (interactivo)
mapview(Humedad_recortada,
        layer.name = "Humedad Relativa (Recortada)",
        col.regions = brewer.pal(9, "GnBu"),
        legend = TRUE)



4.2 Temperatura

La Temperatura es un factor microclimático de alta relevancia para la fenología y sanidad del aguacate. El Modelo Exponencial (o el que haya resultado con el menor loss value) fue identificado como el más adecuado para representar la autocorrelación espacial de esta variable, gracias a su capacidad para capturar la rápida descorrelación a distancias cortas. El Kriging Ordinario se ejecuta utilizando los parámetros de este modelo para generar un Mapa de Predicción de Temperatura y su respectivo Mapa de Desviación Estándar. La cartografía resultante permitirá identificar microclimas de riesgo (zonas más calientes o más frías) que requieren manejo diferenciado.

# 1. Se crea los raster de salida
invisible(capture.output({
geodatos_ko = krige.conv(geo_temperatura, loc=geodatos_grid,
      krige = krige.control(nugget=0,trend.d="cte", 
      trend.l = "cte",cov.pars=c(sigmasq=3.0163, phi=0.000127 )))
}))
# 2. Temperatura_predict: Crea el ráster de las predicciones (Z = geodatos_ko$predict)
pred_matrix = cbind(geodatos_grid, geodatos_ko$predict)
Temperatura_predict = rast(pred_matrix, crs="EPSG:4326") 
names(Temperatura_predict) = "Temperatura_Predicha"

# 3. Temperatura_error: Crea el ráster del error (Z = sqrt(geodatos_ko$krige.var))
error_matrix = cbind(geodatos_grid, sqrt(geodatos_ko$krige.var))
Temperatura_error = rast(error_matrix, crs="EPSG:4326")
names(Temperatura_error) = "Error_StDev"

# 4. Visualizacion de mapas
# 4.1 Mapa 1: Predicción de Temperatura (Kriging)
plot_prediccion = ggplot() +
    # Usar geom_spatraster para dibujar el ráster de terra
    geom_spatraster(data = Temperatura_predict) + 
    # Usar una paleta de color moderna (viridis)
    scale_fill_viridis_c(name = "Temperatura Relativa (%)") +
    # Adicionar contornos para mejor lectura
    geom_spatraster_contour(data = Temperatura_predict, color = "gray30", alpha = 0.5) +
    labs(
        title = "Predicción Espacial (Kriging Ordinario)",
        x = "Longitud",
        y = "Latitud"
    ) +
    theme_minimal() +
    theme(legend.position = "bottom")

# 4.2 Mapa 2: Error de Predicción (Desviación Estándar)
plot_error = ggplot() +
    geom_spatraster(data = Temperatura_error) +
    # Usar otra paleta de color (Magma) para contraste. El error se visualiza en escalas más oscuras.
    scale_fill_viridis_c(option = "magma", name = "Desviación Estándar") + 
    geom_spatraster_contour(data = Temperatura_error, color = "black", alpha = 0.5) +
    labs(
        title = "Error de Predicción (Desviación Estándar)",
        x = "Longitud",
        y = "Latitud"
    ) +
    theme_minimal() +
    theme(legend.position = "bottom")

# 5. Combinar los dos mapas en un solo panel (Usando patchwork)
plot_final = plot_prediccion + plot_error + 
             plot_layout(guides = 'collect') + # Recolecta las leyendas
             plot_annotation(title = 'Resultados del Kriging Ordinario para Temperatura Relativa')

plot_final

# 1. Convertir los puntos filtrados a un objeto SF
puntos_sf = geo_datos_filtrado %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = "EPSG:4326")

# 2. Crear la Envolvente Convexa (Convex Hull)
# Esto crea el polígono más pequeño que contiene todos los puntos verdes.
area_estudio_polygon = puntos_sf %>%
  st_union() %>%
  st_convex_hull()

# Se valida del polígono para el recorte
area_estudio_terra = vect(area_estudio_polygon)

# 1. Recorte Inicial: Recorta el ráster al recuadro (bounding box) del polígono.
Temperatura_crop = crop(Temperatura_predict, area_estudio_terra)

# 2. Aplicar Máscara: Elimina los píxeles que están fuera del polígono (el área blanca).
Temperatura_recortada = mask(Temperatura_crop, area_estudio_terra)

# Visualizar el resultado con mapview (interactivo)
mapview(Temperatura_recortada,
        layer.name = "Temperatura Relativa (Recortada)",
        col.regions = brewer.pal(9, "YlOrRd"),
        legend = TRUE)



4.3 Altura

La Altura, siendo una variable topográfica y fundamentalmente estable, exhibió la estructura espacial más robusta y continua de todo el estudio, con un Efecto Pepita insignificante. Esto permite modelar su variación con la máxima confianza, utilizando el Modelo Exponencial para el Kriging Ordinario. El objetivo de este procedimiento no es solo generar un Mapa de Predicción de Altitud, sino establecer el Modelo Digital del Terreno (MDT) de la finca. Este MDT es la base para futuros análisis geomorfométricos y para interpretar los gradientes observados en las variables climáticas.

# 1. Se crea los raster de salida
invisible(capture.output({
geodatos_ko = krige.conv(geo_altura, loc=geodatos_grid,
      krige = krige.control(nugget=0,trend.d="cte", 
      trend.l = "cte",cov.pars=c(sigmasq=127.3645, phi=0.000639 )))
}))
# 2. Altura_predict: Crea el ráster de las predicciones (Z = geodatos_ko$predict)
pred_matrix = cbind(geodatos_grid, geodatos_ko$predict)
Altura_predict = rast(pred_matrix, crs="EPSG:4326") 
names(Altura_predict) = "Altura_Predicha"

# 3. Altura_error: Crea el ráster del error (Z = sqrt(geodatos_ko$krige.var))
error_matrix = cbind(geodatos_grid, sqrt(geodatos_ko$krige.var))
Altura_error = rast(error_matrix, crs="EPSG:4326")
names(Altura_error) = "Error_StDev"

# 4. Visualizacion de mapas
# 4.1 Mapa 1: Predicción de Altura (Kriging)
plot_prediccion = ggplot() +
    # Usar geom_spatraster para dibujar el ráster de terra
    geom_spatraster(data = Altura_predict) + 
    # Usar una paleta de color moderna (viridis)
    scale_fill_viridis_c(name = "Altura Relativa (%)") +
    # Adicionar contornos para mejor lectura
    geom_spatraster_contour(data = Altura_predict, color = "gray30", alpha = 0.5) +
    labs(
        title = "Predicción Espacial (Kriging Ordinario)",
        x = "Longitud",
        y = "Latitud"
    ) +
    theme_minimal() +
    theme(legend.position = "bottom")

# 4.2 Mapa 2: Error de Predicción (Desviación Estándar)
plot_error = ggplot() +
    geom_spatraster(data = Altura_error) +
    # Usar otra paleta de color (Magma) para contraste. El error se visualiza en escalas más oscuras.
    scale_fill_viridis_c(option = "magma", name = "Desviación Estándar") + 
    geom_spatraster_contour(data = Altura_error, color = "black", alpha = 0.5) +
    labs(
        title = "Error de Predicción (Desviación Estándar)",
        x = "Longitud",
        y = "Latitud"
    ) +
    theme_minimal() +
    theme(legend.position = "bottom")

# 5. Combinar los dos mapas en un solo panel (Usando patchwork)
plot_final = plot_prediccion + plot_error + 
             plot_layout(guides = 'collect') + # Recolecta las leyendas
             plot_annotation(title = 'Resultados del Kriging Ordinario para Altura Relativa')

plot_final

# 1. Convertir los puntos filtrados a un objeto SF
puntos_sf = geo_datos_filtrado %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = "EPSG:4326")

# 2. Crear la Envolvente Convexa (Convex Hull)
# Esto crea el polígono más pequeño que contiene todos los puntos verdes.
area_estudio_polygon = puntos_sf %>%
  st_union() %>%
  st_convex_hull()

# Se valida del polígono para el recorte
area_estudio_terra = vect(area_estudio_polygon)

# 1. Recorte Inicial: Recorta el ráster al recuadro (bounding box) del polígono.
Altura_crop = crop(Altura_predict, area_estudio_terra)

# 2. Aplicar Máscara: Elimina los píxeles que están fuera del polígono (el área blanca).
Altura_recortada = mask(Altura_crop, area_estudio_terra)

# Visualizar el resultado con mapview (interactivo)
mapview(Altura_recortada,
        layer.name = "Altura Relativa (Recortada)",
        col.regions = brewer.pal(9, "YlOrBr"),
        legend = TRUE)



5. Validación cruzada de los modelos por variable

Para asegurar la robustez y la precisión de los modelos predictivos generados, se implementa la técnica de Validación Cruzada (Cross-Validation). Dado que el Modelo Exponencial demostró el mejor ajuste global (o el menor valor de la función de pérdida) para la mayoría de las variables analizadas (Temperatura, Humedad Relativa y Altura), se utiliza este modelo para la validación. Este proceso esencial evalúa la capacidad del modelo para predecir valores en puntos no muestreados, comparando la estimación obtenida mediante Kriging con el valor real medido, lo cual confirma la validez de los parámetros antes de la interpretación final de los mapas.

5.1 Humedad relativa

# 1. Ejecutar la Validación Cruzada (Leave-One-Out)
# 'xvalid' realiza el Kriging en cada punto, excluyéndolo de la estimación.
valida_humedad = xvalid(
    geodata = geo_humedad_relativa,
    model = model_mco_exp_h
)
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
# 2. Resumen Estadístico de los Errores de Predicción
print(summary(valida_humedad))
                Min.    1st Qu.      Median        Mean   3rd Qu.      Max.
errors     -8.998917 -2.2136801 -0.10169526 0.005089921 1.9878688 17.861771
std.errors -2.472844 -0.6078266 -0.02799733 0.000672704 0.5490121  4.945636
                  sd
errors     3.2705069
std.errors 0.8987788
# 3. Gráfico de Diagnóstico (Predicho vs. Observado)
# Este gráfico es clave para la interpretación visual
plot(valida_humedad)



Conclusión de la Validación Cruzada (Modelo de Humedad Relativa)

El análisis de la validación cruzada para este modelo indica que, si bien cumple el requisito de insesgadez, su desempeño y precisión son deficientes, lo que invalida su uso para la predicción final.

  • Desempeño y Precisión

    • Bajo Sesgo: El modelo demuestra ser insesgado al registrar una media de los errores de ≈0.005. Esto indica que, en promedio, el modelo no subestima ni sobreestima los valores de manera sistemática.

    • Baja Precisión: La Desviación Estándar (SD) de los errores es de 3.27 unidades. Este valor es considerablemente alto. Por lo tanto, el error de predicción típico asociado a este modelo puede ser inaceptable para los objetivos de precisión requeridos por el estudio.

  • Incumplimiento de Supuestos

    • Correlación en Residuos (Subestimación Sistémica): Se observó una correlación y tendencia positiva en los residuos al graficarlos contra los valores reales. Esto indica que el modelo subestima sistematicamente los valores altos de la variable (es decir, cuando data>80), generando errores positivos significativos en estos rangos.

    • Heterocedasticidad: La varianza de los errores no es constante (heterocedasticidad), ya que la incertidumbre y la dispersión de los residuos aumentan significativamente en los extremos superiores de la variable predicha.

    • Normalidad: La distribución de los errores es cercana a la normal, aunque la alta dispersión y la asimetría en la cola derecha (por la subestimación de valores altos) son evidentes.

  • Recomendación

Será necesario aplicar una extensión del Kriging. Esto puede incluir la incorporación de una tendencia secundaria o el uso de CoKriging / Kriging con Tendencia Externa con una covariable relevante para modelar y corregir la tendencia lineal observada en los residuos.

5.2 Temperatura

# 1. Ejecutar la Validación Cruzada (Leave-One-Out)
# 'xvalid' realiza el Kriging en cada punto, excluyéndolo de la estimación.
valida_temperatura = xvalid(
    geodata = geo_temperatura,
    model = model_mco_exp_t
)
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
# 2. Resumen Estadístico de los Errores de Predicción
print(summary(valida_temperatura))
                Min.    1st Qu.       Median         Mean   3rd Qu.     Max.
errors     -3.674394 -0.6481378 -0.004007469 -0.009936681 0.6112504 2.784354
std.errors -3.418071 -0.5968083 -0.003748917 -0.004614072 0.5805013 2.801736
                  sd
errors     1.0224049
std.errors 0.9670636
# 3. Gráfico de Diagnóstico (Predicho vs. Observado)
# Este gráfico es clave para la interpretación visual
plot(valida_temperatura)



Conclusión de la Validación Cruzada (Modelo de Temperatura)

El análisis de la validación cruzada indica que este modelo de predicción ha alcanzado un desempeño predictivo excelente, cumpliendo de manera satisfactoria con los supuestos estadísticos fundamentales para la aplicación del Kriging.

  • Desempeño y Robustez del Modelo

    • Alta Precisión: La Desviación Estándar (SD) de los errores es de 1.02 unidades. Este valor bajo se considera excelente, ya que representa la magnitud del error de predicción típico del modelo.

    • Insessgado: El modelo es insesgado, con una media de los errores de ≈−0.01. Esto confirma que el modelo está bien calibrado y no presenta una subestimación o sobreestimación sistémica de la variable.

  • Cumplimiento de Supuestos Estadísticos

    • Homocedasticidad: Se verifica el cumplimiento del supuesto de varianza constante (homocedasticidad). La dispersión de los errores es uniforme a lo largo de todo el rango de valores predichos, lo que asegura que la precisión de la predicción es consistente independientemente del valor de la variable.

    • Normalidad: La distribución de los errores se ajusta de manera muy cercana a una distribución normal (confirmado por el Gráfico Q-Q y el Histograma). Este cumplimiento es crucial, ya que garantiza la confiabilidad de los intervalos de confianza generados durante la predicción.

  • Recomendación Final

El modelo es estadísticamente robusto y preciso para la variable de temperatura. No se identificaron debilidades significativas en el análisis de residuos.

5.3 Altura

# 1. Ejecutar la Validación Cruzada (Leave-One-Out)
# 'xvalid' realiza el Kriging en cada punto, excluyéndolo de la estimación.
valida_altura = xvalid(
    geodata = geo_altura,
    model = model_mco_exp_a
)
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
# 2. Resumen Estadístico de los Errores de Predicción
print(summary(valida_altura))
                 Min.    1st Qu.      Median         Mean   3rd Qu.     Max.
errors     -16.919682 -3.9293200 -0.22394466 0.0039189854 2.7940335 16.19836
std.errors  -2.724152 -0.6326585 -0.03613168 0.0003143015 0.4464078  2.61632
                  sd
errors     5.9409564
std.errors 0.9610738
# 3. Gráfico de Diagnóstico (Predicho vs. Observado)
# Este gráfico es clave para la interpretación visual
plot(valida_altura)



Conclusión de la Validación Cruzada (Modelo de Altura)

El análisis de la validación cruzada indica que el modelo de predicción es insesgado, pero sufre de problemas significativos de precisión, heterocedasticidad y normalidad que comprometen la validez estadística de los resultados del Kriging.

  • Desempeño General del Modelo

    • Bajo Sesgo (Punto Fuerte): El modelo cumple con el requisito de insesgadez, ya que la media de los errores es prácticamente cero (≈0.004) y su distribución está centrada en el histograma. Esto valida que no existe una subestimación o sobreestimación sistemática.

    • Precisión Insuficiente (Punto Débil): La Desviación Estándar (SD) de los errores es de 5.94 unidades, un valor considerado alto. Esta magnitud representa el error de predicción típico del modelo, y su aceptabilidad debe ser evaluada en el contexto específico de la aplicación.

  • Incumplimiento de Supuestos Estadísticos

    • Heterocedasticidad (Problema Crítico): Se observó un patrón de “abanico” en la gráfica de residuos versus valores predichos. Este fenómeno de varianza no constante indica que la precisión del modelo varía a lo largo del rango de la variable, mostrando una mayor incertidumbre y error al predecir los valores bajos.

    • No-Normalidad: El gráfico Q-Q y el histograma de errores revelaron que la distribución de los residuos no sigue una curva normal perfecta, especialmente en las colas (valores extremos). La violación de este supuesto compromete la confiabilidad de los intervalos de confianza calculados a partir del Kriging.

  • Recomendación Inmediata Para corregir los problemas de heterocedasticidad y no-normalidad, se recomienda aplicar una transformación matemática a la variable de interés, como la transformación Box-Cox o logarítmica. Este procedimiento tiene por objetivo estabilizar la varianza (lograr homocedasticidad) y normalizar la distribución de los datos, mejorando la robustez y la precisión del modelo geoestadístico.