import os
import subprocess
import glob

# Definición de rutas
directorio = r"P:\Cartografia\02_Altimetria\023_LIDAR\2025\68Bandas\LEON"
lista_archivos_txt = os.path.join(directorio, "lista_archivos.txt")
archivo_vrt = os.path.join(directorio, "LEON_mosaico.vrt")
archivo_salida = os.path.join(directorio, "LEON.IMG")

# Parámetros espaciales y de datos
epsg = "EPSG:25830"
resolucion = "10 10"
bounding_box = "165000 4655850 358760 4789070"

# --- CONFIGURACIÓN DE VALOR NULO ---
# Dependiendo de cómo se generaron las bandas, el nulo puede ser -9999, "nan", o -3.4028234663852886e+38. 
# Ajusta este valor si el software anterior utilizó otro identificador.
valor_nodata = "-9999" 

# 1. Generar la lista de archivos de entrada
print("Localizando archivos PNOA*.tif...")
archivos_tif = glob.glob(os.path.join(directorio, "PNOA*.tif"))

with open(lista_archivos_txt, "w") as f:
    for tif in archivos_tif:
        f.write(f"{tif}\n")
        
print(f"Se han listado {len(archivos_tif)} archivos.")

# 2. Construir el Mosaico Virtual (VRT)
comando_vrt = [
    "gdalbuildvrt",
    "-vrtnodata", valor_nodata,  # Rellena los huecos sin archivos con el valor nulo
    "-input_file_list", lista_archivos_txt,
    archivo_vrt
]

print("Construyendo Raster Virtual (VRT)...")
subprocess.run(comando_vrt, check=True)

# 3. Traducir el VRT a formato .IMG (HFA) 
comando_translate = [
    "gdal_translate",
    "-of", "HFA",
    "-a_srs", epsg,
    "-te", *bounding_box.split(),
    "-tr", *resolucion.split(),
    "-a_nodata", valor_nodata,   # Define formalmente el valor nulo en el archivo .IMG
    "-co", "COMPRESSED=YES",
    archivo_vrt,
    archivo_salida
]

print("Generando archivo final LEON.IMG... (Este proceso tomará tiempo)")
subprocess.run(comando_translate, check=True)

print("Proceso finalizado con éxito.")