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


# 1. PERCORSI
input_dir = r"/home/massimiliano/Solar4CE/output_dsm"
horizon_dir = r"/home/massimiliano/Solar4CE/horizon"
buildings = r"/home/massimiliano/Solar4CE/input/ELEMENTOCOPERTURA.shp"


dsm_path = os.path.join(input_dir, "dsm_clean.tif")
slope_path = os.path.join(input_dir, "slope.tif")
aspect_path = os.path.join(input_dir, "aspect.tif")

# Setup Ambiente GRASS
os.environ["GRASS_NUM_THREADS"] = "4" # Imposta il numero di core della tua CPU 
gisdb = r"/home/massimiliano/grassdata"
gisbase = r"/usr/lib/grass84"
os.environ['GISBASE'] = gisbase
sys.path.append(os.path.join(gisbase, "etc", "python"))
import grass.script as gs
import grass.script.setup as gsetup


# 2. INIZIALIZZAZIONE
gsetup.init(os.path.join(gisdb, "epsg6708", "PERMANENT"))
gs.run_command("g.mapset", flags="c", mapset="CALCOLO_CATASTO_UDINE", project="epsg6708", dbase=gisdb)

# 3. IMPORTA DSM, SLOPE, ASPECT
print("--- Passo 1: Importazione dati di base ---")
gs.run_command("r.in.gdal", input=dsm_path, output="dsm_interno", flags="o", overwrite=True)
gs.run_command("g.region", raster="dsm_interno")

gs.run_command("r.in.gdal", input=slope_path, output="slope", flags="o", overwrite=True)
gs.run_command("r.in.gdal", input=aspect_path, output="aspect", flags="o", overwrite=True)

# importa
#gs.run_command("v.in.ogr", input=buildings, output="buildings", flags="o", overwrite=True)

#Trasforma il tuo vettoriale degli edifici in un raster. Ogni pixel che cade dentro un poligono degli edifici avrà un valore, il resto sarà NULL
#gs.run_command("v.to.rast", input="buildings", output="build_mask", use="val", overwrite=True)

# Rimuove la maschera se esiste (non dà errore se non c'è)
#gs.run_command("r.mask", flags="r")
#Attiva la maschera. Da questo momento in poi, tutte le operazioni (incluso r.sun) verranno eseguite solo ed esclusivamente all'interno dei perimetri degli edifici. Il resto della mappa verrà ignorato nei calcoli.
#gs.run_command("r.mask", raster="build_mask",  overwrite=True)

gs.run_command("g.region", raster="dsm_interno")

# 4. Generazione horizon maps
prefisso_orizzonti = "horizon"
horizon_step = 60
gs.run_command("r.horizon",
               elevation="dsm_interno",
               step=horizon_step,
               output=prefisso_orizzonti,
               overwrite=True)


# 5. CALCOLO SOLARE ANNUALE 
days = [i+1 for i in range(365)]

rad_maps = []
gs.run_command("g.region", raster="dsm_interno")


#inverno 
gs.run_command("r.mapcalc", expression=f"winter_cbh = 0.40", overwrite=True, quiet=True)
gs.run_command("r.mapcalc", expression=f"winter_cdh = 0.35", overwrite=True, quiet=True)

#primavera
gs.run_command("r.mapcalc", expression=f"spring_cbh = 0.65", overwrite=True, quiet=True)
gs.run_command("r.mapcalc", expression=f"spring_cdh = 0.55", overwrite=True, quiet=True)

#estate 
gs.run_command("r.mapcalc", expression=f"summer_cbh = 0.85", overwrite=True, quiet=True)
gs.run_command("r.mapcalc", expression=f"summer_cdh = 0.75", overwrite=True, quiet=True)


#autunno 
gs.run_command("r.mapcalc", expression=f"fall_cbh = 0.75", overwrite=True, quiet=True)
gs.run_command("r.mapcalc", expression=f"fall_cdh = 0.60", overwrite=True, quiet=True)




print("\n--- Passo 3: Avvio calcolo solare (r.sun) ---")
for d in days:
    out = f"rad_{d}"
    print(f"   Elaborazione giorno {d}...", end=" ", flush=True)
    try:
        c_bh = "summer_cbh"
        c_dh = "summer_cdh"
        linke=3
        if d <= 59 or d > 334:
            c_bh = "winter_cbh"
            c_dh = "winter_cdh"
            linke=2.2
        elif 60 <= d <= 151:
            c_bh = "spring_cbh"
            c_dh = "spring_cdh"
            linke=3
        elif 152 <= d <= 243:
            c_bh = "summer_cbh"
            c_dh = "summer_cdh"
            linke=3.8
        else:
            c_bh = "fall_cbh"
            c_dh = "fall_cdh"
            linke=2.8

        gs.run_command("r.sun",
                       elevation="dsm_interno",
                       slope="slope",
                       aspect="aspect",
                       horizon_basename=prefisso_orizzonti,
                       horizon_step=horizon_step,
                       albedo_value=0.2, 
                       linke_value=linke,
                       coeff_bh=c_bh,
                       coeff_dh=c_dh,
                       glob_rad=out,
                       day=d,
                       overwrite=True)
        rad_maps.append(out)
        print("OK")
    except Exception as e:
        print(f"FALLITO: {e}")

# 6. ELABORAZIONE RASTER FINALI
print("\n--- Passo 4: Generazione output annuali kWh ---")
expr = " + ".join(rad_maps)
gs.mapcalc(f"irradiation = {expr}", overwrite=True)
# Calcolo irraggiamento potenziale Kwh/m2/anno
gs.mapcalc(f"potenziale_annuale = irradiation / 1000.0", overwrite=True)

# 7. SEGMENTAZIONE E VETTORIALIZZAZIONE
print("\n--- Passo 5: Segmentazione e Smoothing ---")
gs.run_command("i.group", group="geometria_tetti", input="slope,aspect", overwrite=True)

gs.run_command("i.segment", group="geometria_tetti", output="segments", threshold=0.05, minsize=20, overwrite=True)

gs.run_command("r.to.vect", input="segments", output="falde", type="area", overwrite=True)

# Smoothing per tetti più belli (metodo Chaiken)
#gs.run_command("v.generalize", input="falde", output="falde_smooth", method="chaiken", threshold=2, overwrite=True)

# 8. STATISTICHE E CALCOLI FV
print("\n--- Passo 6: Calcolo statistiche e Catasto FV ---")
# Zonal Stats (i_average)
gs.run_command("v.rast.stats", map="falde", raster="irradiation", column_prefix="i", method="average")
gs.run_command("v.db.addcolumn", map="falde", columns="area_m2 double precision, kwp double precision, kwh_a double precision, irr_kwh double precision, slope_deg double precision, azimut_deg double precision")

gs.run_command("v.to.db", map="falde", option="area", columns="area_m2", overwrite=True)

# Calcolo kWp: Area * 70% copertura falda * 14% efficienza (standard attuali va messo .20 o .22)
gs.run_command("v.db.update", map="falde", column="kwp", query_column="area_m2 * 0.7 * 0.22")

# Calcolo kWh_a: (Radiazione_annua_wh/m2/anno / 1000) * Area * 70% Percentuale copertura falda * 85% perdite sistema
gs.run_command("v.db.update", map="falde", column="kwh_a",
               query_column="(i_average / 1000) * area_m2 * 0.7 * 0.85 * 0.22")


gs.run_command("v.db.update", map="falde", column="irr_kwh", query_column="i_average/1000")


# aggiungo slope
gs.run_command("v.rast.stats", map="falde", raster="slope", column_prefix="s", method="average")
gs.run_command("v.db.update", map="falde", column="slope_deg", query_column="s_average")

# aggiungo azimut medio corretto (0°=Nord, senso orario)
gs.mapcalc("azimut_standard = if(aspect == 0, 0, (450 - aspect) % 360)", overwrite=True)
gs.run_command("v.rast.stats", map="falde", raster="azimut_standard", column_prefix="a", method="average")
gs.run_command("v.db.update", map="falde", column="azimut_deg", query_column="a_average")

# aggiungo risparmio e incentivi
#gs.run_command("v.db.update", map="falde", column="", query_column="i_average/1000")


# 9. EXPORT
print("\n--- Passo 7: Esportazione finale ---")
# GPKG
catasto_gpkg = os.path.join(input_dir, "catasto_solare_udine.gpkg")
gs.run_command("v.out.ogr", input="falde", output=catasto_gpkg, format="GPKG", overwrite=True)

#SHP
catasto_shp_path = os.path.join(input_dir, "catasto_solare_shapefile")
if not os.path.exists(catasto_shp_path): os.makedirs(catasto_shp_path)
gs.run_command("v.out.ogr", input="falde", output=os.path.join(catasto_shp_path, "catasto_udine.shp"), format="ESRI_Shapefile", overwrite=True)

# Raster
#final_tif = os.path.join(input_dir, "Mappa_Potenziale_Solare_Udine.tif")
#gs.run_command("r.out.gdal", input="potenziale_annuale", output=final_tif, format="GTiff", overwrite=True)

print("\n" + "="*40)
print("LAVORO COMPLETATO CON SUCCESSO!")
print("="*40)