#!/usr/bin/env python3
import os
import sys
import subprocess
import glob
import json

# --- 1. CONFIGURAZIONE PERCORSI ---
laz_folder = r"/home/massimiliano/Solar4CE/laz_tiles"
output_dir = r"/home/massimiliano/Solar4CE/output_dsm"
buildings = r"/home/massimiliano/Solar4CE/input/ELEMENTOCOPERTURA.shp"

# Configurazione GRASS
gisbase = r"/usr/lib/grass84"
gisdb = r"/home/massimiliano/grassdata"
location = "epsg6708"
mapset = "Udinemapset"

# Parametri Solari
resolution = 0.5
linke = 3.0

if not os.path.exists(output_dir):
    os.makedirs(output_dir)

# --- 2. SETUP AMBIENTE GRASS ---
os.environ['GISBASE'] = gisbase
sys.path.append(os.path.join(gisbase, "etc", "python"))

import grass.script as gs
import grass.script.setup as gsetup

# --- 3. INIZIALIZZAZIONE SESSIONE GRASS ---
location_path = os.path.join(gisdb, location)
if not os.path.exists(location_path):
    subprocess.run(["grass", "-c", "EPSG:6708", location_path, "--exec", "g.proj", "-p"], check=True)

gsetup.init(os.path.join(location_path, "PERMANENT"))

if not os.path.exists(os.path.join(location_path, mapset)):
    gs.run_command("g.mapset", flags="c", mapset=mapset, location=location, dbase=gisdb)
else:
    gs.run_command("g.mapset", mapset=mapset, location=location, dbase=gisdb)


# --- 4. FUNZIONE PDAL ---
def build_dsm(laz_files, filter_range, output_path):
    readers = [{"type": "readers.las", "filename": f} for f in laz_files]
    pipeline = {
        "pipeline": readers + [
            {"type": "filters.range", "limits": filter_range},
            {"type": "writers.gdal", "filename": output_path, "resolution": resolution, "output_type": "max"}
        ]
    }
    json_path = output_path.replace(".tif", ".json")
    with open(json_path, "w") as f:
        json.dump(pipeline, f, indent=4)
    subprocess.run(["pdal", "pipeline", json_path], check=True)
    return output_path


# --- 5. ESECUZIONE PDAL ---
laz_files = glob.glob(os.path.join(laz_folder, "*.la[zs]"))

if laz_files:
    path_roofs = os.path.join(output_dir, "dsm_roofs.tif")
    path_full = os.path.join(output_dir, "dsm_full.tif")

    print("\nGenerazione DSM tramite PDAL...")
    build_dsm(laz_files, "Classification[6:6]", path_roofs)
    build_dsm(laz_files, "Classification[2:6]", path_full)

    # --- 6. IMPORT ---
    print("\nGRASS -> Importazione raster...")
    gs.run_command("r.in.gdal", input=path_roofs, output="dsm_roofs", flags="o", overwrite=True)
    gs.run_command("r.in.gdal", input=path_full, output="dsm_full", flags="o", overwrite=True)


    # --- 7. GEOMETRIE (SLOPE & ASPECT) E SALVATAGGIO IMMEDIATO ---
    print("\nElaborazione geometrie del terreno...")
    gs.run_command("g.region", raster="dsm_full", align="dsm_full")

    gs.run_command("r.neighbors", input="dsm_roofs", output="dsm_smooth", method="average", size=3, overwrite=True)
    gs.run_command("r.neighbors", input="dsm_smooth", output="dsm_clean", method="median", size=3, overwrite=True)

    # ---ESPORTAZIONE DSM_CLEAN ---
    print("Salvataggio dsm_clean.tif nella cartella output...")
    gs.run_command("r.out.gdal", input="dsm_clean", output=os.path.join(output_dir, "dsm_clean.tif"), format="GTiff",
                   overwrite=True)

    #print("Calcolo Slope e Aspect...")
    gs.run_command("r.slope.aspect", elevation="dsm_clean", slope="slope", aspect="aspect", overwrite=True)

    # ESPORTAZIONE IMMEDIATA SLOPE E ASPECT
    print("Esportazione Slope e Aspect in formato GeoTIFF...")
    gs.run_command("r.out.gdal", input="slope", output=os.path.join(output_dir, "slope.tif"), format="GTiff",
                   overwrite=True)
    gs.run_command("r.out.gdal", input="aspect", output=os.path.join(output_dir, "aspect.tif"), format="GTiff",
                   overwrite=True)


else:
    print("Nessun file LiDAR (.laz) trovato.")