104 lines
3.8 KiB
Python
104 lines
3.8 KiB
Python
#!/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.")
|