QGIS: confine vettoriale 3D di un DTM

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:

  1. 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;

Immagine 5
Modello Digitale del Terreno

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.

Immagine 7
calcolatore raster: “valaosta@1″=0
Immagine 8
screenshot QGIS

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.

offset
processing -> geo-algoritmo GRASS

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:

PUNTI_PROFILO
processing -> geo-algoritmo SAGA

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.

Immagine 14
DTM – linea di offset_40 – punti (pixel centroid)
pti_prof
screenshot canvas QGIS

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

Immagine 13
finestra di dialogo

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.

Immagine 20

OSGeo4W Shell

Immagine 17
database SpatiaLite

oppure  usando il calcolatore di campi di QGIS e la funzione:

geom_to_wkt( $geometry ) 

Immagine 19
rappresentazione WKT

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:

Immagine 21
finestra di dialogo salva con nome…

aprendo il file con LibreOffice Calc vedremo:

Immagine 23
file – iniziale

dobbiamo trasformarla utilizzando ‘trova e sostituisci’ ed eliminare tutte le altre colonne tranne la A, eliminare la riga 1:

Immagine 24
file – trasformato

per ultimare le modifiche occorre aprire il file con NotePad++ (o similari):

Immagine 25
file -iniziale

apportare le modifiche (tramite trova e sostituisci) per ottenere:

Immagine 27
NotePad++ – finale

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

osgeo4w
OSGeo4W
wwww
OSGeo4W

per verificare che la linea generata è 3D:

Immagine 29
OSGeo4W Shell
postgis
PostGIS da QGIS

Unico modo per calcolare la lunghezza della linea in 3D è attraverso PostGIS:

lungh_postgis
PostGIS da QGIS

lunghezza linea 2D:

lunghezza_2D
SpatiaLite in QGIS
  • 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

  1. avviamo spatialite_gui;
  2. creare un nuovo SQLite DB – my_db;
  3. importare lo shapefile paletti_3D.shp nel database my_db;
  4. 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;
  5. esportare il risultato in shapefile – 3D_linea;
  6. importarlo nel my_db per verificare la dimensione con ‘check geometries’.
sql_3
query
sql_5
verifica 3D

 

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.

Immagine 34
SAGA – visualizzazione 2D/3D

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:

Procedure e geo-algoritmi:

  1. Riclassificazione DTM con calcolatore raster;
  2. poligonizzazione DTM appiattito;
  3. conversione da poligono a linea;
  4. v.parallel -> geo-processing GRASS;
  5. profil from lines -> geo-processing SAGA;
  6. v.drope -> geo-processing GRASS;
  7. salva file in CSV 3D con geometria WKT;
  8. LibreOffice e NotePad++;
  9. creazione file virtuale .vrt;
  10. OSGeo4W -> creazione shapefile 3D e verifica;
  11. ST_3DLength (geom) per calcolo lunghezza linea 3D – PostGIS;
  12. ST_Length (geom) per calcolo lunghezza linea 2D – SpatiaLite;
  13. Blank Flowchart - New Page (1)
    flowchart

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

 

by-nc-sa.eu_


 

Dal gruppo pubblico GIS ITALIA di Facebook:

Immagine 4
screenshot Facebook
Immagine 3
screenshot Facebook
Annunci

Rispondi

Inserisci i tuoi dati qui sotto o clicca su un'icona per effettuare l'accesso:

Logo WordPress.com

Stai commentando usando il tuo account WordPress.com. Chiudi sessione / Modifica )

Foto Twitter

Stai commentando usando il tuo account Twitter. Chiudi sessione / Modifica )

Foto di Facebook

Stai commentando usando il tuo account Facebook. Chiudi sessione / Modifica )

Google+ photo

Stai commentando usando il tuo account Google+. Chiudi sessione / Modifica )

Connessione a %s...