Friday 30 December 2011

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 (x0, y0, z0) e perpendicolare ad una retta con coseni direttori (l, m, n) è (Zwirner, 1983, p. 198):
l (x – x0) + m ( y – y0 ) + n ( z – z0 ) = 0
sviluppabile come:
z = - ( l / n ) (x – x0 ) – ( m / n ) ( y – y0 ) + z0

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 – x0 ) + b ( y – y0 ) + z0
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 ( xn, zn ) e ( xn+1, zn+1 ).


Se i segmenti relativi al DEM e al piano lungo l'asse delle x hanno coefficienti angolari ed intercette ( m1, q) e ( m2, q2 ), l'intersezione ha valore:
xinters = ( q2 – q1 ) / ( m1 – m2 )

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.

Thursday 15 December 2011

'Vector field processings': un plugin Quantum GIS per l'analisi di campi vettoriali 2D

Questo plugin Python per Quantum GIS fornisce una serie di funzionalità per derivare parametri di campi vettoriali 2D, come divergenza, modulo del rotore e linee di flusso. Due esempi di dati su cui sono applicabili queste tecniche sono i flussi glaciali e le colate di frana.

Determinazione di linee di flusso nel ghiacciaio Reeves (Antartide).
Dati periodo 2001-2003, da Biscaro, (2010).

I moduli descritti sono stati realizzati inizialmente come moduli Python o C++ a se stanti, e presentati nei post Un esempio di analisi di campo vettoriale con software open source: il ghiacciaio David (Antartide orientale) (per divergenza e rotore) e Calcolare le linee di flusso in campi vettoriali: due implementazioni open-source in Python e C++ (per le linee di flusso).

In questo post presento alcuni dettagli tecnici del plugin, oltre ai risultati di un test sulla crescita degli errori cumulativi nell'interpolazione di linee di flusso (o più precisamente pathlines), mentre in quello precedente descrivo un esempio di applicazione a dati di flussi glaciali antartici (ghiacciaio Reeves).



Schermata del plugin


I parametri calcolati dal plugin sono elencati nella tabella seguente.

Vector parametersMagnitude gradientsPathlines
magnitudealong x axispathlines from start points
orientationsalong y axis
divergencealong flowlines
curl module

Alcuni operatori su campi vettoriali sono intuitivi, come la magnitudine o l'orientazione del flusso. La loro derivazione è possibile anche con una applicazione delle funzioni standard di map algebra presenti nei software GIS. Più complesse sono le derivazioni del gradiente e del rotore. Si tratta di due operatori fisici che rappresentano rispettivamente la tendenza a espandersi/contrarsi oppure a ruotare di un flusso. Un flusso con divergenza positiva tende ad espandersi, a contrarsi nel caso opposto. Moduli positivi del rotore indicano rotazioni anti-orarie, se negative rotazioni orarie.

Il gradiente spaziale della magnitudine del campo lungo le orientazioni di flusso, come pure lungo le orientazioni degli assi cartesiani, ci informa sulle variazioni di magnitudine sperimentate dai flussi lungo quelle specifiche orientazioni.

Per la derivazione delle linee di flusso sono richieste tecniche più avanzate: si utilizzano metodi numerici per la risoluzione di Ordinary Differential Equations (ODE). Uno fra i più diffusi, grazie all'accuratezza dei risultati prodotti, è il metodo di Runge-Kutta-Fehleberg. Non ripeto qui le formule, presenti in testi di matematica per l'ingegneria (come  Advanced Engineering MathematicsKreyszig, 2006), oltre che in rete (fra cui nel post Calcolare le linee di flusso in campi vettoriali: due implementazioni open-source in Python e C++). Con questi metodi è possibile calcolare le linee di flusso sia per flussi stazionari sia per flussi che variano nel tempo.
L'implementazione nel plugin si riferisce al caso più semplice, i flussi stazionari. Una maggiore disponibilità di dataset GIS vettoriali non stazionari renderebbe utile l'estensione anche a campi variabili nel tempo.


Gli errori cumulativi nella interpolazione di pathlines

Il metodo di Runge-Kutta-Fehleberg (RKF 45) è utilizzato, oltre che per la sua accuratezza, anche perchè consente di stimare l'errore connesso al passo temporale utilizzato. 
Considerando che riducendo o aumentando il passo temporale ridurremmo o aumenteremmo l'errore connesso alla stima dello spostamento, è possibile variare il passo temporale in maniera da ottenere un errore stimato prossimo ma non superiore ad una soglia da noi stabilita. Questo errore stimato si riferisce ad ogni singolo passo di interpolazione.

Ma cosa avviene quando noi interpoliamo per 1000, 2000 o più passi successivi? Come aumentano gli errori totali nella stima, e quindi quanto sono affidabili i nostri risultati quanto interpoliamo su lunghi periodi temporali?
Questo è un aspetto a cui accenno in questa seconda parte, utilizzando un esempio analitico di movimento circolare uniforme per il quale posizioni al variare del tempo sono perfettamente predicibili e quindi le posizioni stimate con l'algoritmo sono confrontabili con quelle esatte.
Quanto segue è un esempio che fornisce delle indicazioni su quale possa essere la crescita nel tempo degli errori cumulativi di interpolazione delle posizioni.
Consideriamo un moto circolare uniforme che avviene su un raggio di 100 m, con una velocità tangenziale costante di 200 m/s. Il periodo è quindi di circa 3.14 s.

Campo vettoriale 'circolare' (visualizzato in base alla magnitudine) e pathlines determinate  dal plugin.
Il punto in alto al centro rappresenta la posizione iniziale.

I parametri utilizzati per il test di interpolazione consistono in step di 0.1 s, un tempo totale di 5000 s, ed un errore massimo ammesso per step di 1e-6 m. L'algoritmo modifica la durata degli steps in 0.025 s, per avere errori costanti di 4.0e-8 m. Per raggiungere il limite temporale richiesto di 5000 s il numero totale di step è uguale a 200 000. 


Confrontando i valori stimati di posizione con quelli teorici, vediamo che in prima approssimazione l'errore cumulativo reale cresce in maniera lineare, proporzionale al numero di step
Errore cumulativo reale nella posizione al variare del numero di passi di interpolazione (linea arancione).
Il valore stimato di errore per singolo step è rappresentato dalla linea orizzontale (rossa). 


I valori del tasso di crescita dell'errore complessivo (1.4e-9 m/step, vedi figura sottostante) sono di circa un ordine di grandezza più piccoli dell'errore stimato per gli step individuali (4.0e-8 m/step). Il tasso di crescita non è costante, ma presenta oscillazioni, relativamente limitate e di origine non chiara. 
Si può supporre che in prima approssimazione l'errore cumulativo sia inferiore al prodotto tra il numero di step effettuati ed il valore medio stimato dell'errore per singolo step.
Ovviamente questa relazione risulta valida per questo specifico caso di analisi e dovrebbe essere confermata o approfondita per altri casi teorici.

Variazioni dell'errore cumulativo / numero di step di interpolazione.



Pagina di download del plugin

Si trova alla url: www.malg.eu/vectorparameters.php
L'installazione è quella usuale per i plugin per Quantum GIS.





Analisi dei flussi glaciali nella zona di disancoraggio del Reeves Glacier (Antartide)

I flussi glaciali sono un esempio di dati vettoriali che possono essere analizzati con il plugin 'Vector field processings' per Quantum GIS, descritto nel post successivo, che permette di determinare divergenza, rotore, gradienti della magnitudine e linee di flusso.

L'esempio qui presentato si riferisce al ghiacciaio antartico Reeves, che sbocca sulla costa del mare di Ross. I dati degli spostamenti sono stati determinati  nella tesi di dottorato in Scienze Polari di Debbie Biscaro a titolo Applicazioni di metodi di telerilevamento per lo studio dei ghiacciai di Terra Nova Bay e della Cook Ice Shelf, Antartide orientale.

I dati sono stati verificati e interpolati per creare un campo di velocità definito dalle sue componenti vx e vy. L'interpolatore utilizzato è la tecnica spline (Regularised Spline with Tension and Smoothing, RSTS, Mitasova & Mitas, 1993, Mitasova et al., 1995), implementata in Grass. I parametri di interpolazione della RSTS sono stati scelti in maniera da ridurre al minimo le statistiche di cross-validation.
Una descrizione più dettagliata della modalità di creazione è presentata in Glacial flow analysis with open source tools: the case of the Reeves Glacier grounding zone, East Antarctica.

I dati interpolati di velocità, come pure il plugin, sono scaricabili dalla pagina: www.malg.eu/vectorparameters.php


Il Reeves Glacier nella zona di confluenza nella Nansen Ice Sheet.
I flussi sono orientati da NW verso SE. I punti rappresentano le misure di spostamento disponibili su un intervallo di 400 giorni nel periodo 2001-2003. La linea di disancoraggio del ghiaccio continentale dal substrato roccioso è rappresentata dalla linea bianca. L'immagine satellitare di base è Landsat (banda pancromatica).

Nelle due immagini sottostanti sono raffigurate i valori della velocità nelle loro componenti x (a sinistra) ed y (destra). Le componenti tendono ad evidenziare la continuità nei flussi lungo la direzione rappresentata dall'asse considerato.

Componenti delle velocità di spostamento lungo l'asse delle x (E-W). L'unità di misura è in m ed il periodo temporale considerato è di 400 giorni (2001-2003).
Componenti delle velocità di spostamento lungo l'asse delle y (N-S). L'unità di misura è in m ed il periodo temporale considerato è di 400 giorni (2001-2003). 
Il primo parametro di interesse è la magnitudine del campo di velocità, che evidenzia i settori con anomalie positive o negative di velocità. Queste sono localizzate soprattutto nella parte nord-occidentale, dove il ghiaccio di origine continentale è ancorato al substrato. Nella porzione orientale il ghiaccio galleggia sulle acque marine e le velocità hanno una maggiore uniformità. Si nota la forte anomalia positiva delle velocità immediatamente a Sud di Andersson Ridge, appaiata ad una anomalia negativa al suo bordo meridionale. Una possibile causa consiste nella presenza di un ostacolo subglaciale che rallenta i flussi in corrispondenza dell'anomalia negativa, costringendo i flussi a nord a scorrere velocemente in una sezione ristretta.

Magnitudine della velocità nella zona di disancoraggio del Reeves Glacier (Antartide).
Gli spostamenti, in m, si riferiscono ad un periodo di 400 giorni (2011-2003).


Il modulo del rotore ed il gradiente sono altri due parametri vettoriali stimabili. Il modulo del rotore evidenzia anomalie longitudinali in corrispondenza delle zone di transizione fra settori con differenti velocità di flusso. L'ostacolo subglaciale precedentemente ipotizzato presenterebbe al suo bordo settentrionale una fascia di anomalie negative (rotazioni orarie) e a quello meridionale una corrispondente fascia di anomalie positive (rotazioni anti-orarie). La lingua del Reeves sarebbe bordata a Nord-Est da una forte anomalia positiva che segnala il limite con il ghiaccio quasi fermo a Nord. A Sud-Ovest di Teall Nunantak  si sviluppa una anomalia negativa.
Anomalie longitudinali ad alta frequenza spaziale sono presenti anche nella zona orientale, dove il ghiaccio è fluttuante. L'esatta natura di queste anomalie (se esse corrispondano a segnali naturali oppure siano artefatti legati alla derivazione ed interpolazione dei dati) è meno semplice da interpretare. Se fossero segnali naturali potrebbero indicare una tendenza alla suddivisione in più sottoelementi longitudinali della lingua, nel momento in cui essa inizia a galleggiare.

Le anomalie nella divergenza ugualmente tendono a concentrarsi in corrispondenza degli ostacoli glaciali emersi e sommersi, ma mentre nel caso del rotore lo sviluppo è laterale a questi e con sviluppo longitudinale, le anomalie di divergenza hanno sviluppo trasversale agli ostacoli e si posizionano a monte o a valle di questi. Un aspetto da notare nel caso della divergenza è che nella zona orientale sono evidenti alternanze di anomalie positive e negative come fasce trasversali ad alta frequenza.

Moduli del  rotore.
Valori positivi indicano rotazioni anti-orarie, positivi rotazioni orarie.   
Divergenza dei flussi glaciali.
Valori positivi indicano espansione, negativi contrazione.  

Il gradiente della velocità lungo le orientazioni di flusso evidenzia le stesse strutture precedentemente descritte, sia quelle a minore frequenza tipiche della zona occidentale, ancorata al substrato roccioso, sia quelle a maggiore frequenza, visibili soprattutto nella porzione orientale, fluttuante nel Mare di Ross.
Anomalie positive marcano la parte frontale della zona con forte anomalia della magnitudine, a Sud di Andersson Ridge, e la parte a valle del probabile ostacolo subglaciale. Quest'ultima corrisponde anche all'inizio di una zona dove il ghiaccio continentale, iniziando a galleggiare, si frattura in blocchi fra loro saldati da ghiaccio marino.


Gradienti spaziali della velocità (magnitudine) lungo le orientazioni di flusso.  Le anomalie positive più sviluppate sono  associate alle morfologie del substrato roccioso a Sud di Andersson Ridge. Nella porzione orientale sono ben evidenti anche fasce di anomalie positive e negative ad elevata frequenza e di natura incerta (se artefatti o segnali naturali). 

Un altro aspetto di interesse è il tragitto percorso da particelle all'interno del campo: il loro nome è pathlines e nel caso di flussi costanti nel tempo esse sono parallele alle streamlines (vedi discussione in Wikipedia).

Il plugin permette di determinare le pathlines per flussi costanti nel tempo, a partire da punti iniziali definiti in un livello puntuale. Il loro tragitto permette di evidenziare le zone in cui i flussi convergono o divergono, per esempio a causa di ostacoli come evidente nell'immagine sottostante. Utilizzando le pathlines è anche possibile stimare lo spostamento nel corso del tempo.
Nell'immagine sottostante si nota che i tracciati delineati sono percorsi in intervalli temporali compresi fra circa 200 e 300 anni, a seconda dei differenti  settori.


Pathlines a partire dai punti specificati in alto a sinistra, assumendo flussi costanti durante l'intervallo considerato (compreso tra 200 e 300 anni).