Metodi di Runge-Kutta per la determinazione di linee di flusso nei GIS

Questo post prosegue sul tema del precedente: la determinazione delle linee di flusso da campi vettoriali in una modalità GIS-friendly. Software come Paraview consentono di effettuare questa determinazione, ma sia l'input sia l'output sono in formati non usuali per i GIS.

Queste tecniche numeriche rientrano nel tema dell'integrazione di ODE (Ordinary Differential Equations). 

Le fonti bibliografiche per la creazione degli algoritmi sono Computational Physics di P.O.J. Scherer, Advanced engineering Mathematics di E. Kreyszig e A Primer on Scientific Programming with Python di H.P. Langtangen. Esistono inoltre molte fonti di informazione accessibili online (vedi p.e. Wikipedia).

Nel precedente post ho presentato uno script Python che calcola le pathlines per campi vettoriali costanti nel tempo. Quello script preliminare si basava sul metodo più semplice fra quelli proposti in letteratura, cioè il metodo di Eulero esplicito (forward Euler), la cui formula è:

u _k+1 = u _k +  Δt f( u _k , t _k)                                                            (eq. 9.23 in Langtangen, 2009)

dove:
 u _k :  valore della variabile allo step k,
u _k+1 : il suo valore allo step successivo
Δt : incremento temporale considerato
f( u _k , t _k) : la funzione che ci fornisce la derivata di u _k al tempo t _k, nel caso specifico la velocità.

Essendo un campo vettoriale con componenti x e y, la formula soprastante verrà applicata sia all'una sia all'altra componente. La derivata viene stimata calcolando l'incremento della variabile durante l'intervallo di tempo considerato. Nell'algoritmo implementato, i valori delle componenti in un particolare punto vengono stimati per interpolazione bilineare.

I risultati ottenibili con questo metodo non sono stabili ed è preferibile utilizzare altri metodi, il più noto dei quali è il metodo di Runge-Kutta, con le sue varianti.

Prima di considerare le altre metodologie, è necessario stabilire come si possa confrontare la qualità dei risultati fra i diversi metodi. In Computational Physics è utilizzato come esempio quello di un campo vettoriale centrale (esempi sono il campo gravitazionale e quello elettrostatico) con orbite circolari (Fig. 1a). Un metodo di interpolazione di linee di flusso in un tale campo dovrebbe generare orbite strettamente circolari. Si può dimostrare sia teoricamente sia sperimentalmente che il metodo di Eulero esplicito (forward) non è stabile e determina delle orbite che man mano hanno un raggio sempre maggiore (Fig. 1b). Metodi con una maggiore stabilità sono l'improved Euler method (midpoint) e soprattutto i metodi di Runge-Kutta del terzo e quarto ordine ed il metodo di Runge-Kutta-Fehlberg (Fig. 1c). Questi producono risultati molto più stabili, in cui le orbite si degradano molto più lentamente con i vari cicli di iterazione.

Fig. 1. A) campo centrale con orbite circolari. B) L'orbita interpolata con il metodo di Eulero esplicito non mantiene un raggio costante ma cresce rapidamente con il tempo. C) Metodi come Runge-Kutta 4 e Runge-Kutta-Fehlberg sono stabili e le orbite determinate rimangono circolari al passare del tempo.
Ridisegnato da Schroeder & Martin,  fig. 1.17, in 'The visualization handbook', Academic Press.  

Dei metodi inseriti nella nuova versione dello script Python (mid-point, Runge-Kutta del terzo e quarto ordine, Runge-Kutta-Fehlberg) inseriamo le formule del più complesso e preciso fra questi, cioé il metodo di Runge-Kutta-Fehlberg (RKF) (riprese da Kreyszig, 2006):
y_n+1 = y_n + h ( γ₁k₁ +  ... +  γ₆k₆)
con vettore di coefficienti:
γ  = [ 16/135    0    6656/12825    28561/56430    -9/50    2/55]
e stimatori:
k₁ = f(x_n, y_n)
k₂ = f(x_n + (1/4) h, y_n + h (1/4) k₁)
k₃ = f(x_n + (3/8) h, y_n + h ( (3/32) k₁ + (9/32) k₂) )
k₄ = f(x_n + (12/13) h, y_n + h ( (1932/2197) k₁ - (7200/2197) k₂ + (7296/2197) k₃) )
k₅ = f(x_n + h, y_n + h ( (439/216) k₁ - (8) k₂ + (3680/513) k₃ - (845/4104) k₄) )
k₆ = f(x_n + (1/2) h, y_n + h ( - (8/27) k₁ + (2) k₂ - (3544/2565) k₃ + (1859/4104) k₄ – (11/40) k₅) )

Confrontando i risultati del metodo di Eulero esplicito con il Runge-Kutta-Fehlberg ottenuti per un campo con orbite circolari con la nuova versione dello script Python, si nota appunto l'instabilità del primo e la stabilità del secondo.

Fig. 2. Confronto tra risultati del metodo di Runge-Kutta-Fehlberg e di Eulero esplicito su un punto iniziale (quadrato)  in un campo centrale con orbite circolari: il circolo interno rappresenta vari cicli calcolati dal primo metodo, la spirale crescente i corrispondenti cicli prodotti dal secondo.

Consideriamo ora l'applicazione di questi due metodi, agli estremi come precisione, su un dataset naturale, già analizzato in post precedenti: i flussi del ghiacciaio David in corrispondenza della sua zona di disancoraggio dal substrato per formare la lingua Drygalski (Antartide). Questo dataset deriva dalla tesi di dottorato in Scienze Polari di Debbie Biscaro: [http://www.mna.it/italiano/Didattica/dott_polare/Tesi_abstract/Biscaro_XXI_2010.pdf]   

Il campo di velocità è stato dedotto confrontando due immagini Landsat fra loro coregistrate del 2001 e del 2003, con una separazione temporale di 400 giorni. Le velocità sono state determinate tramite il software open-source Imcorr: http://nsidc.org/data/velmap/imcorr.html .

Applichiamo lo script usando il metodo di Eulero esplicito e di Runge-Kutta-Fehlberg, confrontando i risultati fra di loro e rispetto alle lineazioni glaciali visibili dalle immagini satellitari Landsat.
Il prodotto tra i due metodi non differiscono in maniera sensibile. Alla scala della mappa sottostante sono indistinguibili (Fig. 3). Fra i punti terminali della traccia più settentrionale in Fig. 3, la differenza è di soli 9.3 m su un percorso totale di circa 15.3 km.
Queste basse differenze sono presumibilmente dovute alle traiettorie fondamentalmente rettilinee seguite dai flussi.

Fig. 3. Mappa delle linee di flusso determinate con il metodo di Runge-Kutta-Fehlberg. A questa scala il risultato è indistinguibile da quello del metodo di Eulero esplicito.  Sfondi: mosaico satellitare LANDSAT e grid della magnitudine delle velocità. Mappa creata con QGIS Copiaco. 

Quali sono le relazioni con le lineazioni di flusso visibili nel mosaico satellitare LANDSAT?
Confrontando le traiettorie da RKF con le lineazioni, si vedono chiaramente delle discordanze dove le traiettorie sono particolarmente variabili (Fig. 4). Questo potrebbe essere dovuto a varie cause: determinazione non ottimale del campo di velocità, inadeguatezza del metodo RKF, variazioni nel tempo del campo di velocità glaciale, tali per cui le streamlines determinate da RKF differiscono dalle pathlines che sono evidenti nelle immagini glaciali.
Non abbiamo indizi per discernere fra queste e altre ipotesi possibili. L'ultima ipotesi però ci sembra verosimile, in quanto è noto che i flussi glaciali possono variare su tempi anche abbastanza brevi ed inevitabilmente le lineazioni glaciali subiscono processi di deformazione che possono farle divergere dalle streamlines "istantanee". Metodi di analisi della deformazione glaciale potrebbero aiutare a valutare meglio questa possibilità.

  

Fig. 4. Mappa delle traiettorie determinate con il metodo di Runge-Kutta-Fehlberg confrontabili con le lineazioni di flusso glaciale visibili nel mosaico stellitare LANDSAT. Mappa creata con QGIS Copiaco.

Lo script Python è ancora in fase di test e inserimento di ulteriori funzionalità. Chi fosse interessato, può richiedere la versione attuale a alberti.m65@gmail.com.

Frank Warmerdam, mente di GDAL/OGR, passa a Google

Nel blog di Frank Warmerdam c'è l'annuncio del suo prossimo passaggio a Google.
Dopo essere stato libero consulente GIS soprattutto sulla sua creatura, GDAL/OGR, Warmerdam andrà ad assumere il ruolo di GIS Data Engineer in Google.
Interessante è anche la sua descrizione del processo di interview con i recruiters Google.
GDAL/OGR è una libreria GIS open source, per dati raster e vettoriale, con un importante ruolo in vari programmi GIS e può essere usata usata anche direttamente negli script Python.
Warmerdam pianifica di dedicare parte del suo 20% free time in Google alla sua creatura... 


http://fwarmerdam.blogspot.com/2011/06/joining-google.html

Rilasciato QuantumGIS 1.7.0, "Wroclaw"

Oggi il team di sviluppo, che comprende anche Paolo Cavallini di Faunalia,  ha distribuito la nuova versione, denominata Wroclaw: http://www.qgis.org/component/content/article/127-qgis-1-7-release.html
Ora può quindi iniziare un processo di verifica di possibili bug da parte di noi utilizzatori, oltre ovviamente all'uso produttivo. Per gli utenti Windows, è disponibile un installer stand-alone: http://www.qgis.org/wiki/Download.

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.