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

# Forzar a lidR a usar 1 solo hilo interno (seguridad para RAM)
set_lidr_threads(1)

# 2. Configuración de Rutas Fijas (USANDO BARRAS NORMALES)
ruta_laz <- "P:/Cartografia/02_Altimetria/023_LIDAR/2025/Todos2/PNOA_Ejemplo.laz"
ruta_mdt <- "T:/MDT2x2/PNOA_Ejemplo.tif" 
ruta_salida <- "P:/Cartografia/02_Altimetria/023_LIDAR/2025/Salida2/PNOA_Ejemplo_Estadisticas.tif"

# 3. Función de métricas (68 variables)
mis_metricas_forestales <- function(X, Y, Z, I, ReturnNumber, Classification) {
  
  # Calcula métricas base predeterminadas (zmax, zmean, zentropy, etc.)[cite: 12]
  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))
  }
  
  # Corrección: Zentropy solo para retornos > 15 cm
  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))
}

# 4. Comprobaciones de seguridad iniciales
if (!file.exists(ruta_laz)) stop(paste("El archivo LAZ no existe en:", ruta_laz))
if (!file.exists(ruta_mdt)) stop(paste("El archivo MDT no existe en:", ruta_mdt))
if (!dir.exists(dirname(ruta_salida))) dir.create(dirname(ruta_salida), recursive = TRUE)

# 5. Procesamiento del archivo
cat("\nIniciando proceso para:", basename(ruta_laz), "\n")
cat("1. Leyendo y filtrando nube (excluyendo clases 6, 7 y 12)...\n")
las <- readLAS(ruta_laz, filter = "-drop_class 6 7 12")

if (is.empty(las)) stop("El archivo está vacío tras aplicar el filtro de clases.")

cat("2. Leyendo MDT y normalizando alturas...\n")
mdt <- rast(ruta_mdt)
las_norm <- normalize_height(las, mdt)

cat("3. Creando malla de rasterización 10x10...\n")
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("4. Calculando 68 bandas estadísticas...\n")
metricas <- pixel_metrics(las_norm, ~mis_metricas_forestales(X, Y, Z, Intensity, ReturnNumber, Classification), res = plantilla_10x10)

cat("5. Escribiendo archivo de salida TIF...\n")
writeRaster(metricas, filename = ruta_salida, overwrite = TRUE)

cat("\n¡PROCESO COMPLETADO!\nEl archivo se ha guardado en:", ruta_salida, "\n")