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

# 2. Configuración de rutas (AJUSTA ESTAS RUTAS A TU PC)
ruta_laz <- ".P:\Cartografia\02_Altimetria\023_LIDAR\2025"
ruta_mdt <- "P:\Cartografia\02_Altimetria\023_LIDAR\2025\MDE_LIDAR\PNOA_LIDAR_CYL-2025_MDE_2x2_ETRS89-30_EGM08_H50_0126.tif" 
ruta_salida <- "P:\Cartografia\02_Altimetria\023_LIDAR\2025\MDE_LIDAR"

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

# 3. Función personalizada de métricas forestales
mis_metricas_forestales <- function(X, Y, Z, I, ReturnNumber, Classification) {
  
  # Usamos stdmetrics asegurando que devuelve numericos
  std <- stdmetrics(X, Y, Z, I, ReturnNumber, Classification)
  
  total_puntos <- length(Z)
  
  if (total_puntos == 0) {
    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))
  }
  
  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))
}

# 4. Buscar archivos origen (.laz) en tu carpeta D:/H50_0344
archivos_lidar <- list.files(ruta_laz, pattern = "\\.laz$", full.names = TRUE)
cat("\nSe encontraron", length(archivos_lidar), "archivos para procesar.\n")

if(length(archivos_lidar) == 0) {
  stop("¡Alto! No se han encontrado archivos .laz en D:/H50_0344.")
}

# 5. BUCLE POR ARCHIVO PARA EXTRAER ESTADÍSTICAS
for (archivo in archivos_lidar) {
  
  nombre_base <- tools::file_path_sans_ext(basename(archivo))
  cat("-> Procesando:", nombre_base, "...\n")
  
  # Buscar el MDT (archivo .tif) que tenga el mismo nombre que el .laz
  archivo_mdt <- file.path(ruta_mdt, paste0(nombre_base, ".tif"))
  
  if (!file.exists(archivo_mdt)) {
    cat("   [Aviso] No se encontró el MDT (", paste0(nombre_base, ".tif"), ") en D:/H50_0344. Saltando...\n")
    next
  }
  
  # Leer el archivo excluyendo clases de ruido (6, 7 y 12)
  las <- readLAS(archivo, filter = "-drop_class 6 7 12")
  
  if (is.empty(las)) {
     cat("   [Aviso] El archivo no tiene puntos tras el filtro. Saltando...\n")
     next
  }
  
  # Normalización al vuelo usando tu MDT
  mdt <- rast(archivo_mdt)
  las_norm <- normalize_height(las, mdt)
  
  # Alineación perfecta (Bounding Box de Castilla y Leó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)
  
  # APLICAR NUESTRA FUNCIÓN FORESTAL
  metricas <- pixel_metrics(las_norm, 
                            ~mis_metricas_forestales(X, Y, Z, Intensity, ReturnNumber, Classification), 
                            res = plantilla_10x10)
  
  # Guardar resultados en D:/H50_0344/Esta
  archivo_salida <- file.path(ruta_salida, paste0(nombre_base, "_Estadisticas.tif"))
  writeRaster(metricas, filename = archivo_salida, overwrite = TRUE)
  cat("   [OK] Archivo guardado correctamente.\n")
  
  # Liberar memoria RAM
  rm(las, mdt, las_norm, plantilla_10x10, metricas)
  gc()
}

cat("\n¡PROCESO COMPLETADO! Revisa la carpeta D:/H50_0344/Esta.\n")