1. Objetivo y alcance

Este documento evalua la estructura espacial de Temperature y construye una superficie interpolada mediante kriging. El flujo es:

  1. Analisis exploratorio de la variable.
  2. Creacion de datos espaciales proyectados.
  3. Diagnostico de tendencia espacial.
  4. Autocorrelacion espacial con semivariograma experimental.
  5. Ajuste y seleccion de un modelo teorico.
  6. Prediccion espacial y validacion cruzada.

El informe esta configurado para que el codigo quede oculto por defecto, pero disponible desde el boton Code de cada bloque. Asi se conserva la trazabilidad tecnica sin interrumpir la lectura principal.

2. Carga y preparacion de datos

datos_raw <- read_excel(path = params$archivo, sheet = params$hoja)

columnas_necesarias <- c(params$lat_col, params$lon_col, params$variable_objetivo)
faltantes <- setdiff(columnas_necesarias, names(datos_raw))
if (length(faltantes) > 0) {
  stop("No se encontraron estas columnas en el archivo: ", paste(faltantes, collapse = ", "))
}

datos_base <- datos_raw %>%
  mutate(
    lon = as.numeric(.data[[params$lon_col]]),
    lat = as.numeric(.data[[params$lat_col]]),
    valor = as.numeric(.data[[params$variable_objetivo]])
  ) %>%
  filter(is.finite(lon), is.finite(lat), is.finite(valor))

datos_geo <- datos_base %>%
  group_by(lon, lat) %>%
  summarise(
    valor = mean(valor, na.rm = TRUE),
    sd_intrapunto = sd(valor, na.rm = TRUE),
    n_observaciones = dplyr::n(),
    .groups = "drop"
  ) %>%
  mutate(
    sd_intrapunto = ifelse(is.na(sd_intrapunto), 0, sd_intrapunto),
    id_punto = row_number()
  )

resumen_base <- tibble::tibble(
  indicador = c(
    "Registros originales",
    "Registros validos",
    "Puntos espaciales unicos",
    "Variable analizada",
    "Promedio espacial",
    "Desviacion estandar espacial",
    "Minimo espacial",
    "Maximo espacial"
  ),
  valor = c(
    nrow(datos_raw),
    nrow(datos_base),
    nrow(datos_geo),
    params$variable_objetivo,
    fmt_num(mean(datos_geo$valor), 2),
    fmt_num(sd(datos_geo$valor), 2),
    fmt_num(min(datos_geo$valor), 2),
    fmt_num(max(datos_geo$valor), 2)
  )
)

kable(resumen_base, caption = "Resumen inicial de la informacion usada en el analisis espacial")
Resumen inicial de la informacion usada en el analisis espacial
indicador valor
Registros originales 20271
Registros validos 20271
Puntos espaciales unicos 1685
Variable analizada Temperature
Promedio espacial 23.06
Desviacion estandar espacial 1.22
Minimo espacial 20.02
Maximo espacial 25.8

Variable objetivo

Temperature
Se analiza la media por coordenada para controlar mediciones repetidas.

Puntos espaciales

1,685
Registros validos originales: 20,271.

Rango observado

20.02 a 25.8
Media espacial: 23.06.

3. Analisis exploratorio

ggplot(datos_geo, aes(x = valor)) +
  geom_histogram(bins = 35, fill = "#2563eb", color = "white", alpha = 0.9) +
  labs(
    title = paste("Distribucion espacial de", params$variable_objetivo),
    x = params$variable_objetivo,
    y = "Frecuencia"
  )

ggplot(datos_geo, aes(x = n_observaciones)) +
  geom_histogram(bins = 30, fill = "#0f766e", color = "white", alpha = 0.9) +
  labs(
    title = "Numero de mediciones por punto espacial",
    x = "Mediciones agregadas por coordenada",
    y = "Frecuencia"
  )

ggplot(datos_geo, aes(x = lon, y = lat, color = valor)) +
  geom_point(size = 1.8, alpha = 0.85) +
  scale_color_viridis_c(option = "C") +
  coord_equal() +
  labs(
    title = paste("Distribucion espacial observada de", params$variable_objetivo),
    x = "Longitud",
    y = "Latitud",
    color = params$variable_objetivo
  )

4. Creacion de geodata

Para medir distancias en metros, las coordenadas geograficas WGS84 se transforman a UTM zona 18N (EPSG:32618). Esta proyeccion es adecuada para el suroccidente de Colombia y evita calcular semivariogramas en grados decimales.

puntos_sf <- st_as_sf(
  datos_geo,
  coords = c("lon", "lat"),
  crs = params$crs_origen,
  remove = FALSE
)

puntos_utm <- st_transform(puntos_sf, params$crs_proyectado)
xy <- st_coordinates(puntos_utm)
puntos_utm$x <- xy[, 1]
puntos_utm$y <- xy[, 2]

puntos_sp <- as(puntos_utm, "Spatial")

bbox_utm <- st_bbox(puntos_utm)
extension_m <- c(
  ancho = as.numeric(bbox_utm["xmax"] - bbox_utm["xmin"]),
  alto = as.numeric(bbox_utm["ymax"] - bbox_utm["ymin"])
)

kable(
  data.frame(
    medida = names(extension_m),
    metros = fmt_num(unname(extension_m), 1)
  ),
  caption = "Extension aproximada del area de estudio en metros"
)
Extension aproximada del area de estudio en metros
medida metros
ancho 11,686.1
alto 8,541.4

5. Diagnostico de tendencia espacial

Antes de interpolar, se prueba si la variable presenta una tendencia global respecto a las coordenadas proyectadas. Si la tendencia es relevante, se usa kriging universal; si no, kriging ordinario.

modelo_tendencia <- lm(valor ~ x + y, data = puntos_utm)
resumen_tendencia <- summary(modelo_tendencia)
anova_tendencia <- anova(modelo_tendencia)

r2_adj <- resumen_tendencia$adj.r.squared
p_tendencia <- pf(
  resumen_tendencia$fstatistic[1],
  resumen_tendencia$fstatistic[2],
  resumen_tendencia$fstatistic[3],
  lower.tail = FALSE
)

usar_tendencia <- is.finite(p_tendencia) && p_tendencia < 0.05 && r2_adj >= 0.05
formula_kriging <- if (usar_tendencia) valor ~ x + y else valor ~ 1
tipo_kriging <- if (usar_tendencia) "Kriging universal" else "Kriging ordinario"

kable(
  data.frame(
    metrica = c("R2 ajustado", "p-valor global", "Decision", "Formula usada"),
    valor = c(
      fmt_num(r2_adj, 4),
      fmt_num(p_tendencia, 5),
      ifelse(usar_tendencia, "Tendencia espacial relevante", "Sin tendencia global fuerte"),
      deparse(formula_kriging)
    )
  ),
  caption = "Diagnostico de tendencia espacial"
)
Diagnostico de tendencia espacial
metrica valor
R2 ajustado 0.5751
p-valor global 0
Decision Tendencia espacial relevante
Formula usada valor ~ x + y 
ggplot(puntos_utm, aes(x = x, y = valor)) +
  geom_point(alpha = 0.45, color = "#2563eb") +
  geom_smooth(method = "lm", color = "#b45309", se = TRUE) +
  labs(
    title = paste("Tendencia oeste-este de", params$variable_objetivo),
    x = "Coordenada X proyectada (m)",
    y = params$variable_objetivo
  )

ggplot(puntos_utm, aes(x = y, y = valor)) +
  geom_point(alpha = 0.45, color = "#2563eb") +
  geom_smooth(method = "lm", color = "#b45309", se = TRUE) +
  labs(
    title = paste("Tendencia sur-norte de", params$variable_objetivo),
    x = "Coordenada Y proyectada (m)",
    y = params$variable_objetivo
  )

Tendencia espacial

Detectada
R2 ajustado = 0.575; p-valor = 0.

Metodo seleccionado

Kriging universal
Formula para semivariograma y prediccion: valor ~ x + y.

6. Autocorrelacion espacial: semivariograma experimental

El semivariograma resume como cambia la diferencia entre observaciones al aumentar la distancia. Si la semivarianza crece y luego se estabiliza, hay estructura espacial util para interpolar.

distancias <- spDists(coordinates(puntos_sp), longlat = FALSE)
max_dist <- max(distancias[upper.tri(distancias)], na.rm = TRUE)
cutoff <- max_dist / 2
width <- cutoff / 15

vg_emp <- variogram(
  formula_kriging,
  data = puntos_sp,
  cutoff = cutoff,
  width = width
)

kable(
  head(vg_emp, 12),
  digits = 4,
  caption = "Primeras clases de distancia del semivariograma experimental"
)
Primeras clases de distancia del semivariograma experimental
np dist gamma dir.hor dir.ver id
364968 78.0606 0.3593 0 0 var1
39487 1177.2090 1.5337 0 0 var1
128751 1312.1178 1.4578 0 0 var1 
ggplot(vg_emp, aes(x = dist, y = gamma)) +
  geom_point(size = 2.4, color = "#2563eb") +
  geom_line(color = "#64748b") +
  labs(
    title = paste("Semivariograma experimental de", params$variable_objetivo),
    x = "Distancia (m)",
    y = "Semivarianza"
  )

7. Ajuste del modelo teorico

Se comparan modelos teoricos comunes: esferico, exponencial, gaussiano y Matern. El mejor modelo se elige por el menor error de ajuste (SSErr) entre los modelos validos.

var_obj <- var(puntos_utm$valor, na.rm = TRUE)
nugget_ini <- var_obj * 0.10
psill_ini <- var_obj * 0.90
range_ini <- cutoff / 3

modelos_candidatos <- list(
  Esferico = vgm(psill = psill_ini, model = "Sph", range = range_ini, nugget = nugget_ini),
  Exponencial = vgm(psill = psill_ini, model = "Exp", range = range_ini, nugget = nugget_ini),
  Gaussiano = vgm(psill = psill_ini, model = "Gau", range = range_ini, nugget = nugget_ini),
  Matern = vgm(psill = psill_ini, model = "Mat", range = range_ini, nugget = nugget_ini, kappa = 0.5)
)

ajustar_modelo <- function(nombre, modelo_inicial) {
  tryCatch({
    fit <- fit.variogram(vg_emp, model = modelo_inicial, fit.method = 6)
    sse <- attr(fit, "SSErr")
    sill_total <- sum(fit$psill, na.rm = TRUE)
    data.frame(
      modelo = nombre,
      codigo = as.character(fit$model[2]),
      nugget = fit$psill[1],
      partial_sill = fit$psill[2],
      sill = sill_total,
      rango = fit$range[2],
      prop_nugget = fit$psill[1] / sill_total,
      SSErr = sse,
      valido = all(is.finite(c(fit$psill, fit$range[2], sse))) &&
        all(fit$psill >= 0) &&
        fit$range[2] > 0,
      mensaje = "",
      stringsAsFactors = FALSE
    )
  }, error = function(e) {
    data.frame(
      modelo = nombre,
      codigo = NA_character_,
      nugget = NA_real_,
      partial_sill = NA_real_,
      sill = NA_real_,
      rango = NA_real_,
      prop_nugget = NA_real_,
      SSErr = Inf,
      valido = FALSE,
      mensaje = conditionMessage(e),
      stringsAsFactors = FALSE
    )
  })
}

tabla_ajustes <- bind_rows(Map(ajustar_modelo, names(modelos_candidatos), modelos_candidatos)) %>%
  arrange(SSErr)

if (!any(tabla_ajustes$valido)) {
  stop("Ningun modelo teorico produjo un ajuste valido. Revise datos, tendencia o parametros iniciales.")
}

mejor_nombre <- tabla_ajustes$modelo[tabla_ajustes$valido][1]
mejor_modelo <- fit.variogram(
  vg_emp,
  model = modelos_candidatos[[mejor_nombre]],
  fit.method = 6
)

tabla_ajustes %>%
  mutate(
    nugget = round(nugget, 4),
    partial_sill = round(partial_sill, 4),
    sill = round(sill, 4),
    rango = round(rango, 2),
    prop_nugget = round(prop_nugget, 3),
    SSErr = round(SSErr, 6)
  ) %>%
  kable(caption = "Comparacion de modelos teoricos de semivariograma")
Comparacion de modelos teoricos de semivariograma
modelo codigo nugget partial_sill sill rango prop_nugget SSErr valido mensaje
Matern Mat 0.0000 1.5146 1.5146 284.16 0.000 0.003635 TRUE
Exponencial Exp 0.0000 1.5146 1.5146 284.16 0.000 0.003635 TRUE
Esferico Sph 0.2166 1.2256 1.4421 852.60 0.150 0.009271 TRUE
Gaussiano Gau 0.0697 1.1306 1.2003 371.29 0.058 0.235550 TRUE 
vg_line <- variogramLine(mejor_modelo, maxdist = cutoff, n = 200)

ggplot() +
  geom_point(data = vg_emp, aes(x = dist, y = gamma), size = 2.3, color = "#2563eb") +
  geom_line(data = vg_line, aes(x = dist, y = gamma), linewidth = 1, color = "#b45309") +
  labs(
    title = paste("Modelo teorico seleccionado:", mejor_nombre),
    subtitle = paste0("SSErr = ", fmt_num(min(tabla_ajustes$SSErr[tabla_ajustes$valido]), 6)),
    x = "Distancia (m)",
    y = "Semivarianza"
  )

Mejor modelo teorico

Matern
Seleccionado por menor SSErr: 0.003635.

Rango espacial

284.2 m
Distancia aproximada hasta la cual la variable conserva autocorrelacion espacial.

Efecto pepita

0%
Proporcion de variabilidad no explicada por estructura espacial a la escala muestreada.

8. Prediccion espacial por kriging

La prediccion se realiza sobre una grilla regular dentro del poligono convexo de los puntos observados. El tamano de celda se ajusta automaticamente si la grilla supera max_grid_points.

crear_grilla <- function(puntos, cellsize, max_points) {
  hull <- st_convex_hull(st_union(puntos))
  cell_actual <- cellsize
  repeat {
    grilla <- st_make_grid(hull, cellsize = cell_actual, what = "centers") %>%
      st_sf(geometry = .) %>%
      st_filter(hull, .predicate = st_within)

    if (nrow(grilla) <= max_points || cell_actual > cellsize * 10) {
      grilla$cellsize_m <- cell_actual
      return(grilla)
    }
    cell_actual <- cell_actual * 1.25
  }
}

grilla_sf <- crear_grilla(puntos_utm, params$grid_cell_m, params$max_grid_points)
coords_grilla <- st_coordinates(grilla_sf)
grilla_sf$x <- coords_grilla[, 1]
grilla_sf$y <- coords_grilla[, 2]
grilla_sp <- as(grilla_sf, "Spatial")

pred_sp <- krige(
  formula = formula_kriging,
  locations = puntos_sp,
  newdata = grilla_sp,
  model = mejor_modelo
)
## [using universal kriging]
pred_sf <- st_as_sf(pred_sp)
pred_df <- cbind(
  as.data.frame(st_coordinates(pred_sf)),
  st_drop_geometry(pred_sf)
)

puntos_df <- cbind(
  as.data.frame(st_coordinates(puntos_utm)),
  st_drop_geometry(puntos_utm)
)

ggplot() +
  geom_tile(data = pred_df, aes(x = X, y = Y, fill = var1.pred)) +
  geom_point(data = puntos_df, aes(x = X, y = Y), size = 0.35, alpha = 0.35, color = "#0f172a") +
  scale_fill_viridis_c(option = "C") +
  coord_equal() +
  labs(
    title = paste("Prediccion espacial de", params$variable_objetivo),
    subtitle = paste(tipo_kriging, "| Modelo:", mejor_nombre),
    x = "X proyectada (m)",
    y = "Y proyectada (m)",
    fill = paste("Pred.", params$variable_objetivo)
  )

ggplot() +
  geom_tile(data = pred_df, aes(x = X, y = Y, fill = var1.var)) +
  geom_point(data = puntos_df, aes(x = X, y = Y), size = 0.35, alpha = 0.35, color = "#0f172a") +
  scale_fill_viridis_c(option = "magma") +
  coord_equal() +
  labs(
    title = "Incertidumbre de la prediccion",
    subtitle = "Varianza de kriging: valores altos indican menor certeza local.",
    x = "X proyectada (m)",
    y = "Y proyectada (m)",
    fill = "Varianza"
  )

Metodo de prediccion

Kriging universal
Grilla usada: 14,205 celdas; celda = 62.5 m.

Prediccion minima

19.73
Valor minimo interpolado para Temperature.

Prediccion maxima

25.06
Valor maximo interpolado para Temperature.

9. Mapa base tipo streets

Este mapa interactivo usa teselas de OpenStreetMap para ubicar los puntos de muestreo y la superficie predicha sobre una referencia de calles. Requiere el paquete leaflet y conexion a internet al momento de abrir el HTML para cargar el mapa base.

10. Validacion cruzada

Se usa validacion cruzada espacial por pliegues (nfold) con el mismo modelo de semivariograma y la misma formula de kriging. La validacion permite evaluar sesgo y precision del interpolador.

cv <- krige.cv(
  formula = formula_kriging,
  locations = puntos_sp,
  model = mejor_modelo,
  nfold = params$cv_folds
)

cv_df <- as.data.frame(cv)
cv_metricas <- data.frame(
  metrica = c(
    "Error medio (ME)",
    "Error absoluto medio (MAE)",
    "Raiz del error cuadratico medio (RMSE)",
    "RMSE estandarizado",
    "Correlacion observado-predicho"
  ),
  valor = c(
    mean(cv_df$residual, na.rm = TRUE),
    mean(abs(cv_df$residual), na.rm = TRUE),
    sqrt(mean(cv_df$residual^2, na.rm = TRUE)),
    sqrt(mean(cv_df$zscore^2, na.rm = TRUE)),
    cor(cv_df$observed, cv_df$var1.pred, use = "complete.obs")
  )
)

kable(cv_metricas, digits = 4, caption = "Metricas de validacion cruzada")
Metricas de validacion cruzada
metrica valor
Error medio (ME) 0.0029
Error absoluto medio (MAE) 0.2269
Raiz del error cuadratico medio (RMSE) 0.3054
RMSE estandarizado 1.8702
Correlacion observado-predicho 0.9684 
ggplot(cv_df, aes(x = observed, y = var1.pred)) +
  geom_point(alpha = 0.55, color = "#2563eb") +
  geom_abline(slope = 1, intercept = 0, color = "#b45309", linewidth = 1) +
  labs(
    title = "Validacion cruzada: observado vs predicho",
    x = paste("Observado", params$variable_objetivo),
    y = paste("Predicho", params$variable_objetivo)
  )

ggplot(cv_df, aes(x = residual)) +
  geom_histogram(bins = 35, fill = "#64748b", color = "white") +
  geom_vline(xintercept = 0, color = "#b45309", linewidth = 1) +
  labs(
    title = "Distribucion de residuales de validacion",
    x = "Residual: observado - predicho",
    y = "Frecuencia"
  )

RMSE

0.305
Menor RMSE indica mejor precision predictiva.

Sesgo medio

0.003
Valores cercanos a cero indican baja tendencia a sobre/subestimar.

Correlacion CV

0.968
Asociacion entre valores observados y predichos en validacion cruzada.

11. Interpretacion final

Modelo recomendado

Matern
Kriging universal con rango espacial de 284.2 m.

Error relativo

5.3%
RMSE comparado con el rango observado de la variable.

Uso sugerido

Zonificacion exploratoria
La interpretacion debe priorizar patrones amplios y no diferencias puntuales minimas.

Lectura metodologica

El flujo aplicado es coherente con una interpolacion geoestadistica clasica: primero se consolido la informacion por coordenada, luego se evaluo la tendencia espacial, despues se construyo el semivariograma experimental y finalmente se ajustaron modelos teoricos candidatos. Esta secuencia evita interpolar directamente sobre registros repetidos en el mismo punto, lo cual podria inflar artificialmente la informacion espacial disponible.

El mejor modelo teorico fue Matern, con rango aproximado de 284.2 m. La variable presenta una tendencia espacial global suficiente para justificar kriging universal. En terminos practicos, esto significa que parte de la variacion puede explicarse por la posicion general dentro del lote o zona de estudio, no solamente por dependencia local entre puntos vecinos.

Lectura espacial

El efecto pepita es bajo. La estructura espacial capturada por el semivariograma es consistente y favorable para interpolar. El rango espacial indica la distancia aproximada hasta la cual dos puntos tienden a parecerse en la variable analizada. Por encima de esa distancia, la relacion espacial se debilita y el modelo depende menos de la vecindad inmediata.

La superficie de prediccion debe leerse junto con la varianza de kriging: las zonas con mayor varianza no necesariamente tienen valores climaticos extremos, sino mayor incertidumbre por distancia a los puntos muestreados o por menor soporte local del semivariograma.

Lectura de validacion

El RMSE fue 0.305, el MAE fue 0.227 y el sesgo medio fue 0.003. El RMSE equivale al 5.3% del rango observado y el MAE al 3.9%. La validacion cruzada muestra un error bajo frente al rango observado de la variable. La superficie interpolada puede interpretarse con buena confianza relativa, especialmente dentro de zonas con mayor densidad de puntos. El sesgo medio esta cerca de cero, por lo que el modelo no muestra una tendencia marcada a sobreestimar o subestimar en promedio.

Conclusion operativa

El resultado es especialmente util para reconocer gradientes y zonas relativas de mayor o menor Temperature. Para decisiones de manejo, conviene cruzar esta superficie con informacion de altitud, pendiente, cobertura, sanidad del cultivo y fechas de medicion, porque la variable climatica puede variar tanto por posicion espacial como por condiciones temporales de muestreo.

Nota final. Este informe se presenta como un analisis enfocado en temperatura. Si posteriormente se requiere estudiar otra variable climatica, es recomendable generar un informe independiente con la misma ruta metodologica para conservar claridad en los resultados.