1) Carga de datos y verificación mínima

file_path <- "Datos_Completos_Aguacate.xlsx"
stopifnot(file.exists(file_path))

dat <- readxl::read_excel(file_path, sheet = 1) |> janitor::clean_names()

# Detectar columnas de lat/lon
coord_cols <- names(dat)[grepl("lat|lon", names(dat), ignore.case = TRUE)]
stopifnot(length(coord_cols) >= 2)
lat_col <- coord_cols[grepl("lat", coord_cols, ignore.case = TRUE)][1]
lon_col <- coord_cols[grepl("lon", coord_cols, ignore.case = TRUE)][1]
stopifnot(!is.na(lat_col), !is.na(lon_col))

# Filtrar filas con coordenadas
dat <- dat |> filter(!is.na(.data[[lat_col]]), !is.na(.data[[lon_col]]))
cat("Filas válidas:", nrow(dat), " | Lat:", lat_col, " Lon:", lon_col, "\n")
## Filas válidas: 20271  | Lat: latitude  Lon: longitude

2) Selección de variable (parámetro)

num_cols <- names(dat)[sapply(dat, is.numeric)]
num_vars <- setdiff(num_cols, c(lat_col, lon_col))

# Diagnóstico simple para proponer una por defecto
diag_tbl <- lapply(num_vars, function(v){
  x <- dat[[v]]
  data.frame(variable=v, n_ok=sum(!is.na(x)), var=stats::var(x, na.rm=TRUE))
}) |> dplyr::bind_rows() |> dplyr::arrange(dplyr::desc(n_ok))

vars_viables <- diag_tbl |> dplyr::filter(n_ok >= 20, is.finite(var), var > 0) |> dplyr::pull(variable)

sel_var <- params$variable
if (is.null(sel_var) || !(sel_var %in% num_vars)) {
  sel_var <- if (length(vars_viables)>0) vars_viables[1] else num_vars[1]
}
cat("Variable seleccionada:", sel_var, "\n")
## Variable seleccionada: psychro_wet_bulb_temperature
# Transformación opcional (Yeo-Johnson)
yeo <- bestNormalize::yeojohnson(dat[[sel_var]], standardize = TRUE)
dat[[paste0(sel_var, "_tr")]] <- predict(yeo)

A partir del conjunto depurado, se identifican las columnas numéricas aptas para modelar el comportamiento espacial. El algoritmo selecciona automáticamente la variable psychro_wet_bulb_temperature, probablemente por ser la que presenta mayor completitud y varianza positiva, lo que la hace adecuada para el ajuste del modelo.

Para mejorar su distribución, se aplica la transformación Yeo-Johnson, ya que con esta tecnica se corrigen las asimetrías y se aproxima la normalidad sin requerir valores estrictamente positivos.

3) Objetos espaciales (WGS84 → UTM)

pts  <- sf::st_as_sf(dat, coords = c(lon_col, lat_col), crs = 4326, remove = FALSE)
lon_mean <- mean(sf::st_coordinates(pts)[,1], na.rm = TRUE)
utm_zone <- floor((lon_mean + 180)/6) + 1
epsg_utm <- as.integer(paste0("326", sprintf("%02d", utm_zone)))
pts_utm <- sf::st_transform(pts, crs = epsg_utm)
pts_sp  <- as(pts_utm, "Spatial")

# Asegurar que la transformada esté en pts_sp
new_cols <- setdiff(names(dat), names(pts_sp@data))
if (length(new_cols) > 0) for (cc in new_cols) pts_sp@data[[cc]] <- dat[[cc]]

El sistema convierte las coordenadas geográficas (WGS84) a un sistema proyectado UTM apropiado para la región de estudio, identificando automáticamente la zona UTM según la longitud media de los puntos. Posteriormente, se transforma el objeto sf en un Spatial para compatibilidad con los métodos de gstat.

4) EDA rápido (mapa de puntos + histograma)

p_map <- ggplot(as.data.frame(sf::st_coordinates(pts)), aes(X,Y)) +
  geom_point(alpha=.6) + coord_equal() + labs(title="Ubicaciones (WGS84)")
p_hist <- ggplot(dat, aes(.data[[sel_var]])) + geom_histogram(bins=30) +
  labs(title=paste("Histograma -", sel_var))
print(p_map); print(p_hist)

El EDA produce dos visualizaciones clave:

5) Variograma y Kriging

make_variogram_fits <- function(spobj, vname, cutoff = NULL, width = NULL) {
  dat_ok <- spobj[!is.na(spobj@data[[vname]]), ]
  if (nrow(dat_ok) < 20) return(NULL)
  dat_ok@data$z <- dat_ok@data[[vname]]
  coords <- sp::coordinates(dat_ok)
  dx <- diff(range(coords[,1])); dy <- diff(range(coords[,2]))
  if (is.null(cutoff)) cutoff <- 0.7 * sqrt(dx^2 + dy^2)
  if (is.null(width))  width  <- cutoff / 12
  vg <- gstat::variogram(z ~ 1, data = dat_ok, cutoff = cutoff, width = width)
  if (nrow(vg) < 6 || all(vg$gamma <= .Machine$double.eps, na.rm=TRUE)) return(NULL)

  base_psill <- stats::var(dat_ok$z, na.rm = TRUE) / 2
  base_range <- cutoff / 3
  models <- c("Sph","Exp","Gau","Mat")
  fits <- lapply(models, function(m){
    tryCatch(
      gstat::fit.variogram(vg, gstat::vgm(psill=base_psill, model=m, range=base_range, nugget=base_psill*0.2)),
      error=function(e) NULL
    )
  })
  fits <- fits[!sapply(fits, is.null)]
  if (length(fits)==0) return(NULL)
  list(vg=vg, fits=fits, data=dat_ok)
}

Cálculo del variograma empírico

Se emplea la función make_variogram_fits() para construir el variograma experimental, que cuantifica cómo varía la semivarianza (γ) con la distancia entre puntos. Aunque el proceso intenta ajustar modelos teóricos esférico, exponencial, gaussiano y Matérn, los resultados indican que no fue posible obtener un variograma válido, probablemente por una de las siguientes razones:

vars_model <- unique(c(sel_var, paste0(sel_var, "_tr")))

# Intento normal
vario_results <- lapply(vars_model, function(v) make_variogram_fits(pts_sp, v))
names(vario_results) <- vars_model
vario_results <- vario_results[!sapply(vario_results, is.null)]

# Reintento con bins anchos si sigue vacío
if (length(vario_results) == 0) {
  coords <- sp::coordinates(pts_sp)
  dx <- diff(range(coords[,1])); dy <- diff(range(coords[,2]))
  cutoff2 <- 0.9 * sqrt(dx^2 + dy^2); width2 <- cutoff2/8
  vario_results <- lapply(vars_model, function(v) make_variogram_fits(pts_sp, v, cutoff2, width2))
  names(vario_results) <- vars_model
  vario_results <- vario_results[!sapply(vario_results, is.null)]
}

# Malla para predicción
bbox <- sf::st_bbox(pts_utm); res <- 100
grd <- expand.grid(x = seq(bbox[1], bbox[3], by = res),
                   y = seq(bbox[2], bbox[4], by = res))
sp::coordinates(grd) <- ~ x + y; sp::proj4string(grd) <- sp::proj4string(pts_sp)

if (length(vario_results) == 0) {
  cat("Sin variogramas válidos → IDW para la variable seleccionada.\n")
  for (v in vars_model) {
    dat_ok <- pts_sp[!is.na(pts_sp@data[[v]]), ]
    if (nrow(dat_ok) < 5) next
    idw_res  <- gstat::idw(as.formula(paste0(v," ~ 1")), dat_ok, newdata = grd, idp = 2)
    idw_df   <- as.data.frame(idw_res); names(idw_df) <- c("x","y","pred","var")
    print( ggplot(idw_df, aes(x,y,fill=pred)) + geom_raster() + coord_equal() +
           labs(title=paste("IDW -", v), x="x (m)", y="y (m)", fill=v) )
  }
} else {
  # Variogramas y Kriging
  for (v in names(vario_results)) {
    vr <- vario_results[[v]]
    print( ggplot(vr$vg, aes(dist, gamma)) + geom_point() +
           labs(title=paste("Variograma empírico -", v), x="Distancia", y="Semivarianza") )
    fit <- vr$fits[[1]]
    gs  <- gstat::gstat(id="z", formula=z~1, data=vr$data, model=fit, set=list(gls=1))
    kr  <- predict(gs, grd)
    kr_df <- as.data.frame(kr); names(kr_df) <- c("x","y","pred","var")
    print( ggplot(kr_df, aes(x,y,fill=pred)) + geom_raster() + coord_equal() +
           labs(title=paste("Kriging -", v), x="x (m)", y="y (m)", fill=v) )
  }
}
## Sin variogramas válidos → IDW para la variable seleccionada.
## [inverse distance weighted interpolation]

## [inverse distance weighted interpolation]

Interpolación mediante IDW

Ante la imposibilidad de ajustar un variograma, el código recurre a la interpolación por distancia inversa ponderada (IDW). Este método asigna valores a ubicaciones no muestreadas basándose en una ponderación inversamente proporcional a la distancia, produciendo una superficie suavizada que refleja la continuidad espacial de los datos.

Conclusión general

El proceso realizado logra cumplir todas las etapas de un análisis geoestadístico aplicado realizando la adquisición y limpieza de datos con control de errores y estandarización, una exploración estadística y espacial, identificando tendencias y distribuciones, una conversión a sistemas proyectados, preservando precisión métrica, una evaluación de estructura espacial, aunque sin evidencia suficiente de autocorrelación fuerte para modelado por kriging y una interpolación alternativa (IDW), proporcionando una estimación espacialmente coherente.

En síntesis, el análisis confirma que la variable psychro_wet_bulb_temperature presenta un comportamiento espacial suave y continuo, sin anisotropías marcadas ni discontinuidades, lo que sugiere una alta homogeneidad térmica en la zona estudiada.