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


# 1. Cargar librerias principales (Sin librerias de paralelismo)
suppressPackageStartupMessages({
  library(lidR)
  library(terra)
})

# Forzar a lidR a usar 1 solo hilo interno para proteger la RAM de 8GB
set_lidr_threads(1)

# 2. Configuracion de Rutas
ruta_laz <- "P:/Cartografia/02_Altimetria/023_LIDAR/2025/Todos2"
ruta_mdt <- "T:/MDT2x2" 
ruta_salida <- "P:/Cartografia/02_Altimetria/023_LIDAR/2025/Salida2"

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 SECUENCIAL BLINDADO...\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
mis_metricas_forestales <- function(X, Y, Z, I, ReturnNumber, Classification) {
  # Calculo de variables predeterminadas (zmax, zmean, zentropy, etc.)[cite: 9]
  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(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. Bucle Secuencial de Procesamiento
contador <- 1

for (archivo in archivos_lidar) {
  
  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 ya existe, saltar
  if (file.exists(archivo_salida)) { 
    contador <- contador + 1
    next 
  }
  
  # El tryCatch evita que la ejecución se detenga si un archivo falla
  tryCatch({
    
    if (!file.exists(archivo_mdt)) { 
      cat(paste(Sys.time(), "-", nombre_base, "- ERROR: MDT NO ENCONTRADO\n"), file = archivo_log_principal, append = TRUE)
      stop("MDT faltante")
    }
    
    las <- readLAS(archivo, filter = "-drop_class 6 7 12")
    if (is.empty(las)) { 
      cat(paste(Sys.time(), "-", nombre_base, "- AVISO: LAS VACIO TRAS FILTRO\n"), file = archivo_log_principal, append = TRUE)
      stop("Archivo vacio")
    }
    
    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)
    
    cat(paste(Sys.time(), "-", nombre_base, "- OK (", contador, "/", total_archivos, ")\n"), file = archivo_log_principal, append = TRUE)
    
  }, error = function(e) {
    if (e$message != "MDT faltante" && e$message != "Archivo vacio") {
      cat(paste(Sys.time(), "-", nombre_base, "- ERROR CRITICO:", e$message, "\n"), file = archivo_log_principal, append = TRUE)
    }
  })
  
  # Limpieza agresiva de memoria al final de cada iteración
  rm(list = intersect(ls(), c("las", "mdt", "las_norm", "plantilla_10x10", "metricas")))
  gc()
  
  contador <- contador + 1
}

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