Calculating the intersections between planes and DEM: a Python implementation
Lo scopo di questo programma è consentire il calcolo in maniera interattiva dell'intersezione tra una superficie planare, con giacitura espressa secondo la convenzione geologica, ed un DEM. Esso può essere usato come ausilio didattico o per investigazioni sulle possibili giaciture di superfici geologiche. La base teorica dell'algoritmo è stata descritta nel precedente post: Intersezioni tra DEM e superfici planari, un tema di interesse in geologia .
Implementazione
Per il suo utilizzo è richiesta una preinstallazione di Python, con i seguenti moduli aggiuntivi (disponibili in ambienti come Python(x,y) ):
- gdal, ogr
- Numpy
- PyQt4
- Matplotlib
Una parte del codice, e precisamente quella per la gestione della finestra di visualizzazione di Matplotlib, si basa su codice dal libro di Sandro Tosi: Matplotlib for Python developers.
Questa versione, 0.1.0, è ancora allo stadio beta e necessita di ulteriori miglioramenti e testing, soprattutto su piattaforme diverse da Windows Vista, dove è stata creata e testata basandosi su Python(x,y). Ogni segnalazione di bug o miglioramenti suggeriti è utile (alberti.m65@gmail.com).
Download del codice dalla pagina gSurf: un programma Python per la determinazione dell'intersezione tra piano e DEM.
Come utilizzare il programma?
Il suo utilizzo, a parte l'installazione di Python con i moduli richiesti, non richiede altro che la decompressione della cartella in locale. Il programma viene lanciato cliccando sul file gSurf.pyw, col che si apre la finestra simboleggiata nella figura sottostante.
 |
| Screenshot dell'applicazione. |
Nella determinazione delle intersezione i parametri fondamentali, oltre al DEM che rappresenta la superficie topografica, sono la posizione di un “punto di origine” del piano, che giace nel piano che andiamo a creare, e la giacitura dello stesso piano, espressa secondo la convenzione geologica di immersione (dip direction) ed inclinazione (dip angle).
Come aiuto per scegliere una orientazione ottimale del piano, può essere visualizzato anche un livello lineare (come shapefile) che rappresenta tracce di faglie.
La sequenza suggerita di utilizzo del programma è:
- dal menu File, scegliere il DEM (shortcut: ctrl+D) ed eventualmente lo shapefile delle tracce di faglie (shortcut: ctrl+P); visualizzarli in mappa (show: on );
- definire il punto origine del piano. Questo può essere fatto sia numericamente sia graficamente, definendo la posizione di un punto in mappa con lo strumento simbolizzato dal cerchio rosso. Può essere comodo costringerne l'elevazione a quella corrispondente del DEM (Fix Z to DEM: on);
- scegliere l'orientazione preferita del piano, settando inizialmente i valori di immersione ed inclinazione con gli slider, e poi regolandoli nel dettaglio dai valori numerici negli spinbox; il bottone 'Refresh map' ricalcola la soluzione e può essere utile dopo avere modificato numericamente la posizione del punto origine;
- una volta ottenuta una soluzione interessante, salvare una immagine (icona Save nella toolbar di Matplotlib, basso al centro) e/o il livello puntuale (formato shapefile, shortcut: ctrl+S) per l'uso in report o in GIS. Nello shapefile in output sono assegnati ad ogni punto la giacitura e la posizione del punto sorgente, per gestire l'unione di più segmenti a differente giacitura nei GIS.
 |
| Tabella degli attributi dello shapefile puntuale delle intersezioni. Visualizzata con QGis. |
Oltre agli strumenti creati appositamente, sono disponibili anche quelli standard della toolbar di navigazione di Matplotlib, per zoomare, spostarsi nella mappa (pan) e per salvare un'immagine della mappa in vari fomati (png, pdf, etc.). Questi strumenti si trovano in basso al centro della finestra principale.
Problematiche note, ancora da risolvere:
- DEM particolarmente estesi possono esaurire la memoria a disposizione del programma
- lo strumento zoom della toolbar di navigazione della mappa può interferire con lo strumento di posizionamento manuale del punto origine (simbolo: cerchio rosso, in basso a destra). Talora sono stati osservati casi in cui quest'ultimo strumento non si disattiva correttamente. Può essere utile chiudere e rilanciare il programma.
- i punti salvati non sono ordinati in una sequenza adatta per creare un livello lineare. Un riordino, tale da consentire una agevole trasformazione in geometria lineare, dovrebbe essere implementato in una successiva edizione.
Una verifica con un metodo utilizzato in geologia
Il risultato consiste in una serie di punti che rappresentano l'intersezione del piano specifico col DEM, come nell'esempio in figura sottostante. In geologia esiste una tecnica manuale per ricavare da una traccia la sua orientazione ed inclinazione con le isoipse (curve di livello) di una superficie topografica. Proviamo ad utilizzarla per verificare la correttezza del risultato del programma.
Come dataset di test, utilizzo un DEM ASTER a 30 m di un settore umbro corrispondente alla zona tra la piana di Foligno a ovest e Cerreto di Spoleto – Monte Coscerno ad est. Questo settore è attaversato dalla linea della Valnerina, una serie di faglie mioceniche dell'orogenesi appenninica, a loro volta tagliate da faglie normali plio-pleistoceniche. Usaremo poi dati geologici per stimare l'orientazione di alcune di queste faglie.
Il piano teorico utilizzato per verificare la correttezza dell'algoritmo ha immersione ed inclinazione (223°, 16°). Un basso angolo di inclinazione facilita la verifica della correttezza tramite l'intersezione con le curve di livello. La posizione di origine è localizzata in un punto a quota elevata, a Nord di M.te Galenne.
 |
| Test di determinazione delle intersezioni per un piano a bassa inclinazione. Il punto di origine del piano è rappresentato dal cerchio rosso in mappa, e l'intersezione derivante dai punti bianchi (che appaiono come una linea in questa immagine). |
Una volta calcolati e salvati i risultati dell'intersezione, questi sono stati importati in GIS assieme al DEM ASTER. Le curve di livello con una spaziatura verticale di 50 m sono state determinate dal DEM e le intersezioni tra le isoipse a quota massima (1150 m slm) e minima (300 m slm) sono state utilizzate per tracciare segmenti che rappresentano la direzione della superficie planare calcolata dal programma. Questi due segmenti sono fra loro paralleli, come previsto per un piano.
Il dislivello di quota tra i segmenti delle quote 1150 e 300 m slm è di 850 m, a fronte di una distanza orizzontale tra i due segmenti di circa 2962 m. Questo corrisponde ad una inclinazione di 16.01°, in ottimo accordo col valore teorico di 16°. Anche l'orientazione stimata dell'immersione, pari a circa 222°, è in accordo più che accettabile con il valore teorico di 223°.
Si può quindi essere fiduciosi sulla qualità dei risultati prodotti dal programma.
 |
| Verifica dell'inclinazione della superficie espressa dall'intersezione, tramite la distanza da segmenti congiungenti i punti di intersezione a eguale quota, massima (1150 m slm) e minima (300 m slm). Visualizzazione con QGis. |
Quale è l'inclinazione locale della Linea della Valnerina?
La linea della Valnerina è un lineamento tettonico regionale, composto da vari segmenti di faglia ad orientazione variabile, almeno in mappa, e in vari casi con vari splay (segmenti di faglia fra loro ravvicinati e localmente subparalleli). E' generalmente considerata un lineamento con movimenti inversi-destri attiva durante l'orogenesi appenninica mio-pliocenica. La sua orientazione media nella zona della Valnerina è circa NNE-SSW ed immersione a W-NW. Faglie normali, a orientazione NW-SE ed immersione SE, taglierebbero questo lineamento, dissezionandolo in pù elementi.
I dati di campagna qui usati derivano dalla mia tesi di dottorato.
Come possiamo interpretare il suo andamento? Come successione di più tratti a differente orientazione, oppure come una superficie con una giacitura semplice, circa planare?
Un test, con punto sorgente posizionato lungo un splay superiore (occidentale) della faglia, indicherebbe che a grande scala la giacitura della Linea della Valnerina, nella zona tra Scheggino e Cerreto di Spoleto, possa essere approssimato da una superficie subplanare ad angolo molto basso (6°) e con immersione verso Ovest (277°). Avrebbe cioè le caratteristiche tipiche dei thrust appenninici a basso angolo di inclinazione e direzione N-S. Se persistesse verso Est, e con questo stesso basso angolo, andrebbe a determinare la presenza di “klippen” del suo tetto (hangingwall), nella zona del M.te Coscerno, fenomeno che non è descritto in letteratura.
Quindi è presumibile che l'inclinazione aumenti verso Est, oppure che la faglia stessa si annulli verso Est nell'arco di pochi km.
 |
| Test di determinazione del piano che approssima a scala complessiva l'orientazione dello splay superiore (occidentale) della Linea della Valnerina nel settore analizzato. L'intersezione teorica è rappresentata dalla "linea " bianca (in realtà una successione di punti. |
Nel settore settentrionale l'accordo tra andamento dell'intersezione teorica e traccia mappata dello splay superiore non è particolarmente buono, il che fa ritenere che esso devi, per direzione e/o inclinazione, dalla giacitura media descritta prima.
Posizionando un punto sorgente lungo lo splay inferiore (orientale), si vede che localmente si ottiene un ottimo fit per una giacitura (270°, 14°), cioè sempre direzione N-S, immersione a W ed angolo basso, anche se circa doppio del precedente trovato per lo splay superiore con accordo generale.
 |
| Soluzione ottimale per un piano passante nella porzione settentrionale della Linea della Valnerina, splay inferiore (orientale). L'intersezione teorica è rappresentata dalla "linea " bianca . |
Nella sua porzione meridionale l'andamento dell'intersezione devia però marcatamente da quello mappato, e andrebbe ad interessare il fianco opposto della Valnerina, dal lato di M.te Coscerno. Appare probabile quindi che in questo settore ci siano variazioni rispetto alla precedente giacitura. In questo settore meridionale l'inclinazione sembra meno facilmente caratterizzabile: valori compatibili vanno da 80°-90° a 10°. Il migliore accordo locale è ottenuto con una giacitura (315°, 37°).
 |
| Soluzione per il piano passante per lo splay inferiore (occidentale) della Linea della Valnerina, nel settore meridionale. L'intersezione teorica è rappresentata dai punti bianchi. |
Riepilogando, visto nel complesso della zona considerata, la Linea della Valnerina sarebbe approssimabile con un piano a direzione N-S a bassissima inclinazione verso W (5-10°). Visto con maggiore dettaglio, sono ipotizzabili segmenti a orientazione lievemente differente, nella zona settentrionale con direzione N-S e inclinazioni sui 15°, in quella meridionale a direzione NE-SE e inclinazione forse maggiore, di circa 35°-40°.
Per creare un modello 3D di medio dettaglio della faglia nell'intorno della sua traccia topografica, occorrerebbe quindi collegare due superfici subplanari con queste differenti orientazioni.
Questo è possibile con software di modellazione 3D, soprattutto se orientati alla geologia (p.e. GoCAD, 3DMove).
Un obiettivo sarebbe quello di renderlo possibile in futuri rilasci di questa stessa applicazione, anche se solo per elaborazioni semplici.
4 comments:
Complimenti per il tool, come geologo apprezzo molto... ;)
grazie Pietro,
l'intento è proprio della utility per la geologia, e spero che col tempo possa migliorarlo ed arricchirlo.
In qualità di Geologo, mi unisco anche io ai complimenti di Pietro!
Di nuovo grazie. Il prossimo passo, come anche suggeritomi da un collega, sarebbe quello di crearne una versione per Quantum GIS. Forse tra una quindicina di giorni.
Post a Comment