# 1. Cargar librerías necesarias
library(lidR)
library(terra)

# 2. Definir rutas (Origen múltiple y Destino único)
carpetas_origen <- c(
  "P:/Cartografia/02_Altimetria/023_LIDAR/2025/Todos0",
  "P:/Cartografia/02_Altimetria/023_LIDAR/2025/Todos1",
  "P:/Cartografia/02_Altimetria/023_LIDAR/2025/Todos2"
)
ruta_salida <- "P:/Cartografia/02_Altimetria/023_LIDAR/Resultados"

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

# 3. Crear el catálogo continuo y aplicar configuración
# Lee las 3 carpetas a la vez para procesar los datos por bloques[cite: 1]
ctg <- readLAScatalog(carpetas_origen)

# Filtramos para quedarnos solo con el suelo (2) y quitar solapes[cite: 1]
opt_filter(ctg) <- "-keep_class 2 -drop_overlap"

# 4. LECTURA ANTIBALAS DEL ARCHIVO CSV
# Leemos el archivo asegurando que el separador es el punto y coma
hojas_25k <- read.csv("25k.csv", sep = ";", strip.white = TRUE) 

# Obligamos a R a renombrar las columnas para evitar caracteres invisibles
colnames(hojas_25k) <- c("X_min", "Y_min", "X_max", "Y_max", "HOJA25K")

cat("Se van a procesar", nrow(hojas_25k), "hojas cartográficas.\n\n")

# 5. BUCLE: Procesar cada hoja del CSV
for (i in 1:nrow(hojas_25k)) {
  
  # Asignamos las variables (ahora es 100% seguro que R las encontrará)
  nombre_hoja <- hojas_25k$HOJA25K[i] 
  xmin <- hojas_25k$X_min[i]
  ymin <- hojas_25k$Y_min[i]
  xmax <- hojas_25k$X_max[i]
  ymax <- hojas_25k$Y_max[i]
  
  cat("-> Recortando y procesando hoja:", nombre_hoja, "...\n")
  
  # A. Recortar el área de la hoja CON UN BÚFER DE 20 METROS
  # Cogemos puntos extra para que los algoritmos calculen bien los bordes
  las_hoja <- clip_rectangle(ctg, 
                             xleft = xmin - 20, ybottom = ymin - 20, 
                             xright = xmax + 20, ytop = ymax + 20)
  
  # Si en esas coordenadas no hay datos LiDAR, saltamos a la siguiente hoja
  if (is.empty(las_hoja)) {
    cat("   [Aviso] No hay datos en esta extensión. Saltando...\n")
    next
  }
  
  # B. Crear MDT y suavizar
  mdt_base <- rasterize_terrain(las_hoja, res = 2, algorithm = tin())
  ventana_suavizado <- matrix(1, nrow = 3, ncol = 3)
  mdt_suavizado <- focal(mdt_base, w = ventana_suavizado, fun = median, na.rm = TRUE)
  
  # C. Recortar el raster a la medida EXACTA de la hoja (eliminando el búfer)
  extension_exacta <- ext(xmin, xmax, ymin, ymax)
  mdt_final <- crop(mdt_suavizado, extension_exacta)
  
  # D. Guardar el archivo definitivo en la carpeta Resultados
  archivo_salida <- file.path(ruta_salida, paste0(nombre_hoja, "_MDT_25k.tif"))
  writeRaster(mdt_final, filename = archivo_salida, overwrite = TRUE)
  cat("   [OK] Guardado en:", archivo_salida, "\n\n")
  
  # Liberar memoria RAM antes de pasar a la siguiente hoja
  rm(las_hoja, mdt_base, mdt_suavizado, mdt_final)
  gc()
}

cat("\n¡TODAS LAS HOJAS 25K HAN SIDO PROCESADAS PERFECTAMENTE!\n")