1 Introducción Actividad

El presente documento desarrolla una exploración de datos espaciales mediante herramientas de geoestadística. La actividad se enfoca en la variable Temperature, registrada en una finca de aguacate, con el propósito de analizar su comportamiento espacial e interpolar la temperatura en zonas no muestreadas mediante kriging.

De acuerdo con las indicaciones del caso, el análisis se restringe únicamente al primer periodo de medición del 01/10/2020, identificado en el campo FORMATTED_DATE_TIME. Para este periodo se espera trabajar con 534 registros, correspondientes a puntos espaciales con coordenadas de latitud y longitud.

2 Carga de paquetes

options(repos = c(CRAN = "https://cloud.r-project.org"))

paquetes <- c(
  "readxl", "dplyr", "tidyr", "stringr", "lubridate", "janitor",
  "ggplot2", "knitr", "kableExtra", "purrr", "tibble",
  "sf", "sp", "gstat", "spdep", "viridis", "patchwork"
)

paquetes_faltantes <- paquetes[!vapply(paquetes, requireNamespace, logical(1), quietly = TRUE)]

if (length(paquetes_faltantes) > 0) {
  install.packages(paquetes_faltantes, dependencies = TRUE)
}

invisible(lapply(paquetes, library, character.only = TRUE))

3 Carga de datos

Se toma como fuente de datos el archivo Datos_Completos_Aguacate.xlsx relacionado en la actividad.

archivo_excel <- "Datos_Completos_Aguacate.xlsx"

if (!file.exists(archivo_excel)) {
  stop("No se encontró el archivo Datos_Completos_Aguacate.xlsx.")
}

datos_raw <- readxl::read_excel(archivo_excel, sheet = 1) |>
  janitor::clean_names()

dim(datos_raw)
## [1] 20271    21
tibble(
  numero = seq_along(names(datos_raw)),
  variable = names(datos_raw)
) |>
  knitr::kable(caption = "Variables disponibles en la base de datos") |>
  kableExtra::kable_styling(full_width = FALSE)
Variables disponibles en la base de datos
numero variable
1 id_arbol
2 latitude
3 longitude
4 formatted_date_time
5 psychro_wet_bulb_temperature
6 station_pressure
7 relative_humidity
8 crosswind
9 temperature
10 barometric_pressure
11 headwind
12 direction_true
13 direction_mag
14 wind_speed
15 heat_stress_index
16 altitude
17 dew_point
18 density_altitude
19 wind_chill
20 estado_fenologico_predominante
21 frutos_afectados

4 Limpieza de fecha y selección del primer periodo

La columna original FORMATTED_DATE_TIME contiene fecha y hora en formato de texto. Para cumplir con la indicación del caso, se extrae únicamente la fecha y se filtra el periodo 01/10/2020.

extraer_fecha <- function(x) {
  x <- as.character(x)
  x <- stringr::str_replace_all(x, "\\u00A0", " ")
  x <- stringr::str_squish(x)
  fecha_txt <- stringr::str_extract(x, "\\d{1,2}/\\d{1,2}/\\d{4}")
  lubridate::dmy(fecha_txt)
}

datos <- datos_raw |>
  mutate(
    fecha_medicion = extraer_fecha(formatted_date_time),
    latitude = as.numeric(latitude),
    longitude = as.numeric(longitude),
    temperature = as.numeric(temperature)
  )

fecha_objetivo <- lubridate::dmy("01/10/2020")

datos_periodo <- datos |>
  filter(fecha_medicion == fecha_objetivo) |>
  filter(!is.na(latitude), !is.na(longitude), !is.na(temperature))

validacion_periodo <- tibble(
  fecha_analizada = as.character(fecha_objetivo),
  registros_periodo = nrow(datos_periodo),
  registros_esperados = 534,
  cumple_cantidad_esperada = nrow(datos_periodo) == 534,
  coordenadas_unicas = dplyr::n_distinct(
  paste(datos_periodo$latitude, datos_periodo$longitude)),
  temperaturas_faltantes = sum(is.na(datos_periodo$temperature))
)

validacion_periodo |>
  knitr::kable(caption = "Validación del periodo seleccionado") |>
  kableExtra::kable_styling(full_width = FALSE)
Validación del periodo seleccionado
fecha_analizada registros_periodo registros_esperados cumple_cantidad_esperada coordenadas_unicas temperaturas_faltantes
2020-10-01 534 534 TRUE 534 0 
if (nrow(datos_periodo) != 534) {
  stop("La selección del periodo 01/10/2020 no generó 534 registros. Revise el formato de la columna FORMATTED_DATE_TIME.")
}

El análisis se realizará exclusivamente con los 534 registros del periodo 01/10/2020. Esto garantiza que no se mezclen diferentes momentos de medición y que la interpolación espacial represente un único instante de observación.

5 Análisis exploratorio de datos

5.1 Resumen descriptivo de Temperature

resumen_temperature <- datos_periodo |>
  summarise(
    n = n(),
    minimo = min(temperature, na.rm = TRUE),
    q1 = quantile(temperature, 0.25, na.rm = TRUE),
    media = mean(temperature, na.rm = TRUE),
    mediana = median(temperature, na.rm = TRUE),
    q3 = quantile(temperature, 0.75, na.rm = TRUE),
    maximo = max(temperature, na.rm = TRUE),
    desviacion_estandar = sd(temperature, na.rm = TRUE),
    coeficiente_variacion = sd(temperature, na.rm = TRUE) / mean(temperature, na.rm = TRUE)
  )

resumen_temperature |>
  knitr::kable(caption = "Resumen descriptivo de la variable Temperature") |>
  kableExtra::kable_styling(full_width = FALSE)
Resumen descriptivo de la variable Temperature
n minimo q1 media mediana q3 maximo desviacion_estandar coeficiente_variacion
534 22.2 24.5 25.82903 25.8 27.175 29.7 1.771073 0.0685691

5.2 Distribución de la temperatura

ggplot(datos_periodo, aes(x = temperature)) +
  geom_histogram(bins = 25, color = "white") +
  labs(
    title = "Distribución de la variable Temperature",
    subtitle = "Periodo analizado: 01/10/2020",
    x = "Temperature",
    y = "Frecuencia"
  ) +
  theme_minimal()

La distribución de la variable Temperature muestra valores aproximados entre 22 °C y 30 °C, con mayor concentración de observaciones entre 24 °C y 28 °C. El histograma evidencia una concentración importante alrededor de los 26 °C, lo cual sugiere que la temperatura presenta una distribución relativamente concentrada, aunque con variabilidad suficiente para justificar un análisis espacial.

ggplot(datos_periodo, aes(y = temperature)) +
  geom_boxplot() +
  labs(
    title = "Diagrama de caja de Temperature",
    subtitle = "Evaluación inicial de dispersión y posibles valores atípicos",
    y = "Temperature"
  ) +
  theme_minimal()

El diagrama de caja confirma que la mediana de la temperatura se ubica aproximadamente en 25.8 °C. La mayor parte de los datos se encuentra entre aproximadamente 24.5 °C y 27.2 °C. No se observan valores atípicos extremos marcados en el boxplot, lo que indica que los datos son consistentes para continuar con el análisis geoestadístico sin necesidad de eliminar observaciones por comportamiento extremo evidente.

5.3 Exploración espacial inicial en coordenadas geográficas

ggplot(datos_periodo, aes(x = longitude, y = latitude, color = temperature)) +
  geom_point(size = 2, alpha = 0.9) +
  scale_color_viridis_c(option = "C") +
  coord_equal() +
  labs(
    title = "Distribución espacial de Temperature",
    subtitle = "Puntos observados en coordenadas geográficas",
    x = "Longitud",
    y = "Latitud",
    color = "Temperature"
  ) +
  theme_minimal()

El mapa de puntos en coordenadas geográficas permite observar que la temperatura no se distribuye de manera completamente aleatoria en el espacio. Existen zonas con temperaturas más altas, representadas por colores amarillos y naranjas, y zonas con temperaturas más bajas, representadas por tonos morados.

Esta primera visualización sugiere la existencia de una posible estructura espacial, ya que puntos cercanos tienden a presentar valores similares de temperatura.

6 Creación del geodata

Para el análisis geoestadístico se crea un objeto espacial utilizando las coordenadas de latitud y longitud con la librería sf. Inicialmente se declaran las coordenadas en el sistema geográfico WGS84, EPSG:4326. Posteriormente se transforman a un sistema proyectado en metros, UTM zona 18N, EPSG:32618, apropiado para calcular distancias espaciales.

datos_sf <- datos_periodo |>
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326, remove = FALSE)

datos_utm <- datos_sf |>
  st_transform(32618)

coords_utm <- st_coordinates(datos_utm)

x_centro <- mean(coords_utm[, 1])
y_centro <- mean(coords_utm[, 2])

datos_utm <- datos_utm |>
  mutate(
    x = coords_utm[, 1],
    y = coords_utm[, 2],
    x_c = x - x_centro,
    y_c = y - y_centro
  )

datos_sp <- as(datos_utm, "Spatial")

class(datos_sf)
## [1] "sf"         "tbl_df"     "tbl"        "data.frame"
class(datos_utm)
## [1] "sf"         "tbl_df"     "tbl"        "data.frame"
class(datos_sp)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
datos_utm |>
  st_drop_geometry() |>
  summarise(
    min_x = min(x),
    max_x = max(x),
    rango_x_metros = max(x) - min(x),
    min_y = min(y),
    max_y = max(y),
    rango_y_metros = max(y) - min(y)
  ) |>
  knitr::kable(caption = "Resumen de coordenadas proyectadas en UTM") |>
  kableExtra::kable_styling(full_width = FALSE)
Resumen de coordenadas proyectadas en UTM
min_x max_x rango_x_metros min_y max_y rango_y_metros
309656.1 309832.4 176.3104 264519.2 264688.6 169.4073

El resumen de coordenadas proyectadas muestra que el área de estudio tiene una extensión aproximada de:

176.31 metros en el eje X. 169.41 metros en el eje Y.

Esto indica que el análisis se realiza sobre una zona relativamente pequeña y bien delimitada, lo cual es apropiado para evaluar variabilidad espacial local dentro del cultivo o zona de monitoreo.

7 Evaluación de tendencia espacial

La geoestadística clásica asume estacionariedad de segundo orden o, en términos prácticos, que la media no cambia sistemáticamente en el espacio. Por ello, antes de calcular el semivariograma y aplicar kriging, se evalúa si la temperatura presenta tendencia respecto a las coordenadas espaciales.

7.1 Tendencia visual respecto a coordenadas X y Y

La evaluación de tendencia se realizó analizando la relación entre Temperature y las coordenadas centradas x_c y y_c.

ggplot(st_drop_geometry(datos_utm), aes(x = x_c, y = temperature)) +
  geom_point(alpha = 0.7) +
  geom_smooth(method = "lm", se = TRUE) +
  labs(
    title = "Evaluación de tendencia de Temperature respecto al eje X",
    x = "Coordenada X centrada, metros",
    y = "Temperature"
  ) +
  theme_minimal()

ggplot(st_drop_geometry(datos_utm), aes(x = y_c, y = temperature)) +
  geom_point(alpha = 0.7) +
  geom_smooth(method = "lm", se = TRUE) +
  labs(
    title = "Evaluación de tendencia de Temperature respecto al eje Y",
    x = "Coordenada Y centrada, metros",
    y = "Temperature"
  ) +
  theme_minimal()

Los gráficos de tendencia muestran una relación negativa entre la temperatura y ambos ejes espaciales. Esto significa que, a medida que aumentan las coordenadas X o Y, la temperatura tiende a disminuir ligeramente.

7.2 Modelo de tendencia lineal

modelo_tendencia <- lm(temperature ~ x_c + y_c, data = st_drop_geometry(datos_utm))
resumen_tendencia <- summary(modelo_tendencia)
resumen_tendencia
## 
## Call:
## lm(formula = temperature ~ x_c + y_c, data = st_drop_geometry(datos_utm))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0761 -1.1765 -0.0286  1.1375  4.1201 
## 
## Coefficients:
##              Estimate Std. Error t value             Pr(>|t|)    
## (Intercept) 25.829026   0.073548 351.187 < 0.0000000000000002 ***
## x_c         -0.007194   0.002194  -3.279             0.001108 ** 
## y_c         -0.006433   0.001928  -3.337             0.000905 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.7 on 531 degrees of freedom
## Multiple R-squared:  0.08257,    Adjusted R-squared:  0.07911 
## F-statistic: 23.89 on 2 and 531 DF,  p-value: 0.0000000001157
coef_tendencia <- as.data.frame(coef(resumen_tendencia)) |>
  tibble::rownames_to_column("termino") |>
  rename(
    estimacion = Estimate,
    error_estandar = `Std. Error`,
    valor_t = `t value`,
    valor_p = `Pr(>|t|)`
  )

coef_tendencia |>
  knitr::kable(caption = "Coeficientes del modelo de tendencia espacial") |>
  kableExtra::kable_styling(full_width = FALSE)
Coeficientes del modelo de tendencia espacial
termino estimacion error_estandar valor_t valor_p
(Intercept) 25.8290262 0.0735477 351.187242 0.0000000
x_c -0.0071940 0.0021937 -3.279460 0.0011082
y_c -0.0064329 0.0019275 -3.337378 0.0009050 
p_tendencia <- coef_tendencia |>
  filter(termino %in% c("x_c", "y_c")) |>
  pull(valor_p)

hay_tendencia <- any(p_tendencia < 0.05, na.rm = TRUE)

metodo_kriging <- ifelse(
  hay_tendencia,
  "Kriging universal, porque se identificó tendencia espacial significativa en función de las coordenadas.",
  "Kriging ordinario, porque no se identificó tendencia espacial significativa en función de las coordenadas."
)

metodo_kriging
## [1] "Kriging universal, porque se identificó tendencia espacial significativa en función de las coordenadas."

El modelo lineal de tendencia presentó los siguientes resultados principales:

  • Para x_c, el coeficiente fue negativo y estadísticamente significativo, con un valor p aproximado de 0.0011.
  • Para y_c, el coeficiente también fue negativo y estadísticamente significativo, con un valor p aproximado de 0.0009.

El modelo tuvo un R² aproximado de 0.0826, es decir, las coordenadas explican cerca del 8.26 % de la variabilidad de la temperatura.

Aunque el porcentaje explicado no es alto, los coeficientes son estadísticamente significativos. Por tanto, se concluye que existe una tendencia espacial significativa en la variable Temperature.

Debido a esta tendencia, la metodología más adecuada para la interpolación no es el kriging ordinario, sino el kriging universal, ya que este permite incorporar la tendencia espacial asociada a las coordenadas.

La decisión metodológica para la interpolación será: Kriging universal, porque se identificó tendencia espacial significativa en función de las coordenadas.

8 Autocorrelación espacial

La autocorrelación espacial permite evaluar si valores cercanos en el espacio tienden a ser más similares entre sí que valores lejanos. Para ello se utiliza el estadístico Moran’s I, construyendo una matriz de vecinos con los 6 puntos más cercanos.

coords <- cbind(datos_utm$x, datos_utm$y)

vecinos_knn <- spdep::knearneigh(coords, k = 6)
vecinos_nb <- spdep::knn2nb(vecinos_knn)
vecinos_w <- spdep::nb2listw(vecinos_nb, style = "W")

moran_temperature <- spdep::moran.test(datos_utm$temperature, vecinos_w)
moran_temperature
## 
##  Moran I test under randomisation
## 
## data:  datos_utm$temperature  
## weights: vecinos_w    
## 
## Moran I statistic standard deviate = 24.824, p-value <
## 0.00000000000000022
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.5997716415     -0.0018761726      0.0005874067
spdep::moran.plot(
  datos_utm$temperature,
  vecinos_w,
  labels = FALSE,
  main = "Diagrama de Moran para Temperature"
)

El estadístico de Moran’s I obtenido fue aproximadamente:

  • Moran’s I = 0.5998

con un valor p menor a 0.05, prácticamente cercano a cero.

Este resultado indica una autocorrelación espacial positiva fuerte en la variable Temperature. En otras palabras, los puntos cercanos espacialmente tienden a tener temperaturas similares. Esto confirma que la temperatura no se comporta como una variable espacialmente aleatoria, sino que presenta dependencia espacial.

9 Semivariograma empírico

El semivariograma mide cómo cambia la semivarianza entre observaciones a medida que aumenta la distancia espacial. Es una herramienta central para decidir si existe estructura espacial y para parametrizar el kriging.

distancias <- as.matrix(dist(coords))
distancia_maxima <- max(distancias)
cutoff_variograma <- distancia_maxima * 0.5
width_variograma <- cutoff_variograma / 15

formula_geo <- if (hay_tendencia) {
  temperature ~ x_c + y_c
} else {
  temperature ~ 1
}

formula_geo
## temperature ~ x_c + y_c
vg_empirico <- gstat::variogram(
  formula_geo,
  data = datos_sp,
  cutoff = cutoff_variograma,
  width = width_variograma
)

vg_empirico |>
  head(10) |>
  knitr::kable(caption = "Primeros valores del semivariograma empírico") |>
  kableExtra::kable_styling(full_width = FALSE)
Primeros valores del semivariograma empírico
np dist gamma dir.hor dir.ver id
1086 5.750208 1.131413 0 0 var1
3339 11.310900 1.919209 0 0 var1
5296 18.298557 2.484671 0 0 var1
6983 25.473851 2.744942 0 0 var1
8075 32.712558 2.802815 0 0 var1
8893 39.889848 2.735930 0 0 var1
9576 47.091770 2.783918 0 0 var1
10017 54.343093 2.958415 0 0 var1
9946 61.593082 3.131483 0 0 var1
9809 68.805459 3.065673 0 0 var1 
ggplot(vg_empirico, aes(x = dist, y = gamma)) +
  geom_point(aes(size = np), alpha = 0.8) +
  labs(
    title = "Semivariograma empírico de Temperature",
    subtitle = "El tamaño del punto representa el número de pares por intervalo de distancia",
    x = "Distancia, metros",
    y = "Semivarianza",
    size = "Pares"
  ) +
  theme_minimal()

El semivariograma empírico muestra que la semivarianza aumenta conforme aumenta la distancia entre los puntos. Esto significa que las observaciones cercanas tienen valores de temperatura más parecidos, mientras que las observaciones más alejadas presentan diferencias mayores.

En los primeros intervalos de distancia, la semivarianza es menor. Por ejemplo, a distancias cercanas a 5.75 metros, la semivarianza es aproximadamente 1.13. Luego aumenta progresivamente hasta ubicarse alrededor de valores cercanos a 3.0 en distancias mayores.

Este comportamiento es coherente con la teoría geoestadística: cuando existe autocorrelación espacial, los puntos cercanos son más similares y, por tanto, presentan menor semivarianza. A medida que la distancia aumenta, la similitud disminuye y la semivarianza crece hasta aproximarse a una meseta.

10 Semivariograma direccional

El semivariograma direccional permite evaluar si el comportamiento espacial de la temperatura cambia según la dirección analizada. En la gráfica se comparan las direcciones 0°, 45°, 90° y 135°.

vg_direccional <- gstat::variogram(
  formula_geo,
  data = datos_sp,
  alpha = c(0, 45, 90, 135),
  cutoff = cutoff_variograma,
  width = width_variograma
)

head(vg_direccional)
ggplot(vg_direccional, aes(x = dist, y = gamma, color = factor(dir.hor))) +
  geom_point(alpha = 0.7) +
  geom_line() +
  labs(
    title = "Semivariograma direccional de Temperature",
    subtitle = "Evaluación exploratoria de anisotropía",
    x = "Distancia, metros",
    y = "Semivarianza",
    color = "Dirección"
  ) +
  theme_minimal()

Se observa que las curvas no son completamente iguales entre direcciones. Algunas direcciones presentan mayor semivarianza que otras, especialmente la dirección de 135°, que muestra mayor variabilidad en algunos intervalos de distancia. Esto sugiere una posible anisotropía, es decir, que la dependencia espacial podría variar según la dirección.

Sin embargo, para efectos del análisis principal, se trabajó con un modelo isotrópico general.

11 Ajuste de modelos teóricos de semivariograma

Se compararon tres modelos teóricos de semivariograma:

  • Esférico: crecimiento gradual hasta alcanzar una meseta.
  • Exponencial: crecimiento rápido inicial y acercamiento progresivo a la meseta.
  • Gaussiano: comportamiento más suave cerca del origen.

El mejor modelo se selecciona con base en el menor error de ajuste reportado por fit.variogram, usando el atributo SSErr.

var_temperature <- var(datos_utm$temperature, na.rm = TRUE)
rango_inicial <- max(vg_empirico$dist, na.rm = TRUE) / 3

ajustar_modelo <- function(modelo) {
  modelo_inicial <- gstat::vgm(
    psill = var_temperature * 0.8,
    model = modelo,
    range = rango_inicial,
    nugget = var_temperature * 0.2
  )
  
  ajuste <- try(
    suppressWarnings(gstat::fit.variogram(vg_empirico, modelo_inicial)),
    silent = TRUE
  )
  
  if (inherits(ajuste, "try-error")) {
    return(tibble(
      modelo = modelo,
      nugget = NA_real_,
      psill = NA_real_,
      rango = NA_real_,
      sse = Inf,
      ajuste = list(NULL)
    ))
  }
  
  tibble(
    modelo = modelo,
    nugget = ajuste$psill[1],
    psill = ajuste$psill[2],
    rango = ajuste$range[2],
    sse = attr(ajuste, "SSErr"),
    ajuste = list(ajuste)
  )
}

modelos_evaluados <- c("Sph", "Exp", "Gau")
ajustes_modelos <- purrr::map_dfr(modelos_evaluados, ajustar_modelo)

ajustes_resumen <- ajustes_modelos |>
  mutate(
    modelo_nombre = case_when(
      modelo == "Sph" ~ "Esférico",
      modelo == "Exp" ~ "Exponencial",
      modelo == "Gau" ~ "Gaussiano",
      TRUE ~ modelo
    )
  ) |>
  select(modelo_nombre, modelo, nugget, psill, rango, sse) |>
  arrange(sse)

ajustes_resumen |>
  knitr::kable(caption = "Comparación de modelos teóricos de semivariograma") |>
  kableExtra::kable_styling(full_width = FALSE)
Comparación de modelos teóricos de semivariograma
modelo_nombre modelo nugget psill rango sse
Exponencial Exp 0.0000000 3.006653 11.37984 0.7687701
Esférico Sph 0.3473595 2.509836 26.42126 1.0470685
Gaussiano Gau 0.8478602 1.997911 13.83694 1.0985083 
ajustes_validos <- ajustes_modelos |>
  filter(is.finite(sse), !is.na(psill), psill > 0, !is.na(rango), rango > 0) |>
  arrange(sse)

if (nrow(ajustes_validos) == 0) {
  stop("No fue posible ajustar un modelo teórico de semivariograma válido. Revise la variable, las coordenadas o los parámetros iniciales.")
}

mejor_modelo <- ajustes_validos |> slice(1)
modelo_variograma_final <- mejor_modelo$ajuste[[1]]
nombre_mejor_modelo <- case_when(
  mejor_modelo$modelo == "Sph" ~ "Esférico",
  mejor_modelo$modelo == "Exp" ~ "Exponencial",
  mejor_modelo$modelo == "Gau" ~ "Gaussiano",
  TRUE ~ mejor_modelo$modelo
)

modelo_variograma_final
nombre_mejor_modelo
## [1] "Exponencial"

Los resultados fueron:

  • Exponencial: SSE aproximado de 0.7688.
  • Esférico: SSE aproximado de 1.0471.
  • Gaussiano: SSE aproximado de 1.0985.

El modelo con menor error fue el modelo exponencial, por lo cual se seleccionó como el modelo teórico más adecuado para representar la estructura espacial de la temperatura.

lineas_modelos <- ajustes_validos |>
  mutate(modelo_nombre = case_when(
    modelo == "Sph" ~ "Esférico",
    modelo == "Exp" ~ "Exponencial",
    modelo == "Gau" ~ "Gaussiano",
    TRUE ~ modelo
  )) |>
  mutate(linea = purrr::map(ajuste, ~gstat::variogramLine(.x, maxdist = max(vg_empirico$dist), n = 100))) |>
  select(modelo_nombre, linea) |>
  tidyr::unnest(linea)

ggplot() +
  geom_point(data = vg_empirico, aes(x = dist, y = gamma, size = np), alpha = 0.7) +
  geom_line(data = lineas_modelos, aes(x = dist, y = gamma, color = modelo_nombre), size = 1) +
  labs(
    title = "Ajuste de modelos teóricos al semivariograma empírico",
    subtitle = paste("Modelo seleccionado:", nombre_mejor_modelo),
    x = "Distancia, metros",
    y = "Semivarianza",
    size = "Pares",
    color = "Modelo"
  ) +
  theme_minimal()

El modelo exponencial ajustado tuvo los siguientes parámetros principales:

  • Nugget: 0.
  • Sill parcial: aproximadamente 3.0067.
  • Rango: aproximadamente 11.38 metros.

En términos prácticos, el modelo indica que existe dependencia espacial a distancias cortas y que la variabilidad de la temperatura se estabiliza conforme aumenta la distancia entre puntos.

12 Predicción espacial mediante kriging

La predicción espacial se realiza sobre una malla regular de 4.921 puntos construida dentro del polígono convexo que cubre los puntos observados. La resolución de la malla se define en metros, aprovechando que las coordenadas fueron proyectadas a UTM.

poligono_finca <- datos_utm |>
  st_union() |>
  st_convex_hull()

grilla_sf <- st_sf(
  geometry = st_make_grid(
    poligono_finca,
    cellsize = 2,
    what = "centers"
  )
)

st_crs(grilla_sf) <- st_crs(datos_utm)

grilla_sf <- grilla_sf |>
  st_filter(poligono_finca, .predicate = st_intersects)

coords_grilla <- st_coordinates(grilla_sf)

grilla_sf <- grilla_sf |>
  mutate(
    x = coords_grilla[, 1],
    y = coords_grilla[, 2],
    x_c = x - x_centro,
    y_c = y - y_centro
  )

grilla_sp <- as(grilla_sf, "Spatial")

nrow(grilla_sf)
## [1] 4921
resultado_kriging <- gstat::krige(
  formula = formula_geo,
  locations = datos_sp,
  newdata = grilla_sp,
  model = modelo_variograma_final
)
## [using universal kriging]
kriging_df <- as.data.frame(resultado_kriging)

head(kriging_df)

13 Mapa de predicción espacial de Temperature

ggplot() +
  geom_raster(
    data = kriging_df,
    aes(x = coords.x1, y = coords.x2, fill = var1.pred)
  ) +
  geom_point(
    data = st_drop_geometry(datos_utm),
    aes(x = x, y = y),
    size = 0.8,
    alpha = 0.6
  ) +
  scale_fill_viridis_c(option = "C") +
  coord_equal() +
  labs(
    title = "Predicción espacial de Temperature mediante kriging",
    subtitle = paste("Periodo: 01/10/2020 | Modelo de semivariograma:", nombre_mejor_modelo),
    x = "Coordenada X UTM, metros",
    y = "Coordenada Y UTM, metros",
    fill = "Temperature\npredicha"
  ) +
  theme_minimal()

El mapa muestra la distribución espacial estimada de la temperatura para las zonas donde no se tenían mediciones directas.

En el mapa se observan zonas con temperaturas más altas, representadas en tonos amarillos y naranjas, y zonas con temperaturas más bajas, representadas en tonos morados. También se evidencian franjas o patrones espaciales que muestran variaciones locales de temperatura dentro del área analizada.

El resultado permite visualizar cómo se comporta la temperatura en toda la superficie del área de estudio, no solo en los puntos donde se tomaron mediciones.

14 Mapa de incertidumbre de la predicción

La varianza de kriging representa la incertidumbre asociada a la predicción. En general, la incertidumbre tiende a ser menor cerca de puntos observados y mayor en zonas con menor densidad de muestreo.

ggplot() +
  geom_raster(
    data = kriging_df,
    aes(x = coords.x1, y = coords.x2, fill = var1.var)
  ) +
  geom_point(
    data = st_drop_geometry(datos_utm),
    aes(x = x, y = y),
    size = 0.8,
    alpha = 0.6
  ) +
  scale_fill_viridis_c(option = "B") +
  coord_equal() +
  labs(
    title = "Varianza de kriging para Temperature",
    subtitle = "Mapa de incertidumbre de la interpolación espacial",
    x = "Coordenada X UTM, metros",
    y = "Coordenada Y UTM, metros",
    fill = "Varianza\nde kriging"
  ) +
  theme_minimal()

Las zonas con menor varianza corresponden a áreas donde el modelo tiene mayor confianza en la predicción, generalmente porque están cerca de puntos observados o en zonas con alta densidad de muestreo.

Las zonas con mayor varianza se ubican principalmente hacia los bordes del área de estudio, especialmente en sectores donde hay menor soporte de puntos cercanos. Esto es esperable, ya que el kriging suele ser menos preciso en los límites del área analizada y en zonas con menor cantidad de observaciones cercanas.

Este mapa es importante porque permite diferenciar entre las zonas donde la predicción es más confiable y aquellas donde debe interpretarse con mayor precaución.

15 Validación cruzada del kriging

La validación cruzada permite evaluar el desempeño predictivo del modelo geoestadístico. Se retira cada observación, se predice su valor con las restantes y se comparan los valores observados contra los predichos.

cv_kriging <- gstat::krige.cv(
  formula = formula_geo,
  locations = datos_sp,
  model = modelo_variograma_final,
  nfold = nrow(datos_sp)
)

cv_df <- as.data.frame(cv_kriging)

metricas_cv <- cv_df |>
  summarise(
    me = mean(residual, na.rm = TRUE),
    mae = mean(abs(residual), na.rm = TRUE),
    rmse = sqrt(mean(residual^2, na.rm = TRUE)),
    correlacion_obs_pred = cor(observed, var1.pred, use = "complete.obs")
  )

metricas_cv |>
  knitr::kable(caption = "Métricas de validación cruzada del kriging") |>
  kableExtra::kable_styling(full_width = FALSE)
Métricas de validación cruzada del kriging
me mae rmse correlacion_obs_pred
-0.0108918 0.786249 1.014373 0.8196698 
ggplot(cv_df, aes(x = observed, y = var1.pred)) +
  geom_point(alpha = 0.7) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
  labs(
    title = "Validación cruzada: valores observados vs. predichos",
    x = "Temperature observada",
    y = "Temperature predicha"
  ) +
  theme_minimal()

La validación cruzada permitió evaluar el desempeño predictivo del modelo geoestadístico. Las métricas obtenidas fueron:

  • ME: aproximadamente -0.0109.
  • MAE: aproximadamente 0.7862.
  • RMSE: aproximadamente 1.0144.
  • Correlación entre observado y predicho: aproximadamente 0.8197.

El valor de ME está muy cerca de cero, lo cual indica que el modelo no presenta un sesgo sistemático importante. Es decir, en promedio, el modelo no sobreestima ni subestima fuertemente la temperatura.

El MAE indica que el error absoluto promedio de predicción es de aproximadamente 0.79 °C, mientras que el RMSE indica un error típico cercano a 1.01 °C. Estos valores son razonables para una variable ambiental como la temperatura, especialmente considerando la variabilidad espacial observada.

La correlación de aproximadamente 0.82 entre valores observados y predichos indica una buena capacidad predictiva del modelo. El gráfico de observados contra predichos muestra que muchos puntos se ubican cerca de la línea diagonal, lo que confirma que las predicciones guardan una relación fuerte con los valores reales.

ggplot(cv_df, aes(x = residual)) +
  geom_histogram(bins = 25, color = "white") +
  labs(
    title = "Distribución de residuos de validación cruzada",
    x = "Residuo: observado - predicho",
    y = "Frecuencia"
  ) +
  theme_minimal()

El histograma de residuos muestra que los errores se concentran alrededor de cero, lo cual refuerza la idea de que el modelo tiene un comportamiento balanceado y no presenta un sesgo marcado.

16 Interpretación de resultados

  1. El análisis geoestadístico realizado sobre la variable Temperature para el periodo 01/10/2020 permitió identificar una estructura espacial clara en los datos. La variable presentó autocorrelación espacial positiva, tendencia significativa respecto a las coordenadas y un patrón de dependencia espacial evidenciado en el semivariograma empírico.

  2. El modelo teórico que mejor representó la estructura espacial fue el modelo exponencial, al presentar el menor error de ajuste frente a los modelos esférico y gaussiano. Debido a la existencia de tendencia espacial significativa, se aplicó kriging universal para realizar la interpolación de la temperatura.

  3. Los mapas generados permiten observar la distribución espacial estimada de la temperatura y la incertidumbre asociada a las predicciones. La validación cruzada mostró un desempeño adecuado, con error promedio cercano a cero, RMSE aproximado de 1.01 °C y correlación cercana a 0.82 entre valores observados y predichos.

El análisis demuestra que la temperatura en el área de estudio presenta dependencia espacial suficiente para ser modelada e interpolada mediante kriging.

17 Conclusiones

  • El ejercicio permitió aplicar el flujo completo de un análisis geoestadístico sobre datos espaciales de una finca de aguacate. La selección del primer periodo evitó mezclar mediciones tomadas en diferentes momentos, lo cual es importante porque la temperatura puede variar temporalmente y afectar la interpretación espacial.

  • La creación del geodata y la transformación a coordenadas proyectadas fueron pasos fundamentales para calcular distancias reales entre puntos. Posteriormente, el análisis de tendencia permitió decidir si era más apropiado aplicar kriging ordinario o kriging universal.

  • El semivariograma empírico y el ajuste de modelos teóricos permitieron caracterizar la estructura de dependencia espacial de la temperatura. Finalmente, el kriging generó una superficie continua de predicción, útil para visualizar zonas más cálidas o más frías dentro del área analizada.