In questo articolo cercherò di rispondere al seguente problema:
Un amico ha un esteso appezzamento di terreno in alta montagna e vuole realizzare una recinzione per proteggerlo da eventuali intrusioni non autorizzate; durante un aperitivo mi chiede come calcolare la lunghezza di filo spinato e il numero di paletti che deve acquistare senza dover spendere un patrimonio.
Questo è un tipico problema da risolvere con approccio GIS 3D o meglio 2.5D, in questo articolo si utilizzeranno solo software Open Source.
Partiamo da una ipotesi:
- abbiamo il DTM perfettamente tagliato per l’appezzamento di terreno.
Procedimento
- importiamo il DTM in QGIS;
- determiniamo il confine del DTM in proiezione (in 2D);
- determiniamo i punti (in 3D) dei paletti;
- generiamo la linea 3D che congiunge tutti i paletti.
Importiamo il DTM in QGIS;
determiniamo il confine del DTM in proiezione 2D:
la tecnica utilizzata è quella della poligonizzazione del raster dopo averlo riclassificato con l’ausilio del calcolatore raster di QGIS (vedi altro articolo su Linea di costa); dopo aver appiattito il DTM si passa alla poligonizzazione che darà come output un vettore poligonale 2D con limiti, appunto, il DTM. Successivamente si trasforma il poligono in linea ottenendo il confine 2D dell’area di impronta del DTM.
Al vettore linea 2D verrà applicato un offset interno 40 m in modo tale che la linea ricada all’interno dei pixel del DTM: questo è importante e si capirà meglio dopo.
determiniamo i punti (in 3D) dei paletti:
Utilizziamo la linea di offset_40, appena creata, per generare i punti (pixel centroid) lungo la linea attraverso l’uso del processing:
questo geo-algoritmo genera dei punti 2D (pixel centroid) seguendo la linea che interseca il raster (offset_40), i pixel del DTM hanno dimensione 80×80 quindi la distanza tra i punti sarà di 80 m. Se volessimo una distanza minore tra i punti occorrerebbe utilizzare un DTM con pixel più piccoli.
Nella figura di sopra sono evidenziati il primo e l’ultimo punto, nella tabella degli attributi sono presenti vari attributi DIST: che è la distanza in 2D; DIST SURF: che è la distanza in 3D e poi la tripletta di coordinate. Si noterà che le due distanze (2D e 3D) sono diverse.
A questo punto siamo pronti per trasformare i punti 2D in 3D utilizzando il DTM e un comando GRASS che ritroviamo in processing: v.drape
per verificare che i punti sono effettivamente 3D si possono seguire varie strade, per esempio scrivere ogrinfo pippo.shp in Osgeo4W Shell; oppure importare lo shapefile in un database spatialite.
OSGeo4W Shell
oppure usando il calcolatore di campi di QGIS e la funzione:
geom_to_wkt( $geometry )
generiamo la linea 3D che congiunge tutti i paletti:
le triplette di coordinate le ricomponiamo in una linea 3D in WKT, concatenando le triplette, che in un CSV spaziale potrebbe essere così (paletti_3d.csv):
WKT,id
“LINESTRING Z (345033.07756361004 4221015.5546970163 50,330734.94084110629 4207332.313714385 30,331946.96878985036 4199477.8127947021 50,344424.26383262046 4197025.7027728558 20)”,1
per iniziare a realizzare il file paletti_3d.csv basti esportare il file paletti_3D.shp in formato CSV, sempre da QGIS:
aprendo il file con LibreOffice Calc vedremo:
dobbiamo trasformarla utilizzando ‘trova e sostituisci’ ed eliminare tutte le altre colonne tranne la A, eliminare la riga 1:
per ultimare le modifiche occorre aprire il file con NotePad++ (o similari):
apportare le modifiche (tramite trova e sostituisci) per ottenere:
finalmente abbiamo creato il primo file paletti_3d.
Creare un file virtuale paletti_3D.vrt cosi fatto:
<OGRVRTDataSource>
<OGRVRTLayer name=”paletti_3d“>
<SrcDataSource>paletti_3d.csv</SrcDataSource>
<GeometryType>wkbLineString25D</GeometryType>
<LayerSRS>PROJCS[inserire codice OGR WKT http://spatialreference.org del SRS]</LayerSRS>
<GeometryField encoding=”WKT” field=’WKT’ > </GeometryField >
</OGRVRTLayer>
</OGRVRTDataSource>
infine creare lo shapefile 3D cosi:
ogr2ogr -overwrite linea_3D.shp paletti_3D.vrt
o meglio con tutto il percorso:
ogr2ogr -overwrite C:\Users\Salvatore\Desktop\linea_3D\work\linea_3D.shp C:\Users\Salvatore\Desktop\linea_3D\work\paletti_3d.vrt
per verificare che la linea generata è 3D:
Unico modo per calcolare la lunghezza della linea in 3D è attraverso PostGIS:
lunghezza linea 2D:
- La lunghezza della linea 2D è 104,145 km
- La lunghezza della linea 3D è 117,827 km
Vi è una differenza di circa 14 km cioè un errore di circa l’11% sulla lunghezza reale.
Vediamo come ottenere la linea 3D utilizzando SpatiaLite
SpatiaLite
- avviamo spatialite_gui;
- creare un nuovo SQLite DB – my_db;
- importare lo shapefile paletti_3D.shp nel database my_db;
- eseguire la seguente query:
SELECT LINE_ID, MakeLine(MakePointZ(CAST(X AS float),CAST(Y AS float),CAST(Z AS float),32632)) AS geom FROM paletti_3D GROUP BY LINE_ID;
- esportare il risultato in shapefile – 3D_linea;
- importarlo nel my_db per verificare la dimensione con ‘check geometries’.
Concludendo:
Questo risultato dipende molto dalla dimensione dei pixel del DTM e dal numero di punti generati per ottenere la linea in 3D. Ho fatto un’altra prova con distanza dei punti molto minore e il risultato è 129 km.
Note finali: L’articolo potrebbe sembrare complicato, difficile da seguire ma non è cosi. Per ottenere il risultato finale occorre circa mezz’oretta o meno. Tengo molto a sottolineare che è stato sviluppato tutto con software Open Source (cioè senza spese per acquisto di licenze software); con sw proprietari (acquistare la licenza costo migliaia di euro) ci si impiegherebbe molto meno ma non si capirebbe il modo con cui arriva alla soluzione! I sw OS sono straordinari dal punto di vista didattico, questo esempio è un caso lampante di teoria e tecnica allo stato puro.
Ringraziamenti:
Ringrazio il collega Ing. Antonio Vinci per l’input sulla problematica.
Ringrazio l’amico Andrea Borruso per avermi guidato sul formato WKT , CSV 3D e SpatiaLite.
Software Open Source utilizzati:
- QGIS Lyon 2.12.2;
- LibreOffice 5;
- NotePad++;
- OSGeo4W shell;
- SpatiaLite; SpatiaLite_gui;
- PostgreSQL/PostGIS;
- SAGA.
Procedure e geo-algoritmi:
- Riclassificazione DTM con calcolatore raster;
- poligonizzazione DTM appiattito;
- conversione da poligono a linea;
- v.parallel -> geo-processing GRASS;
- profil from lines -> geo-processing SAGA;
- v.drope -> geo-processing GRASS;
- salva file in CSV 3D con geometria WKT;
- LibreOffice e NotePad++;
- creazione file virtuale .vrt;
- OSGeo4W -> creazione shapefile 3D e verifica;
- ST_3DLength (geom) per calcolo lunghezza linea 3D – PostGIS;
- ST_Length (geom) per calcolo lunghezza linea 2D – SpatiaLite;
Dati e progetto: nella cartella troverete tutti i file più un modello Sextante.
Il Modello Digitale del Terreno è stato scaricato dal PCN tramite il servizio di download WCS in QGIS.
demo spatialite
Riferimenti:
OGC WKT:
http://spatialreference.org/ref/epsg/wgs-84-utm-zone-32n/ogcwkt/ http://spatialreference.org/ref/epsg/wgs-84-utm-zone-32n/
Virtual Format
http://www.gdal.org/drv_vrt.html
wikipedia
https://it.wikipedia.org/wiki/Well-Known_Text#Oggetti_geometrici
Dal gruppo pubblico GIS ITALIA di Facebook:
pigrecoinfinito 