# 1. Cargar la librería
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_210.txt" 

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 != ""] 

hojas_df <- data.frame(Nombre=character(), Xmin=numeric(), Ymin=numeric(), 
                       Xmax=numeric(), Ymax=numeric(), stringsAsFactors=FALSE)

cat("Leyendo las coordenadas de", archivo_txt, "...\n")
for(linea in lineas){
  linea_limpia <- gsub("[:;]", ",", linea)
  partes <- unlist(strsplit(linea_limpia, ","))
  
  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]))
    ))
  }
}

# 4. Crear el Mosaico Virtual (VRT) específico para este listado
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")
}

# Se nombra el VRT con el sufijo 210 para evitar conflictos
vrt_path <- file.path(dir_in, "catalogo_global_210.vrt")
cat("Creando mosaico virtual propio:", vrt_path, "...\n")
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 del listado 210...\n")
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")
  
  ext_50k_exacta <- ext(hojas_df$Xmin[i], hojas_df$Xmax[i], hojas_df$Ymin[i], hojas_df$Ymax[i])
  
  ext_buffer <- ext(
    hojas_df$Xmin[i] - 20, 
    hojas_df$Xmax[i] + 20, 
    hojas_df$Ymin[i] - 20, 
    hojas_df$Ymax[i] + 20
  )
  
  mdt_hoja_ampliada <- crop(vrt_global, ext_buffer)
  mdt_suavizado <- focal(mdt_hoja_ampliada, w = ventana_suavizado, fun = mean, na.rm = TRUE)
  mdt_final <- crop(mdt_suavizado, ext_50k_exacta)
  
  ruta_archivo_salida <- file.path(dir_out, paste0("MDT_50k_", nombre_hoja, ".tif"))
  writeRaster(mdt_final, filename = ruta_archivo_salida, overwrite = TRUE)
  
  rm(mdt_hoja_ampliada, mdt_suavizado, mdt_final)
  gc()
}

cat("\n¡Proceso del listado 210 completado!\n")