# Sentencia para procesarlo
# ·C:\Users\ita-rodricos\AppData\Local\Programs\R\R-4.4.1\bin\x64\Rscript.exe "procesado_masivo.R"


# 1. Cargar librerías 
suppressPackageStartupMessages({
  library(lidR)
  library(terra)
  library(doParallel)
  library(foreach)
})

# 2. Configuración de Rutas
ruta_laz <- "P:/Cartografia/02_Altimetria/023_LIDAR/2025/Todos1"
ruta_mdt <- "T:/MDT2x2" 
ruta_salida <- "P:/Cartografia/02_Altimetria/023_LIDAR/2025/Salida1"

if(!dir.exists(ruta_salida)) {
  dir.create(ruta_salida, recursive = TRUE)
}

# 3. Configurar el clúster de procesamiento paralelo
# ADVERTENCIA: Con 8GB de RAM, NO superar los 3 clústeres. 
# Si el PC se congela, baja este número a 2.
num_cores <- 3
cl <- makeCluster(num_cores)
registerDoParallel(cl)

cat("\n[ INICIO ] Iniciando procesamiento en paralelo con", num_cores, "hilos...\n")

# 4. Leer lista de archivos
archivos_lidar <- list.files(ruta_laz, pattern = "\\.laz$", full.names = TRUE)
total_archivos <- length(archivos_lidar)
cat("Se han detectado", total_archivos, "archivos para procesar.\n\n")

# Función de métricas (SIN Banda 83)
mis_metricas_forestales <- function(X, Y, Z, I, ReturnNumber, Classification) {
  std <- lidR::stdmetrics(X, Y, Z, I, ReturnNumber, Classification)
  total_puntos <- length(Z)
  
  if (total_puntos == 0) {
    std$zentropy <- as.numeric(NA) 
    nuevas_metricas <- list(
      B70_Cob_05m_1R=as.numeric(NA), B71_Cob_2m_1R=as.numeric(NA), B72_Cob_5m_1R=as.numeric(NA),
      B73_Cob_015_05m=as.numeric(NA), B74_Cob_05_15m=as.numeric(NA), B75_Cob_15_3m=as.numeric(NA),
      B76_Cob_3_5m=as.numeric(NA), B77_Cob_5m=as.numeric(NA), B78_Base_Copa=as.numeric(NA), 
      B79_Clase3=as.numeric(NA), B80_Clase4=as.numeric(NA), B81_Clase5=as.numeric(NA)
    )
    return(c(std, nuevas_metricas))
  }
  
  z_sobre_15cm <- Z[Z > 0.15]
  if (length(z_sobre_15cm) > 0) {
    std$zentropy <- as.numeric(entropy(z_sobre_15cm))
  } else {
    std$zentropy <- as.numeric(NA)
  }
  
  idx_1R <- ReturnNumber == 1
  z_1R <- Z[idx_1R]
  total_1R <- length(z_1R)
  
  b70 <- as.numeric(if(total_1R > 0) (sum(z_1R > 0.5) / total_1R) * 100 else 0)
  b71 <- as.numeric(if(total_1R > 0) (sum(z_1R > 2.0) / total_1R) * 100 else 0)
  b72 <- as.numeric(if(total_1R > 0) (sum(z_1R > 5.0) / total_1R) * 100 else 0)
  b73 <- as.numeric((sum(Z >= 0.15 & Z <= 0.5) / total_puntos) * 100)
  b74 <- as.numeric((sum(Z > 0.5 & Z <= 1.5) / total_puntos) * 100)
  b75 <- as.numeric((sum(Z > 1.5 & Z <= 3.0) / total_puntos) * 100)
  b76 <- as.numeric((sum(Z > 3.0 & Z <= 5.0) / total_puntos) * 100)
  b77 <- as.numeric((sum(Z > 5.0) / total_puntos) * 100)
  
  p95 <- quantile(Z, probs = 0.95, na.rm = TRUE)
  if (!is.na(p95) && p95 > 4) {
    z_mayor_2 <- Z[Z > 2.0]
    b78 <- as.numeric(if(length(z_mayor_2) > 0) quantile(z_mayor_2, probs = 0.20, na.rm = TRUE) else NA)
  } else {
    b78 <- as.numeric(NA)
  }
  
  b79 <- as.numeric((sum(Classification == 3) / total_puntos) * 100)
  b80 <- as.numeric((sum(Classification == 4) / total_puntos) * 100)
  b81 <- as.numeric((sum(Classification == 5) / total_puntos) * 100)
  
  nuevas_metricas <- list(
    B70_Cob_05m_1R = b70, B71_Cob_2m_1R = b71, B72_Cob_5m_1R = b72,
    B73_Cob_015_05m = b73, B74_Cob_05_15m = b74, B75_Cob_15_3m = b75,
    B76_Cob_3_5m = b76, B77_Cob_5m = b77, B78_Base_Copa = b78, 
    B79_Clase3 = b79, B80_Clase4 = b80, B81_Clase5 = b81
  )
  
  return(c(std, nuevas_metricas))
}

# 5. BUCLE PARALELO (Reparte los archivos entre los núcleos asignados)
# El paquete foreach aísla la memoria de cada núcleo automáticamente
foreach(archivo = archivos_lidar, .packages = c("lidR", "terra")) %dopar% {
  
  nombre_base <- tools::file_path_sans_ext(basename(archivo))
  archivo_mdt <- file.path(ruta_mdt, paste0(nombre_base, ".tif"))
  archivo_salida <- file.path(ruta_salida, paste0(nombre_base, "_Estadisticas.tif"))
  
  # Si el archivo ya se procesó en el pasado, lo salta (ideal si se apaga el PC y quieres reanudar)
  if (file.exists(archivo_salida)) { return(NULL) }
  if (!file.exists(archivo_mdt)) { return(NULL) }
  
  las <- readLAS(archivo, filter = "-drop_class 6 7 12")
  if (is.empty(las)) { return(NULL) }
  
  mdt <- rast(archivo_mdt)
  las_norm <- normalize_height(las, mdt)
  
  origen_x <- 165000
  origen_y <- 4415300
  bbox <- ext(las_norm)
  xmin_alineado <- origen_x + floor((bbox$xmin - origen_x) / 10) * 10
  xmax_alineado <- origen_x + ceiling((bbox$xmax - origen_x) / 10) * 10
  ymin_alineado <- origen_y + floor((bbox$ymin - origen_y) / 10) * 10
  ymax_alineado <- origen_y + ceiling((bbox$ymax - origen_y) / 10) * 10
  
  plantilla_10x10 <- rast(xmin=xmin_alineado, xmax=xmax_alineado, 
                          ymin=ymin_alineado, ymax=ymax_alineado, res=10)
  
  metricas <- pixel_metrics(las_norm, 
                            ~mis_metricas_forestales(X, Y, Z, Intensity, ReturnNumber, Classification), 
                            res = plantilla_10x10)
  
  writeRaster(metricas, filename = archivo_salida, overwrite = TRUE)
  
  # Forzar limpieza agresiva de RAM del hilo de trabajo
  rm(las, mdt, las_norm, plantilla_10x10, metricas)
  gc()
}

# Apagar el clúster al finalizar
stopCluster(cl)
cat("\n[ FIN ] Proceso masivo completado.\n")