# SCRIPT DE DIAGNÓSTICO (Modo Secuencial - 3 archivos)
suppressPackageStartupMessages({
  library(lidR)
  library(terra)
})

ruta_laz <- "P:/Cartografia/02_Altimetria/023_LIDAR/2025/Todos1"
ruta_mdt <- "T:/MDT2x2" 
ruta_salida <- "P:/Cartografia/02_Altimetria/023_LIDAR/2025/Salida1"

# Función de métricas
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))
}

archivos_lidar <- list.files(ruta_laz, pattern = "\\.laz$", full.names = TRUE)
cat("\nTotal archivos encontrados:", length(archivos_lidar), "\n")

# Tomamos solo los 3 primeros archivos para la prueba
archivos_prueba <- archivos_lidar[1:3]

for (archivo in archivos_prueba) {
  nombre_base <- tools::file_path_sans_ext(basename(archivo))
  cat("\n--------------------------------------------------\n")
  cat("-> PROCESANDO:", nombre_base, "\n")
  
  archivo_mdt <- file.path(ruta_mdt, paste0(nombre_base, ".tif"))
  archivo_salida <- file.path(ruta_salida, paste0(nombre_base, "_Estadisticas.tif"))
  
  if (!file.exists(archivo_mdt)) {
    cat("   [ERROR] No se encuentra el MDT en", archivo_mdt, "\n")
    next
  }
  
  cat("   Paso 1: Leyendo LAS (ignorando clases 6, 7 y 12)...\n")
  las <- readLAS(archivo, filter = "-drop_class 6 7 12")
  if (is.empty(las)) {
     cat("   [AVISO] El archivo LAS está vacío tras filtrar.\n")
     next
  }
  
  cat("   Paso 2: Leyendo MDT y normalizando alturas...\n")
  mdt <- rast(archivo_mdt)
  las_norm <- normalize_height(las, mdt)
  
  cat("   Paso 3: Creando malla de alineacion 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("   Paso 4: Calculando metricas estándar y personalizadas...\n")
  # Recordatorio: lidR::stdmetrics calcula las variables predeterminadas (zmax, zmean, area, etc.)[cite: 2]
  metricas <- pixel_metrics(las_norm, 
                            ~mis_metricas_forestales(X, Y, Z, Intensity, ReturnNumber, Classification), 
                            res = plantilla_10x10)
  
  cat("   Paso 5: Guardando TIF en disco de red...\n")
  writeRaster(metricas, filename = archivo_salida, overwrite = TRUE)
  
  cat("   [EXITO] Archivo finalizado correctamente.\n")
  
  rm(las, mdt, las_norm, plantilla_10x10, metricas)
  gc()
}
cat("\n[ FIN ] Prueba de diagnostico finalizada.\n")