Skip navigation.
Home
Istituto Scienze della Terra - SUPSI

Altezze edifici con Grass

Contesto generale

Per la progettazione all’interno degli spazi urbani, in generale da parte degli architetti e pianificatori, è sempre più importante prendere in considerazione aspetti legati ai conflitti tra edifici ed elementi circostanti (viabilità, parcheggi, alberi, giardini, arredo urbano, spazi pubblici o quant’altro). In particolare modelli e dati tridimensionali permettono una miglior elaborazione di progetti per nuove costruzioni e/o riqualificazione di edifici. Indici volumetrici, diritti di vista, ombreggiature, questioni legate al rumore e tutto ciò che può portare ad un’utilizzazione razionale degli spazi necessita la conoscenza di un’informazione preziosa che è l’altezza degli edifici.

La grande maggioranza delle nostre città non dispone ancora di un Catasto 3D e rilevare edificio per edificio in modo puntuale sarebbe un’impresa alquanto dispendiosa in termini di fatica e di tempo.

Per questo motivo, avendo a disposizione il modello di terreno (DTM) ed il modello di superficie (DOM), vi proponiamo una procedura semplice e veloce per acquisire la quota di edifici di una zona urbana relativamente ampia. Il tutto attraverso un esempio pratico e con l’utilizzo del GIS OpenSource Grass 6.4 .
 

Immagine estratto di Berlino in 3D

 

Dati input

Per effettuare il lavoro si hanno a disposizione :
 

  • DTM : il modello numerico del terreno con risoluzione 2m (fornitoci in formato raster Esri Grid poi riconvertito/esportato in ascii Esri)
  • DOM : il modello di superficie con risoluzione 2m (fornitoci non grigliato in formato shape file 3D Esri come nuvola di punti. Punti comunque già filtrati (outliers, algoritmi, …) quindi approvati ufficialmente)
  • Edifici del catasto ufficiale (formato Esri shape file)

I dati citati sono di gestione e competenza dell’Ufficio cantonale TI “Sezione delle bonifiche e del catasto”.

 

Principio Lidar per l’acquisizione DTM e DOM

 

 

Tappe di lavoro

 

Tappa 1) Importare in Grass 6.4 il file ascii (formato in codifica Esri) del DTM 2 m di nostro interesse. Comando da utilizzare r.in.arc

Tappa 2) Con g.region settare la regione basandosi sul raster appena importato (risoluzione sia 2D che 3D uguale a 2 m).

Tappa 3) Importare i punti DOM con il comando v.in.ogr

Visto che i punti input sono degli shape 3D non ho quindi informazioni di quota come attributo e bisognerà con v.in.ogr attivare il flag -z che consente di importare anche la geometria tridimensionale.

Tappa 4) L’idea è quella di generare una superficie a partire dai punti DOM. Generalmente per creare un modello continuo si passa attraverso il calcolo di un TIN e successiva rasterizzazione (in ArcGIS per esempio 1. Create TIN 2. TIN to Raster). Con Grass potremmo pensare d’utilizzare dapprima il comando v.delaunay od eventualmente v.voronoi ed in seguito un v.to.rast per rasterizzare le superfici. Purtroppo però v.delaunay crea solo i triangoli (aree o linee se flag -l) ai quali ho un’informazione di quota legata al centroide ma nessun andamento tridimensionale. I triangoli sono piatti (2D) ed una rasterizzazione non sortirebbe nessun risultato.

Dobbiamo forzatamente applicare un’interpolazione a partire dai punti DOM 3D. Usiamo il comando v.surf.bspline, l’unico con l’opzione layer number = 0 (dove viene utilizzata la geometria z di quota dei punti come valore d'interpolazione per generare il raster). Otteniamo il raster del DOM con risoluzione 2 m.

 

Punti DOM 3D ai vertici e Triangolazione di Delaunay (quest’ultima non va applicata!)

 

 

Tappa 5) Importare lo shp degli edifici con il solito v.in.ogr

Attenzione con la versione Grass 6.4 può capitare ci siano dei problemi nell'importare la tabella degli attributi dello shp vettore. Grass importa quindi esclusivamente solo la geometria e bisognerà quindi ricostituire una tabella, utilizzando v.db.addtable, per implementare poi gli attributi di quota derivanti dal DTM-DOM (vedi tappe successive).

Utilizzando il solito v.in.ogr, ma trattandosi di poligoni, selezionare "Do not create attribute table" in "Attributes" senò se cerco d’importare anche gli attributi non importa niente. Come detto con questa opzione vengono però importati solo i poligoni (geometria) senza tavola attributiva.

 

DOM raster ed edifici. Si possono già notare gli edifici che presentano quote assolute più alte (colori tendenti al viola)


 

 

Tappa 6) Con v.db.addtable creare ed associare una tabella attributiva agli edifici (solo se v.in.ogr non importa in modo completo uno shape). Esempio :

v.db.addtable map=edifici@ex_Campari_Campus table=edifici layer=1 'columns=cat integer,h_edificio DOUBLE
PRECISION,h_slm DOUBLE PRECISION'

con columns= posso crere insieme alla tavola anche gli attributi che la compongono specificandone il type. In questo caso oltre all'attributo obbligatorio cat integer ne ho creati altri due di tipo double precision.

Tappa 7) Si tratta ora di calcolare la differenza tra DOM e DTM (DOM - DTM) per generare una mappa della differenza di quota rispetto al terreno. Operazione che va eseguita con un semplice comando di r.mapcalc

r.mapcalc 'Delta_DOMDTM=dom2m_vsurfbspline-dtm2_camp'

Il risultato conterrà i dati che serviranno per l'acquisizione dell'altezza edifici. Nell’immagine di dettaglio si possono distinguere interpretando la colorazione del raster gli edifici con tetti a falde o tetti piatti nonché la presenza di vegetazione (alberi, ecc.)

 

 Raster Delta DOM-DTM dopo calcolo con r.mapcalc + layer edifici

 

 

Dopo cambiamento colorazione, dettaglio DOM + edifici con interpretazione visiva dell’immagine (comunque ogni pixel informa sulla quota rispetto al terreno DTM)

 

 

 

Tappa 8) Non resta altro che acquisire i dati relativi le quote di edifici sul layer edifici utilizzando il comando v.rast.stats che calcola statistiche su vettori (es. edifici nel nostro caso) associate ad un raster. Questo comando crea tutte le statistiche (min, max, mean, ecc.) senza possibilità di scelta.

Per ogni edificio (quindi per ogni entità poligono) acquisiremo una statistica completa sui valori dei pixel tra i quali delta_max che rappresenta la quota massima (al punto più alto, colmo) dell’edificio o delta_mean (media dei valori), ecc. , ecc. Tutte le informazioni sono importanti per capire la realtà sul terreno che tuttavia andrebbe controllata se si riscontrano incongruenze o se si desiderano maggiori dettagli o precisione su alcuni edifici in particolare.

Importante il prefisso della colonna (es. delta, dopo il calcolo delta_min, delta_max, ecc.) che consente all'utilizzatore dei dati di capire da quale raster é stata fatta l'acquisizione (Delta_DOMDTM per noi).

 


Interrogazione del layer edifici dopo v.rast.stats


 

 

 

 

Tappa 9) I risultati (ovvero il vettore edifici aggiornato con i valori di quota) può essere esportato per esempio in formato Esri (v.out.ogr) per poi essere utilizzato nel generare mappe e stampe (più conviviali in ArcGIS o programmi Desktop), come pure da base per l'utilizzo in altri software CAD (es. ArchiCad) per l'elaborazione di 3D più attrattivi dal punto di vista grafico, visto oltretutto che questo genere di dati é fortemente richiesto nel settore della progettazione architettonica

 

 

 

Ortofoto con indicazione sugli edifici della quota massima registrata per ogni edificio

 

 

 

Riassunto con linee di comando in Grass
 

Tappa 1)

r.in.arc input=/home/ubuntu/share/DTM_DOM_raster_export_ascii/dtm2_camp.asc output=dtm2_camp type=FCELL mult=1.0 --overwrite

 

Tappa 2)

g.region rast=dtm2_camp@ex_Campari_Campus res=2 res3=2 --overwrite

 

Tappa 3)

v.in.ogr -z -o dsn=/home/ubuntu/share/DTM_DOM_raster_export_ascii/shp_DOM/DOM.shp output=DOM_punti3d min_area=0.0001 snap=-1


Tappa 4)

v.surf.bspline input=DOM_punti3d@ex_Campari_Campus raster=dom2m_vsurfbspline sie=4 sin=4 type=bilinear lambda_i=1 layer=0 --overwrite

 

Tappa 5)

v.in.ogr -t –o dsn=/home/ubuntu/share/DTM_DOM_raster_export_ascii/edifici_altezze/edifici_altezze.shp output=edifici min_area=0.0001 snap=-1 --overwrite

 

Tappa 6)

v.db.addtable map=edifici@ex_Campari_Campus table=edifici layer=1 'columns=cat integer,h_edificio DOUBLE PRECISION,h_slm DOUBLE PRECISION'

 

Tappa 7)

r.mapcalc 'Delta_DOMDTM=dom2m_vsurfbspline-dtm2_camp'

 

Tappa 8)

v.rast.stats -e vector=edifici@ex_Campari_Campus layer=1 raster=Delta_DOMDTM@ex_Campari_Campus colprefix=delta percentile=90

 

 

Ing. Alessio Spataro