#Forma de lanzarlo.
#C:\Users\ita-rodricos\AppData\Local\Programs\R\R-4.4.1\bin\x64\Rscript.exe "procesado_masivo_todos1-3.169-4nucleos.R"

# 1. Cargar librerias
suppressPackageStartupMessages({
  library(lidR)
  library(terra)
  library(parallel)
})

# 2. Configuracion 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)
}

archivo_log_principal <- file.path(ruta_salida, "REGISTRO_PRINCIPAL.txt")
cat(paste(Sys.time(), "- Iniciando proceso...\n"), file = archivo_log_principal)

# 3. Leer lista de archivos
archivos_lidar <- list.files(ruta_laz, pattern = "\\.laz$", full.names = TRUE)
total_archivos <- length(archivos_lidar)
cat(paste("Archivos a procesar:", total_archivos, "\n"), file = archivo_log_principal, append = TRUE)

# 4. Funcion de metricas (Debe estar autocontenida para el nuevo cluster)
mis_metricas_forestales <- function(X, Y, Z, I, ReturnNumber, Classification) {
  # lidR::stdmetrics calcula variables predeterminadas como zmax, zmean, zentropy, etc.
  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 por defecto captura la entropia de toda la distribucion de alturas[cite: 2]
    # Aqui lo forzamos solo para retornos > 15 cm.
    std$zentropy <- as.numeric(lidR::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 <- stats::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) stats::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. Definir la funcion trabajadora para un solo archivo
procesar_un_archivo <- function(archivo) {
  
  # Variables requeridas en el subproceso
  ruta_mdt_local <- "T:/MDT2x2"
  ruta_salida_local <- "P:/Cartografia/02_Altimetria/023_LIDAR/2025/Salida1"
  
  nombre_base <- tools::file_path_sans_ext(basename(archivo))
  archivo_mdt <- file.path(ruta_mdt_local, paste0(nombre_base, ".tif"))
  archivo_salida <- file.path(ruta_salida_local, paste0(nombre_base, "_Estadisticas.tif"))
  archivo_log <- file.path(ruta_salida_local, paste0("LOG_", nombre_base, ".txt"))
  
  # Evitar re-procesos
  if (file.exists(archivo_salida)) { return(NULL) }
  
  # Iniciar registro individual
  cat("Iniciando archivo\n", file = archivo_log)
  
  if (!file.exists(archivo_mdt)) { 
    cat("MDT no encontrado\n", file = archivo_log, append = TRUE)
    return(NULL) 
  }
  
  # Cargar dependencias
  require(lidR)
  require(terra)
  set_lidr_threads(1)
  
  cat("Leyendo LAS...\n", file = archivo_log, append = TRUE)
  las <- readLAS(archivo, filter = "-drop_class 6 7 12")
  
  if (is.empty(las)) { 
    cat("LAS vacio\n", file = archivo_log, append = TRUE)
    return(NULL) 
  }
  
  cat("Normalizando...\n", file = archivo_log, append = TRUE)
  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)
  
  cat("Calculando metricas...\n", file = archivo_log, append = TRUE)
  metricas <- pixel_metrics(las_norm, 
                            ~mis_metricas_forestales(X, Y, Z, Intensity, ReturnNumber, Classification), 
                            res = plantilla_10x10)
  
  cat("Escribiendo TIF...\n", file = archivo_log, append = TRUE)
  writeRaster(metricas, filename = archivo_salida, overwrite = TRUE)
  
  cat("FINALIZADO OK\n", file = archivo_log, append = TRUE)
  
  rm(las, mdt, las_norm, plantilla_10x10, metricas)
  gc()
  return(NULL)
}

# 6. Lanzar cluster nativo (Seguro contra cuelgues)
cat(paste(Sys.time(), "- Creando Cluster...\n"), file = archivo_log_principal, append = TRUE)

# PSOCK es mas robusto en Windows
num_cores <- 4
cl <- makePSOCKcluster(num_cores)

# Exportar variables y funciones al cluster
clusterExport(cl, varlist = c("mis_metricas_forestales", "procesar_un_archivo"))

cat(paste(Sys.time(), "- Repartiendo trabajo a los nodos...\n"), file = archivo_log_principal, append = TRUE)

# Load Balancing: Asigna los archivos a medida que los nucleos quedan libres
parLapplyLB(cl, archivos_lidar, procesar_un_archivo)

stopCluster(cl)
cat(paste(Sys.time(), "- PROCESO COMPLETADO\n"), file = archivo_log_principal, append = TRUE)