# 1. Cargar la librería necesaria
library(terra)

# 2. Definir los directorios y archivos
dir_in <- "T:/MDT2x2"
dir_out <- "P:/Cartografia/02_Altimetria/023_LIDAR/2025/MDT_50k"
archivo_txt <- "hojas50k.txt" # Asegúrate de que esté en tu directorio de trabajo o pon la ruta completa

# Crear la carpeta de salida si no existe
if(!dir.exists(dir_out)) {
  dir.create(dir_out, recursive = TRUE)
}

# 3. Leer y procesar el archivo de texto de las hojas 50k
lineas <- readLines(archivo_txt)
lineas <- lineas[lineas != ""] # Eliminar líneas vacías

# Crear una tabla vacía para guardar las coordenadas
hojas_df <- data.frame(Nombre=character(), Xmin=numeric(), Ymin=numeric(), 
                       Xmax=numeric(), Ymax=numeric(), stringsAsFactors=FALSE)

cat("Leyendo las coordenadas de las hojas 50k...\n")
for(linea in lineas){
  # Reemplazar cualquier dos puntos o punto y coma por comas para estandarizar
  linea_limpia <- gsub("[:;]", ",", linea)
  partes <- unlist(strsplit(linea_limpia, ","))
  
  # Si la línea tiene los 5 datos esperados, se añade a la tabla
  if(length(partes) >= 5){
    hojas_df <- rbind(hojas_df, data.frame(
      Nombre = trimws(partes[1]),
      Xmin   = as.numeric(trimws(partes[2])),
      Ymin   = as.numeric(trimws(partes[3])),
      Xmax   = as.numeric(trimws(partes[4])),
      Ymax   = as.numeric(trimws(partes[5]))
    ))
  } else {
    cat("   [Aviso] Revisar formato de la línea (faltan separadores):", linea, "\n")
  }
}

# 4. Crear el Mosaico Virtual (VRT) de todos los archivos 1x1 km
cat("\nBuscando archivos 1x1 en", dir_in, "...\n")
archivos_1x1 <- list.files(dir_in, pattern = "\\.tif$", full.names = TRUE, ignore.case = TRUE)

if(length(archivos_1x1) == 0) {
  stop("¡Alto! No se encontraron archivos .tif en T:/MDT2x2")
}
cat("Se encontraron", length(archivos_1x1), "archivos. Creando catálogo virtual...\n")

# Se crea el catálogo virtual. Es súper rápido y apenas ocupa espacio.
vrt_path <- file.path(dir_in, "catalogo_global.vrt")
vrt_global <- vrt(archivos_1x1, vrt_path, overwrite = TRUE)

# 5. Bucle de Procesamiento por Hoja 50k
cat("\nIniciando el procesamiento de", nrow(hojas_df), "hojas 50k...\n")

# Definir la ventana de suavizado (3x3 para difuminar las uniones 1x1)
ventana_suavizado <- matrix(1/9, nrow = 3, ncol = 3)

for (i in 1:nrow(hojas_df)) {
  nombre_hoja <- hojas_df$Nombre[i]
  cat("-> Procesando hoja 50k:", nombre_hoja, "...\n")
  
  # Definir la caja (extensión) exacta de la hoja 50.000
  ext_50k_exacta <- ext(hojas_df$Xmin[i], hojas_df$Xmax[i], hojas_df$Ymin[i], hojas_df$Ymax[i])
  
  # Añadir un BUFFER de 20 metros alrededor de la caja exacta
  # Esto garantiza que el suavizado de los bordes funcione perfectamente
  ext_buffer <- ext(
    hojas_df$Xmin[i] - 20, 
    hojas_df$Xmax[i] + 20, 
    hojas_df$Ymin[i] - 20, 
    hojas_df$Ymax[i] + 20
  )
  
  # A) Extraer la porción del VRT global con ese buffer de seguridad
  mdt_hoja_ampliada <- crop(vrt_global, ext_buffer)
  
  # B) Homogeneizar y suavizar todas las uniones de los 1x1 km internos
  mdt_suavizado <- focal(mdt_hoja_ampliada, w = ventana_suavizado, fun = mean, na.rm = TRUE)
  
  # C) Recortar el margen de seguridad para dejar la hoja a su medida 50k estricta
  mdt_final <- crop(mdt_suavizado, ext_50k_exacta)
  
  # D) Guardar el resultado
  ruta_archivo_salida <- file.path(dir_out, paste0("MDT_50k_", nombre_hoja, ".tif"))
  writeRaster(mdt_final, filename = ruta_archivo_salida, overwrite = TRUE)
}

cat("\n¡Fusión y homogeneización de Hojas 50k completadas con éxito!\n")