Pierre Raybaut rilascia WinPython, ambiente portabile per Python 2.7 e 3.3

Raybaut è lo sviluppatore di Python(x,y) e di packages come guiqwt, guidata, spyder e userconfig, tutti rigorosamente OS.
In questi giorni ha rilasciato una nuova versione di WinPython, ambiente alternativo a Python(x,y) in quanto è portabile, ma conserva l'orientamento per l'analisi e la gestione di dati scientifici.
Presenta infatti packages come NumPy, SciPy, Matplotlib e PyQwt, PyQt, numexpr, Pandas, scikit-image, scikit-learn, sympy, tables, xlrd e xlwt ed anche Cython.
E' possibile utilizzare IPython ed IPython Notebook, interessante per testare e divulgare codice Python.

La sua portabilità consente di avere installazioni plurime di Python che non interferiscono, p.e. è possibile mantenere Python 2.7 per la produzione ed avere a disposizione  Python 3.3 per testare la nuova serie.

Link:
https://code.google.com/p/winpython/

Plotting profiles from GPX files within Quantum GIS

A new version of qProf, a Quantum GIS plugin, allows the plotting of profiles by reading the track points stored inside a GPX file. This new tool complements the previous implementation of profile calculation (elevation and slope) from DEMs and line tracks.

Screenshot for a GPX profile produced by qProf plugin v. 0.1.4.
Upper window: elevations, lower w.: slopes.
Data (implicit) courtesy of S. Peduzzi.

Data in source GPX file should contain point tracks, with point coordinates in latitude and longitude, plus elevation and time. Elevations are generally referred to the WGS84 datum, not to the orthometric height.

The algorithm converts the latitude-longitude coordinates in ECEF Cartesian coordinates, avoiding the need for any planar, 2D projection. The 3D distances between consecutive track points are easily calculated, and with some trigonometry formulas, slope and 2D horizontal distances between points are then derived by assuming a vertical offset between consecutive points equal to the elevation differences derived from the GPX file.

The resulting profile can be saved as .csv files, 2D point shapefiles, or as 3D line shapefiles.

This plugin can be easily installed using the QGis plugin manager, or by downloading it from the QGis plugin repository and extracting it into your personal Python plugin folder.

Ricuperare giaciture geologiche da immagini satellitari e DEM: nuova versione del plugin qgSurf per QGis

Presento una nuova versione sperimentale di qgSurf, la 0.2.0, che ha come intento quello di permettere, all'interno di Quantum Gis, di utilizzare immagini satellitari, come quelle liberamente disponibili in GoogleEarth, accompagnandole a DEM, per ricavare in maniera interattiva, per trial-and-error, la giacitura di piani geologici .

Schermata del plugin, con un esempio di analisi su M.te Alpi (Basilicata), da immagine Google, e dati DEM Aster a 30 m di risoluzione.

Come detto la versione è sperimentale (ogni segnalazione di baco è benvenuta -> alberti.m65@gmail.com) ed è stata testata fondamentalmente sulla versione Windows di Quantum GIS 1.8.0. A differenza della versione precedente, l'interazione con la mappa avviene all'interno del canvas di Quantum GIS, il che permette di utilizzare come base cartografica qualsiasi livello vettoriale o raster, come appunto le immagini satellitari disponibili dai servizi di GoogleEarth, Bing, etc, caricabili in Quantum GIS attraverso il plugin OpenLayers.

Finestra del plugin qgSurf v. 0.2.0, tab 'Geographic data'.

Brevemente, la successione di utilizzo di qgSurgf v. 0.2.0 è la seguente:
0) occorre settare come sistema proiettivo del progetto QGis quello del DEM utilizzato per la determinazione;
1) dal tab 'Geographic data' si definisce il DEM sul quale viene calcolata l'intersezione del piano. Questo permette quindi di definire nella mappa il punto sorgente (source point) che contiene il piano geologico;
2) dal tab 'Geological data' si settano i parametri (immersione ed inclinazione) del piano, e col bottone 'Calculate intersection' vengono calcolati i punti di intersezione teorici tra DEM e piano geologico. Questi parametri possono essere variati sino a trovare un 'best-fit' con un lineamento (stratificazione, faglia, etc.) visibile in immagine satellitare.
3) nel tab 'Output' troviamo gli strumenti per salvare le intersezioni, come punti o linee, in shapefile.

Ovviamente se si dispone di misurazioni strutturali di terreno è possibile testare la loro corrispondenza con gli andamenti visibili in immagini satellitari/aeree o in cartografie geologiche georiferite.

Il plugin può essere scaricato direttamente all'interno di Quantum GIS, oppure da http://plugins.qgis.org/plugins/qgSurf/ per essere poi decompresso nella cartella dei plugin Python di Quantum GIS.

qProf: un nuovo plugin QuantumGIS per la creazione di profili topografici

Questo è il primo rilascio di QProf, un plugin per la creazione di  profili di elevazione e di pendenza 2D e 3D, a partire da uno o più DEM. E' possibile visualizzare più profili contemporaneamente. I risultati possono essere salvati come dati (shapefiles puntuali e lineari, csv) o immagini. L'idea iniziale nasce da Marco Zanieri e l'implementazione è di Mauro Alberti. Il plugin è installabile direttamente come nuovo plugin in Quantum GIS.

Esempio di profili prodotti da qProf da DEM SRTM2 ed Aster di una zona della Valnerina (Umbria): in alto i profili topografici,  in basso il corrispondente profilo delle pendenze lungo la traccia: valori positivi indicano salite, negativi discese.

Tool analoghi in Quantum GIS

In Quantum GIS sono disponibili almeno altri due moduli che consentono l'elaborazione di profili topografici: Profile Tool di Boris Jurgiel, Etienne Tourigny e Patrice Verchere, e Profile from Line di Ricardo Garcia Silva. Il primo utilizza come traccia di input un percorso "dinamico", definito dall'utente al momento dell'analisi, su uno o più DEM. L'output grafico e tabellare riguarda le elevazioni ma non le pendenze lungo il profilo. Il secondo determina le elevazioni, a spaziatura costante, a partire da un livello lineare e produce come output un livello puntuale con le elevazioni determinate a partire dallo specifico DEM scelto dall'utente. Non viene creato un output grafico. Entrambi i plugin estraggono le elevazioni dai DEM per nearest neighbor.


Descrizione del plugin

Schermata del plugin.

1 - Path shapefile
La traccia del profilo è costituita da una linea singola e continua  (anche se con andamento variabile), contenuta in uno shapefile.

2 -Use DEMs
Come fonte delle elevazioni possono essere scelti i raster DEM che sono già caricati nel progetto QGIS corrente. Ovviamente in ogni raster ci deve essere una sola banda di dati; in caso contrario il modulo fornirà un messaggio di errore. Raster di notevoli dimensioni possono generare errori di memoria del modulo (MemoryError).

3 - Max sampling distance
L'utente definisce la distanza massima di campionamento lungo la traccia del profilo. Questo plugin mantiene tutti i punti di campionamento definiti nel profilo, sia l'inizio, sia quelli interni, sia quello terminale. La spaziatura risultante può essere quindi più fitta, oltre che irregolare,  della max sampling distance.
Il passo di campionamento determinerà la risoluzione spaziale del profilo, incidendo anche, ovviamente sulla distanza cumulativa 3D. Quest'ultima sarà tanto maggiore quanto minore la distanza (fino al limite di risoluzione del DEM), dato che la topografia ha carattere frattale (e.g. Turcotte, 1997).
Si suggerisce di impostare un valore circa analogo a quello della risoluzione dei DEM utilizzati: per esempio per un DEM ASTER si può inserire un valore di 30 m, per un DEM SRTM2 a 90 m un valore di 90 o superiore.

4 - Calculate profile
I dati DEM vengono letti da Python tramite il binding a GDAL. I valori di elevazione per i punti di campionamento sono interpolati per bilinear interpolation a partire dai centri delle quattro celle che racchiudono il punto da interpolare.
Quando viene modificato un parametro di analisi, come p.e. i DEM usati per il profilo, occorre ricalcolare il profilo, altrimenti verrà visualizzato il profilo coi parametri precedentemente calcolati.

5 - Plot

La visualizzazione dei profili si basa sul modulo Matplotlib, che fornisce gli usuali strumenti di zoom, pan, e salvataggio in formato grafico della visualizzazione (SVG, PDF, jpg, etc.).
Lo zoom non funziona in maniera sincronizzata fra i due grafici che rappresentano le elevazioni e le pendenze. Questo aspetto potrebbe essere risolto nelle prossime versioni.

6 - Save full results as
Tutti i dati vengono salvati come informazioni puntuali in formato csv oppure come uno shapefile 2D. Il valore di inclinazione attributo ad ogni punto di campionamento si riferisce al segmento che lo precede immediatamente. Da notare che se il profilo scelta ha il vertice iniziale al di fuori di un DEM utilizzato, la distanza cumulativa 3D non sarà calcolata per nessun punto.

7 - Export 3D line from DEM
E' possibile salvare le informazioni come dati lineari 3D per uno specifico raster. Il formato di output è sempre shapefile.


Possibili modifiche future
Alcuni spunti per possibili futuri miglioramenti del plugin come tipologie di analisi sono elencati nei link sottostanti.
Esempi di integrazione di Grass con R, per il calcolo della dimensione frattale lungo profili topografici, sono riportati in Raster profile along arbitrary line segments del California Soil Resource Lab. Esempi molto interessanti (ed avanzati) di integrazione tra profili topografici ed informazioni geologiche sono presentati in Grass Gis et R: superposition des couleurs d'un raster quelconque sur un profil topographique élaboré à partir d'un MNT di di Martin Laloux, come pure Python: utilisation des couches vectorielles et matricielles dans une perspective géologique, sans logiciel SIG, sempre di Martin Laloux.

E' prevista anche la possibilità di personalizzare la vestizione del profilo, con scelta del rapporto di scala e delle caratteristiche grafiche.

Bibliografia

Turcotte, D. L., 1997. Fractals and Chaos in Geology and Geophysics. Cambridge University Press, 398 pagine.

Pierre Raybaut rilascia WinPython versione 2.7.3.1

WinPython è una distribuzione scientifica Python per Windows, portabile e sia a 32 sia a 64bit, a differenza di Python(x,y), sempre dello stesso autore.
Si tratta di un progetto open source, interessante per chi vuole disporre di una installazione che può essere spostata di cartella, oppure sul network o su un drive esterno (p.e. una chiavetta).
Oltre a moduli classici come Numpy, Scipy, IPython, guiqwt e Matplotlib, contiene anche Spyder, un ambiente interattivo per sperimentare codice Python, e una suite di Qt: QtDesigner, QtAssistant, QtLinguist e QtDemo.
Il sito è: http://code.google.com/p/winpython/

Andamenti a scala metrica vs. scala kilometrica di faglie. Un esempio di analisi geostrutturale e GIS sulla Linea della Valnerina (Umbria)

Una delle principali fonti di informazioni per stimare l'orientazione di una faglia è la giacitura di mesofaglie e specchi di faglia localizzate in corrispondenza di affioramenti della faglia stessa. La scala su cui queste misure possono essere considerate indicative è dell'ordine di decine-centinaia di m. Su dimensioni maggiori, per esempio km o decine di km, l'andamento di una faglia è suggerito dalla sua traccia topografica della faglia. Conciliare queste due tipologie di informazioni può non essere banale.

In questo post  analizzo i segmenti  compresi fra Schioppo e Meggiano (Umbria) della Linea della Valnerina, un lineamento tettonico situato ad Est della linea Olevano-Antrodoco. Contrasto i dati meso-strutturali acquisiti per la tesi di dottorato con le informazioni derivate dalle tracce topografiche. Per questo secondo aspetto utilizzo dati topografici Aster, sui quali applico la metodologia di intersezione tra DEM e piano geologico (vedi gSurf) che ho presentato in vari post precedenti ed implementato sia come plugin per Quantum GIS, sia come tool Python a se stante.

Nello schema geologico sottostante si vede come la zona studiata sia all'interno di un triangolo definito dalle città di Foligno, Terni e Visso. Subito a SE della zona studiata affiorano sovrascorrimenti col tipico andamento appenninico, circa N-S o NNW-SSE.

Qual'è un problema geologico-strutturale di questo lineamento? La causa dell'andamento obliquo (NNE-SSW) della Linea della Valnerina rispetto ai sovrascorrimenti che sono presenti in altre zone, come per esempio a SSE della linea stessa (prima di incontrare ad oriente la Olevano-Antrodoco).

Oltre all'obliquità, anche la sua inclinazione che localmente può essere molto elevata: per esempio nel classico affioramento di Schioppo si misurano valori di inclinazione di circa 70°, atipici per le faglie inverse. Da queste caratteristiche alcuni Autori, principalmente Decandia negli anni '70 del secolo scorso, hanno proposto che la Linea della Valnerina rappresenti un lineamento normale pre-miocenico coinvolto e riattivato durante l'orogenesi appenninica.

Quello che vado a comparare in questo post è la corrispondenza tra giaciture delle meso-faglie ed andamento topografico della Linea della Valnerina. Sono circa analoghi? Entrambi supportano l'osservazione di andamenti obliqui rispetto al trend usuale dei sovrascorrimenti miocenici, con inclinazioni sensibilmente maggiori di questi ultimi?

I dati meso-strutturali raccolti in quattro stazioni lungo due segmenti della Linea della Valnerina - il segmento di Schioppo-Tassinare (stazioni 106 e 105) a Sud e il segmento di Grotti a Nord (stazioni 104 e 102) - sono presentati nella figura sottostante. Tutti i movimenti rappresentati sono diretti verso NE, di tipo destro-inverso per le stazioni 102, 105 e parte della 106, destri-normali per la 104 e la maggioranza della 106. Nel caso della stazione 106  (affioramento di Schioppo), la maggior parte delle mesofaglie sono ad alto angolo. Nel complesso, i dati mesostrutturali suggeriscono che questi due segmenti della Linea della Valnerina costituiscono dei segmenti obliqui a media-alta inclinazione. Rampe oblique, secondo un'altra terminologia.

Anallizzando più in dettaglio, la zona compresa tra le stazioni 105 e 104, corrisponde ad una biforcazione della Linea della Valnerina, con un segmento circa N-S che si sviluppa a Nord della stazione 105, ed un segmento ad andamento più articolato che comprende le stazioni 104 e 102.

Proviamo ora a determinare le intersezioni topografiche attese utilizzando i valori medi delle mesofaglie osservate nelle quattro stazioni considerate, calcolati e plottati con FaultKin e Stereonet di Rick Allmendinger. Quanto si accordano con gli andamenti topografici mappati?

Nelle mappe sottostanti le tracce teoriche delle faglie sono rappresentate dalle linee bianche, mentre le tracce mappate da quelle rosse. Negli stereonet, i cerchi grandi rappresentano i piani di meso-faglia con strie, i punti neri i poli ai piani ed il triangolo al centro del cerchio blu ll'orientazione del polo medio.

Stazione 106

Per la stazione 106, abbiamo un piano medio con immersione 295.5° ed inclinazione di 70.8°.

Il fit tra traccia mappata e quella teorica è buono sino a circa 1 km a sud della stazione 105, dove l'orientazione della traccia mappata devia nettamente e l'andamento è più articolato, il che suggerisce anche minori inclinazioni dei segmenti di faglia. 

Stazione 105

La stazione 105 ha un piano medio con immersione 271.5° ed inclinazione di 42.5°.

Il fit tra traccia mappata e teorica è limitato e solo locale. Probabilmente questa stazione si trova in corrispondenza di una zona di variazione della orientazione della Linea della Valnerina. 

  

Se si cerca di fittare la traccia topografica della faglia a Nord della stazione 105, si ottiene un valore di circa 259° di immersione e 30° di inclinazione. Questo concorda con l'ipotesi predetta che si stia osservando una variazione sia di direzione sia di inclinazione della faglia, che assume un andamento più 'appenninico' ed una inclinazione più tipica di sovrascorrimenti.  

Stazione 104

La stazione 104 ha piano medio con immersione 39° ed inclinazione 34.4°.

Come si nota nell'immagine sottostante, la traccia teorica derivata dalle mesofaglie non si accorda affatto con quella mappata. Si può quindi ritenere che le giaciture osservate nella stazione strutturale 104 abbiano una validità solo locale, inferiore al km. La faglia potrebbe avere subito un piegamento post-cinematico, forse causato dalla faglia normale pliocenica, ad andamento NNW-SSE, mappata subito ad W della stazione 104.

Stazione 102

Infine la stazione 102 ha valori simili a quelli della 105: immersione di 282° ed inclinazione di 48.3°.

Come si vede nella mappa sottostante, il fit locale è relativamente buono su distanze di alcuni km, ma poi devia nettamente dalla traccia mappata circa 1 km a NE della stazione 104. 

Un fit più esteso, dalla stazione 102 sino alla 105, si ottiene considerando valori di immersione di 284° ed inclinazione di 8°, quindi mantenendo la stessa immersione ma riducendo nettamente l'inclinazione rispetto a quella osservata per le mesofaglie.

A Sud della stazione 105 questo fit cessa, in quanto la traccia teorica ha un andamento che differisce nettamente da quella reale (cf. zona di Schioppo), inoltre sarebbero attesi dei klippen del tetto della faglia nella zona a N-E delle stazioni 105-106. E' quindi da ritenere che nella zona a Sud della stazione 105 la faglia aumenti di inclinazione sino a valori di 70°, che sono compatibili sia con la traccia topografica locale sia con le mesostrutture osservate.

Modello interpretativo dei segmenti

Una possibile struttura dei segmenti analizzati vedrebbe una faglia ad alto angolo nella zona di Schioppo (stazione 106), con inclinazioni di circa 70°, che tende ad assumere a Nord, circa in corrispondenza della stazione 105, valori nettamente più bassi, attorno a 40°. Il segmento a orientazione N-S che si sviluppa a Nord della stazione 105 potrebbe rappresentare uno splay a più alto angolo che si origina nel tetto  della Linea della Valnerina, e che avrebbe inclinazioni sull'ordine di 30°, oltre ad una direzione NNE-SSW.

La prosecuzione verso Est della Linea della Valnerina assumerebbe inclinazioni più ridotte, in media attorno a 10° o poco più, con locali tratti a differente orientazione od inclinazione.

New version of QGIS plugin for DEM/plane intersection analysis

To determine the intersections between a DEM and  a geological plane, one can use a Quantum GIS plugin, qgSurf, that is now available also via the Python plugin manager of QGIS, or it can be downloaded from: http://plugins.qgis.org/plugins/qgSurf/version/0.1.2/

The calculated intersections can now be saved as linear or point shapefiles.
More details at:

• A Quantum GIS plug-in for the determination of plane topographic traces and gSurf: a Python program for the interactive determination of intersections between DEM and planes
  

R da Python: disponibile nuovo installer Rpy2 per Windows

Usare R da Python permette di usufruire di molte funzionalità di analisi e visualizzazione di notevole qualità. Brad Breisfeld ha rilasciato un nuovo installer del modulo Rpy2 per Python in Windows (32 bit), che permette di evitare possibili problemi di configurazione.
A questo scopo è opportuno seguire la sequenza preliminare di configurazione, qui descritta:
https://bitbucket.org/breisfeld/rpy2_w32_fix/src/137d6103a0c2/INSTALL_WINDOWS

L'installer è scaricabile dal link:
https://bitbucket.org/breisfeld/rpy2_w32_fix/issue/1/binary-installer-for-win32

Il progetto di Breisfeld è un fork del lavoro di Laurent Gautier, il creatore di Rpy2.
Documentazione su questo progetto, che consente di sfruttare le funzionalità di R direttamente in Python è presente all'indirizzo:
http://rpy.sourceforge.net/rpy2/doc-2.2/html/introduction.html#getting-started

Intersezioni DEM-piani con gSurf: nuova opzione di salvataggio come linee

Una nuova versione di gSurf stand-alone (richiede Python) è disponibile: consente di salvare le intersezioni tra un piano geologico ed un DEM anche come linee, oltre che come punti. Anche se sperimentale, i risultati del salvataggio come linee sono in generale completi ed equivalenti a quelli come punti.
La versione, testata sia su Linux (Ubuntu) sia su Windows Vista, è disponibile su GitHub oltre che come files compressi alla pagina http://www.malg.eu/gsurf_it.php.

A Quantum GIS plug-in for the determination of plane topographic traces

     Topographic traces of planes can be derived by means of gSurf [1], a Python application that, given a DEM and a plane, determines the intersections (as points) between these two surfaces. This tool is now available also as a plug-in for Quantum GIS. Why Quantum GIS? Because it is a user-friendly, open source and free GIS software, that allows a simple and rapid integration between its core functionalities and Python- and Qt-based modules. 

     This plug-in can be installed by copying the plug-in folder into the standard Python plug-ins folder for Quantum GIS. Since it presents the same interface, tools and theoretical bases as the stand-alone Python application, the reader is referred to previous posts for additional information [1, 2, 3]. It was tested in Quantum GIS 1.7.3 in both Windows Vista and Ubuntu Lucid Lynx (10.4 LTS) and can be freely downloaded here (license: GPL v. 3). 

Fig. 1. Screenshot of the plugin window.

    Next I describe an example of application and result validation using data from the Valnerina zone between S. Anatolia di Narco and Cerreto di Spoleto (Umbria, Central Italy), where Plio-Pleistocene normal faults dissect the previous Miocene thrust-and-fold structure. Geological data derive from my PhD thesis. The example regards a normal fault, for which a structural measure and the topographic trace are available (Fig. 2). The structural measure value is (dip direction, dip angle) = (227°, 58°). However, it cannot be considered representative of the general orientation of this fault, since the computed trace (white line in Fig. 2) does not fit well the mapped fault trace. 

Fig. 2.  Fault traces superposed onto DEM. White dots (mimicking a line) represent the computed intersections for a plane with measured values of dip direction and angle (227°, 58°). The red circle represents the measure location.

     A value of (219°, 61°) seems more appropriate to fit the mapped trace, even with some deviation due to local erroneous fault trace mapping or fault surface deviations from the general orientation (Fig. 3). 

Fig. 3. A plane-solution with a better fitting to the mapped fault trace.
Plane orientation is (dip direction, dip angle) = (210°, 61°).

     In order to demonstrate the correctness of the computed theoretical result, a method routinely used in geology has been applied to the intersections. The highest and lowest intersections between the results and the contour lines derived from the DEM has been extracted, and their vertical and horizontal separations computed. The vertical separation is equal to (1100 – 650) m = 450 m, while the horizontal separation is equal to about 245 m (Fig. 4). This values indicate a dip angle of 61.4°, to be confronted with the theoretical value of 61°. The orientations of the two intersection lines, that are parallel as expected, suggest a dip direction of about 220°, to be compared with the theoretical value of 219°.

We can be therefore confident of the correctness of the application results. 

Fig. 4. Screenshot showing part of the two lines representing the intersection directions with height of 1100 m asl (blue line, top) and 650 m asl (blue line, bottom-left). Their separation is approximately of 245.5 m.

Related posts

[1] gSurf: una applicazione Python per calcolare interattivamente l'intersezione fra piani e DEM.
http://gisoftw.blogspot.com/2012/01/gsurf-una-applicazione-python-per.html 

[2] Calculating the intersections between planes and DEM: a Python implementation
http://gisoftw.blogspot.com/2012/01/calculating-intersections-between.html

[3] Intersezioni tra DEM e superfici planari, un tema di interesse in geologia
http://gisoftw.blogspot.com/2011/12/intersezioni-tra-dem-e-superfici.html

Calculating the intersections between planes and DEM: a Python implementation

I present a new release (vers. 0.1.1) of a Python application that calculates the theoretical intersections between a plane and a DEM (in a raster format). This application could be useful in geological investigations, i.e., for comparing structural measures of geological surfaces with their mapped (or yet to map) topographic traces.

When we define a plane attitude by its dip direction and angle, as well as by a single point point lying in the plane, we can analytically determine the intersection points between these two surfaces.

The theoretical basis of the algorithm are described in the post (in Italian): Intersezioni tra DEM e superfici planari, un tema di interesse in geologia. A case study is presented in:
gSurf: una applicazione Python per calcolare interattivamente l'intersezione fra piani e DEM.

This application allows to save the computed intersections as a point shapefile, for further processing in GIS software, and also images of the map.

A screen shot of the application is below.

To run this program,  the following  Python modules are required:

• numpy
• gdal
• PyQt4
• matplotlib

It has been tested on Windows Vista and Ubuntu Lucid Lynx.

It can be downloaded from the link gSurf: a Python program for the interactive determination of intersections between DEM and planes, where additional information is also available.

gSurf: una applicazione Python per calcolare interattivamente l'intersezione fra piani e DEM

Nota del 29/01/2012: una nuova versione di questo tool è presentata nel post:
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

L'applicazione è stata implementata in Python 2.7.2 con Eclipse + Pydev e viene rilasciata con licenza GPL v. 3.

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 è:

1. 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 );
2. 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);
3. 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;
4. 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.

Intersezioni tra DEM e superfici planari, un tema di interesse in geologia

Nota del 29/01/2012: la più recente versione di questo tool è presentata nel post:
Calculating the intersections between planes and DEM: a Python implementation  


Dato un punto nello spazio 3D, con una orientazione associata (espressa p.e. come inclinazione ed immersione, secondo la convenzione geologica), qual'è l'intersezione della superficie col DEM, ovvero l'andamento della traccia della superficie sul DEM?

Questo problema è di interesse in geologia strutturale, in quanto spesso disponiamo di misure di superfici in stazioni strutturali, misurate come orientazione di un piano.
Possiamo essere interessati a conoscere quale sia la “proiezione” del piano sul DEM, nell'assunzione che la sua orientazione rimanga immutata in un certo intorno della misura.
Se noi non conosciamo l'andamento effettivo della superficie, questa “proiezione” può essere utile per ipotizzare la posizione della traccia nelle vicinanze della misura. Se invece, grazie ad osservazioni di campagna, ne conosciamo l'andamento, possiamo verificare quanto la nostra misura planare sia rappresentativa dell'andamento complessivo della superficie e se nel variarla riusciamo a ottenere un maggiore accordo tra misura e dato di campagna.

Oltre a queste finalità immediate, ce ne possono essere altre: questa determinazione può essere il primo passo per segmentare un volume roccioso, la cui superficie superiore è espressa dal DEM, in più sottoelementi, con superfici di delimitazione rappresentate da più piani o superfici più complesse (p.e. funzioni trigonometriche).
Inoltre, la determinazione delle intersezioni per superfici teoriche permette di disporre di dataset teorici di test per algoritmi che effettuano il procedimento inverso, cioé stimano in maniera (semi)automatica l'andamento di superfici naturali a partire da osservazioni puntuali su una superficie topografica.


Modello analitico

La misura planare rappresenta il caso geometrico più semplice, quindi più facilmente implementabile in un algoritmo e concettualmente più semplice. Questo aiuta la comprensione dell'andamento reale della superficie investigata.
Determinare l'intersezione tra DEM e superficie non è immediato, in quanto di tratta di un problema con dimensionalità superiore a 2, soglia oltre la quale i software GIS attuali sono poco attrezzati.

Il modello concettuale qui usato per calcolare le intersezioni tra DEM e piano nell'algoritmo implementato in Python si basa su semplici considerazioni di geometria analitica.
L'equazione di un piano passante per un punto (x₀, y₀, z₀) e perpendicolare ad una retta con coseni direttori (l, m, n) è (Zwirner, 1983, p. 198):
l (x – x₀) + m ( y – y₀ ) + n ( z – z₀ ) = 0
sviluppabile come:
z = - ( l / n ) (x – x₀ ) – ( m / n ) ( y – y₀ ) + z₀

I coseni direttori di una retta con inclinazione δ e orientazione θ sono (Groshong, 1999, pp. 67-8; vedi anche Cox e Hart, p. 167, eq. 5.1-3):
l = sin θ cos δ
m = cos θ cos δ
n = - sin δ

Considerando le relazioni tra un piano con inclinazione δ_P e immersione θ_P, e la retta ad esso perpendicolare e che punta verso l'alto, si ha:
θ = θ_P
δ = δ_P - 90°
per cui i coseni direttori predetti, espressi rispetto alla giacitura del piano, hanno valori:
l = sin θ_P sin δ_P
m = cos θ_P sin δ_P
n = cos δ_P

Ritornando all'equazione del piano, abbiamo quindi:
z = a ( x – x₀ ) + b ( y – y₀ ) + z₀
con:
a = - l / n = - sin θ_P tan δ_P
b = - m / n = - cos θ_P tan δ_P

I due coefficienti a, b rappresentano le derivate parziali di z rispetto a x e a y.
Essendo θ_P e δ_P costanti nel caso del piano, anche a e b lo sono.

Queste sono le basi analitiche per esprimere un piano in base alla sua giacitura geologica. Ma qual'è la metodologia per il calcolo dell'intersezione tra DEM e superficie planare?
Una procedura semplice può essere quella di calcolare, per ogni cella (e precisamente il suo punto centrale), l'intersezione tra due linee che rappresentano l'andamento locale del DEM e della superficie, sia in direzione parallela all'asse delle x, sia in quella dell'asse delle y.

Consideriamo il caso delle intersezioni lungo sezioni parallele all'asse delle x. Per quelle parallele all'asse delle y le considerazioni sono analoghe.
L'algoritmo determina le equazioni dei segmenti che esprimono l'andamento locale del DEM in base alle coordinate ( x_n, z_n ) e ( x_n+1, z_n+1 ).

Se i segmenti relativi al DEM e al piano lungo l'asse delle x hanno coefficienti angolari ed intercette ( m₁, q₁ ) e ( m₂, q₂ ), l'intersezione ha valore:
x_inters = ( q₂ – q₁ ) / ( m₁ – m₂ )

mentre il valore y è quello della particolare riga considerata. Da ( x, y ) si determina poi z.
Un valore calcolato di intersezione rappresenta una intersezione reale se il suo valore delle x ricade nell'intervallo tra il centro della cella considerata e quello della cella successiva lungo la direzione dell'asse delle x.

Implementazione

Questa prima versione è stata implementata con Python e richiede i moduli numpy e gdal. I DEM corrispondono a matrici di elevazione, il che permette di elaborarli agevolmente in Python con numpy.
L'algoritmo creato non è ottimizzato ed effettua i calcoli necessari per la determinazione dei valori delle intersezioni descritti di seguito per ogni cella del DEM, anche quando questa non presenta nessuna intersezione con la superficie planare. Ma anche in questa versione preliminare, i tempi di esecuzione sono comunque ridotti  grazie alla efficiente gestione delle operazioni su array da parte di numpy.
La versione attuale dell'algoritmo è a linee di codice, senza interfaccia grafica, ed ancora in fase di testing e miglioramento. Verrà rilasciata come software open source ad uno stadio successivo.


Esempio di applicazione su una faglia normale pliocenica della Valnerina (Umbria)

Presento un esempio di applicazione di questo modulo su dati geologici della zona della Valnerina, Umbria (mappati da me nella tesi di dottorato).

La sequenza stratigrafica affiorante è quella umbro-marchigiana, con sedimenti marini che vanno dal Dogger-Malm (Calcari Diasprini) al Miocene inferiore (Bisciaro). Nel Miocene-Pliocene questi sedimenti vengono piegati e fagliati, in corrispondenza della fase compressiva che ha originato la catena appenninica. Nel tardo Pliocene-Pleistocene inizia una fase distensiva, con formazione di faglie normali ad orientazione NW-SE. E' proprio su una di queste faglie che valutiamo la corrispondenza tra misurazioni di terreno e andamento cartografico, cioé traccia della superficie con la topografia.

Consideriamo una faglia normale per la quale si dispone di una misura strutturale, con immersione 227° ed inclinazione 58°. L'andamento mappato della faglia è evidenziato in figura.

Le strie circa E-W nello shaded relief sono dovute ad artefatti tipico del DEM Aster.

Il risultato con i valori della misura non concorda con la traccia mappata della faglia. L'immersione complessiva sembra essere di circa 10° inferiore a quella nel sito.

Testando una nuova intersezione, usando una immersione di 217° e mantenendo l'inclinazione a 58°, si ottiene una traccia che ricalca maggiormente l'andamento complessivo, anche se sembra avere una variabilità locale superiore a quella mappata. Come è noto, tanto più un piano è verticale, tanto più la sua traccia sulla superficie topografica sarà rettilinea. E' possibile quindi testare un piano con la stessa immersione ma inclinazione maggiore.

Con valori di immersione di 217° ed inclinazione 70° si ottiene una traccia della superficie che rispetta in maniera soddisfacente l'andamento mappato della faglia. Discordanze locali possono essere dovute a variazioni della superficie, o anche a interpolazioni erronee della traccia stessa in zone con assenza di dati. Nella porzione sud-orientale della traccia è comunque evidente che si ha una discordanza di circa 20° nell'immersione. Si può quindi ipotizzare che in questo settore la faglia varii la sua orientazione complessiva, oppure che si tratti di un segmento indipendente dalla faglia principale. Questo suggerisce le zone in cui possono essere utili ulteriori controlli di terreno.

Bibliografia

Cox, A., Hart, R.B., 1990. La tettocnica delle palcche. Meccanismi e modalità. Zanichelli, Bologna. 383 pp.
Groshong, R. H. jr., 1999. 3-D Structural Geology. Springer Verlag, Berlin. 324 pp.
Zwirner, G. 1983. Istituzioni di matematiche. Parte seconda. Ed. CEDAM, Padova. 486 pp.