Friday 3 June 2011

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.



No comments: