'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 parameters	Magnitude gradients	Pathlines
magnitude	along x axis	pathlines from start points
orientations	along y axis
divergence	along 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 Mathematics, Kreyszig, 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.

Generare linee di flusso nei GIS

Descrivere e analizzare gli spostamenti di particelle e volumi nel tempo è un compito che possiamo incontrare in ambiti naturalistici, geologici, fisici ed ingegneristici, fra gli altri. Esempi sono i flussi glaciali ed atmosferici, i movimenti tettonici e gli spostamenti di animali, individui e mezzi.

Un caso relativamente semplice, che illustro qui, è quello dei flussi glaciali. Essi possono essere rappresentati in 2 dimensioni ed essere considerati relativamente costanti in tempi brevi, sull'ordine di settimane e mesi (anche se possono talora essere presenti variazioni brusche).

Con i GIS la velocità di spostamento può essere rappresentata da un set di grid, ciascuno dei quali contiene le informazioni di velocità lungo una componente cartesiana, o con formati scientifici più avanzati, anche da un singolo dataset (penso a netCDF e HDF, anche se non li ho testati personalmente).

Il calcolo vero e proprio degli spostamenti in un intervallo di tempo prescelto non è però banale. In base alle mie conoscenze, non esistono funzionalità apposite nei GIS come si possono invece trovare in software di visualizzazione scientifica come ParaView e MayaVi, che permettono di calcolare le streamlines di un campo vettoriale.

Streamlines, pathlines (e anche streaklines) sono rappresentazioni dei flussi che coincidono fra loro nel caso di flussi costanti nel tempo, caso che stiamo considerando (vedi in Wikipedia: http://en.wikipedia.org/wiki/Streamlines,_streaklines,_and_pathlines )

Anche nel caso di ParaView, la funzione di determinazione delle streamlines non consente di precisare un intervallo di tempo da considerare. Inoltre la stessa manipolazione del risultato per importarlo nei GIS non è banale.

Una determinazione di pathlines per flussi costanti con la variabile tempo esplicita può inoltre rappresentare il primo passo per determinare le timelines, la cui visualizzazione può essere molto utile nell'analisi delle variazioni spaziali dei flussi.

Con Python ho creato una versione preliminare (autonoma da specifici software GIS, ma con input ed output direttamente leggibili da QuantumGIS, Saga, etc.), ancora soggetta a verifiche e migliorie, che permette di determinare le pathlines di punti (rappresentati in shapefiles) in base ad un campo vettoriale (conservato in 2 grid delle componenti vx e vy, in formati Arc/Info ASCII).

Chi fosse interessato può richiederla ed eventualmente collaborare nel testing e miglioramento – alberti.m65@gmail.com

Cosa richiede la determinazione delle pathlines di un campo vettoriale costante?

Innanzitutto la stima delle componenti di velocità per punti arbitrari all'interno del dominio del grid.

In questa prima versione questa stima viene effettuata con una interpolazione bilineare applicata sia alle componenti della velocità lungo l'asse delle x, sia delle y.

Secondariamente, una metodologia per interpolare il movimento nel tempo. Si tratta di un problema equivalente alla risoluzione di equazioni differenziali ordinarie (ODE). Uno dei metodi più usati, ben descritto in letteratura e che fornisce ottimi risultati, è il metodo di Runga-Kutte, con varie opzioni (vedi: http://it.wikipedia.org/wiki/Metodi_di_Runge-Kutta ). Questo è il metodo implementato in Paraview.

Esistono metodi più semplici e che forniscono risultati sub-ottimali, come per esempio il Forward Euler method (vedi p.e. Langtangen, 2009, p. 513), che banalmente estrapola il nuovo valore in base alla derivata temporale nel punto di partenza (corrispondente alla velocità di spostamento, nel nostro caso), senza tenere conto dei valori delle derivate in punti intermedi (come invece vale nel caso del metodo di Runga-Kutte). In questa prima versione ho scelto il metodo più semplice, il Forward Euler method.

Dati di test ed esempio utilizzati derivano dai risultati prodotti da D. Biscaro nella sua tesi di dottorato sui flussi glaciali nella Terra Nova Bay (Antartide). La sensatezza dei risultati può essere verificata anche con le lineazioni di flusso evidenti dalle immagini telerilevate (Landsat) della stessa zona, che riguarda la zona dove il David Glacier si disancora dal substrato a formare la lingua Drygalski.

Nell'immagine sottostante una mappa a falsi colori delle velocità di flusso (blu: basse, rosso: elevate) è sovrapposta ad un mosaico satellitare Landsat. I punti per i quali verranno calcolate le pathlines sono rappresentati dai quadrati a sinistra (in giallo).

Nell'immagine sottostante i percorsi derivati dallo script creato evidenziano le variazioni di orientazioni e velocità nei flussi che concordano con quanto evidenziato sia dalle lineazioni glaciali sia dalle magnitudini delle velocità. I percorsi totali si riferiscono ad un periodo totale di circa 82 anni, con spaziatura dei punti corrispondente a 400 giorni. In questa immagini i singoli punti sono distinti solo nel caso dei flussi nella porzione in basso, che presenta velocità maggiori. 

Questa prima implementazione sembra quindi essere efficare nel derivare le pathlines di un flusso costante nel tempo, anche se utilizza un metodo relativamente semplice come il Forward Euler. Possibili successive implementazioni potrebbero basarsi su metodi più precisi, come il Runge-Kutta, oltre a consentire anche la determinazione delle timelines. 

Bibliografia

Langtangen, H.P., 2009. A primer on scientific programming with Python. Springer.