Initial commit with the version 1.0 for Udine
This commit is contained in:
Massimiliano Raciti
commit
44cad2cba2
103
generate_baseline.py
Normal file
103
generate_baseline.py
Normal file
@ -0,0 +1,103 @@
|
||||
#!/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.")
|
||||
199
rsun.py
Normal file
199
rsun.py
Normal file
@ -0,0 +1,199 @@
|
||||
#!/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)
|
||||
Loading…
x
Reference in New Issue
Block a user