Un esempio di analisi di campo vettoriale con software open source: il ghiacciaio David (Antartide orientale)

Post di Mauro Alberti alberti.m65@gmail.com e Debbie Biscaro debbiemail@libero.it

I campi vettoriali possono essere utili per descrivere flussi glaciali, idrici o spostamenti di volumi rocciosi. Purtroppo questo tipo di dati non sono ancora trattati esplicitamente dai software GIS. Con qualche escamotage possono essere processati parzialmente, a differenza della ricchezza di funzioni di software come Paraview, VisIt o Matlab. Linguaggi come Python, ricchi di funzioni numeriche e con interfacce a librerie GIS standard, consentono di ovviare in parte a queste limitazioni e lavorare direttamente su formati dati tipici dei GIS. I risultati possono essere visualizzati negli stessi software GIS o in Paraview e simili.

In questo post presentiamo uno script Python per il calcolo dei parametri vettoriali (rotore, divergenza,  accelerazione lungo le linee di flusso) di un campo 2D e lo applichiamo ad un caso naturale di flussi glaciali in Antartide. Lo script fa parte di un set di tool in Python per l'analisi di flussi glaciali e non necessita di specifici software GIS.


Parametri vettoriali

Divergenza e rotore sono due operatori vettoriali molto utilizzati in fisica e ingegneria. Questi due operatori corrispondono all'applicazione su di un campo vettoriale dell'operatore differenziale nabla:

La divergenza è un valore scalare derivato dal prodotto scalare tra nabla e il campo vettoriale nell'intorno di un punto (x,y,z) derivabile:

Valori positivi di divergenza indicano un flusso che tende ad espandersi (pensiamo all'aria riscaldata che si dilata), mentre valori negativi un flusso che tende a convergere e a comprimersi.

Il rotore è un vettore derivante dal prodotto vettoriale tra nabla e campo vettoriale:

Esso esprime la componente rotazionale dei flussi.

Per dati bidimensionali, come avviene frequentemente per i GIS, le due formule sovrastanti si semplificano. Per la divergenza si considereranno le sole derivate parziali rispetto a x e y. Per il rotore si avrà un vettore costantemente verticale (parallelo a k), il cui modulo sarà dato dall'ultima componente della formula precedente. L'unico parametro di interesse sarà quindi il modulo del rotore. Il segno associato al modulo sarà positivo per rotazioni in senso anti-orario e negativo per rotazioni in senso orario. Valori nulli del modulo indicheranno flussi che non ruotano.

Un altro parametro di possibile interesse è l'entità delle accelerazioni lungo le linee di flusso. Per calcolare questo parametro, utilizziamo una formula modificata per la determinazione della pendenza direzionale da un DEM (Neteler & Mitasova, 2008, eq. A.27):

dove alpha rappresenta l'orientazione locale della linea di flusso che è facilmente derivabile da vx e vy. Il risultato è un campo scalare.


Implementazione e uso dello script Python

Script: curl_div_acc_03.py

Con questo script è possibile leggere i valori delle componenti cartesiane di un campo vettoriale 2D da due grid e determinare divergenza, rotore ed accelerazioni lungo i flussi sfruttando la funzione gradient presente nel package numpy che calcola le derivate parziali lungo gli assi cartesiani. I dati vengono esportati sia nel formato Arc/Info ASCII grid, in maniera da poterli importare in GIS come Quantum GIS o Saga, sia nel formato VTK per processamenti in Paraview o VisIt.

Lo script Python è stato sviluppato e testato con Python 2.6 all'interno dell'ambiente Python (x,y) per Windows. Richiede numpy come modulo Python (presente di default in Python (x,y)). Viene rilasciato con licenza GPL v. 3.

Un esempio di lancio dello script dal command prompt è:
  python curl_div_acc_03.py par.txt
dove par.txt è un file testuale che contiene i parametri di input ed output dell'analisi in sei differenti righe, esempio:
velxclp.asc // ascii grid with x cartesian components of vector field
velycl.asc // ascii grid with y cartesian components of vector field
div.asc // ascii grid with divergence values (output)
curl.asc // ascii grid with curl magnitudes (output)
acc.asc // ascii grid with acceleration values (output)
result.vtk // vtk file storing all input and output data (output)

Il suo input è costituito da due grid che conservano le componenti x e y del campo vettoriale, in formato ARC/INFO ASCII . Produce in output tre grid in formato ARC/INFO ASCII, che contengono la divergenza, il modulo del rotore e l'accelerazione lungo le linee di flusso, oltre ad un file in formato VTK, in cui oltre al campo vettoriale originario, sono conservate divergenza, modulo del rotore e accelerazione.


Esempio di analisi: la parte terminale del ghiacciaio David (Terra Vittoria settentrionale, Antartide)

Consideriamo ora un esempio di analisi di un flusso vettoriale, in cui consideriamo le velocità di flusso del ghiacciaio David, in Antartide, in corrispondenza della sua zona di disancoraggio dal substrato roccioso.

I dati dei flussi derivano da Biscaro (2010), tesi di dottorato della quale potete leggere l'abstract (in inglese).
I dati di partenza sono costituiti da una coppia di immagini Landsat, rispettivamente del 2001 e 2003 (intervallo temporale di 400 giorni fra le due), fra loro co-registrate.
I flussi sono stati derivati con il software libero IMCORR (Scambos et al., 1992), che si basa su tecniche di cross-correlation tra una coppia di immagini co-registrate per derivare gli spostamenti glaciali relativi.
I risultati prodotti da IMCORR sono stati validati e ne sono stati rimossi, per quanto possibile, gli errori di co-registrazione delle due immagini. Nei risultati sono quindi presenti zone prive di valori.
Per il calcolo dei parametri vettoriali di questi flussi, lo script Python richiede grid continui.  E' stato quindi ricostruito un campo vettoriale continuo interpolando separatamente le due componenti cartesiane del campo. Il metodo di interpolazione scelto è il Regularised Spline with Tension and Smoothing (RSTS) di Mitasova & Mitas (1993) e Mitasova et al. (1995), implementato in GRASS.  Essendo una tecnica spline, appare adatta per superfici relativamente liscie e prive di asperità. Consente di incorporare un fattore di smoothing che regola il grado con cui i valori originali vengono rispettati nell'interpolazione, utile per misure che incorporano errori.
I parametri di interpolazione scelti per questo esempio sono quelli risultati ottimali per la vicina zona del Reeves Glacier (analizzata nell'articolo sottoposto): tensione uguale a 100 e smoothing di 0.1. Come si vedrà dall'analisi dei parametri vettoriali risultanti, la tecnica ed i parametri usati per questa zona non sono completamente ottimali e dovrebbero essere riconsiderati per migliorare la qualità delle interpolazioni risultanti.
I due grid interpolati delle velocità lungo gli assi x e y hanno costituito il dataset di input per lo script Python. I risultati sono stati analizzati e visualizzati con Mirone, Quantum GIS e ParaView.

Presentiamo ora alcuni aspetti dei risultati che riteniamo di interesse.

La semplice visualizzazione delle magnitudini delle velocità evidenzia la presenza  nella porzione orientale di un ostacolo roccioso sub-glaciale che suddivide e devia i flussi del ghiacciaio David in un ramo principale meridionale ed uno secondario settentrionale. I flussi, orientati da oriente verso occidente (verso destra nelle figure sottostanti) proseguono con elevate velocità nel corpo principale del ghiacciaio che va a costituire la lingua glaciale Drygalski.

Mappa delle velocità interpolate, sovrapposte a mosaico satellitare Landsat. Le velocità più elevate delineano l'andamento del ghiacciaio David, che scorre da Ovest verso Est (verso la destra in figura). Le velocità sono in metri per un periodo di 400 giorni (2001-2003). I cerchi rappresentano le misure originarie prodotte da IMCORR  e con valori ritenuti validi. Figura creata con Mirone e Quantum GIS.

La magnitudine del rotore evidenzia zone di forti rotazioni nelle direzioni dei flussi, evidenziate anche dalle lineazioni di flusso glaciale. Si ricorda che valori positivi indicano rotazioni anti-orarie, mentre quelli negativi rotazioni in senso orario. Poichè queste variazioni sono legate all'andamento degli ostacoli nel substrato roccioso, la mappa della magnitudine del rotore esalta le principali variazioni di quest'ultimo. Esse risultano collocate in corrispondenza dell'ostacolo a sinistra, e della fascia sub-orizzontale che limita il bordo sinistro del ghiacciaio David. Nel resto della zona, i flussi tendono ad esere irrotazionali.

Mappa della magnitudine del rotore, sovrapposta ad immagine satellitare Landsat.
Creata con Mirone e Quantum GIS.

Analizzando le mappe delle divergenze e accelerazioni si riconoscono artefatti che suggeriscono una qualità non ottimale delle interpolazioni effettuate.
Osserviamo infatti che alcune variazioni nei parametri ricavati sono localizzate in corrispondenza dei margini delle zone con dati validi, il che fa ritenere che le variazioni naturali vengano “attratte” verso questi limiti artificiali. In altri casi, la tecnica spline potrebbe avere originato zone con valori estremi dei parametri vettoriali, sia massimi sia minimi, che non sembrano essere richiesti dai dati originari.

La mappa delle accelerazioni lungo i flussi è illustrata nella figura seguente. Notiamo nella porzione orientale la stretta corrispondenza tra variazioni dell'accelerazione e e limiti dei dati di velocità validi. Altre variazioni rappresentano probabilmente effetti naturali. Si tratta per esempio della fascia di variazioni lungo il bordo sinistro del ghiacciaio David e di quelle nella porzione nord-occidentale della mappa.

Mappa delle accelerazioni, sovrapposta al mosaico satellitare Landsat.
I cerchi rappresentano le misure originarie. Creata con Mirone e Quantum GIS.

La zona orientale presenta molti artefatti, come è visibile nella mappa sottostante della divergenza. Le brusche variazioni di valori tendano ad essere generate e/o attratte verso i limiti nel campionamento. Dove le osservazioni sono discontinue (porzione superiore della figura sottostante) l'interpolazione sembra avere prodotto delle coppie artificiali massimi-minimi.

Particolare della zona orientale per il parametro divergenza. Si riconosce l'influenza della distribuzione dei punti di misura disponibile sui valori di divergenza risultanti. Creata con Quantum GIS.

Gli operatori che si basano su differenze sono molto sensibili a rumori nel segnale analizzato. Come si vede dagli esempi precedenti, tendono anche ad enfatizzare strutture legate al processamento dei dati, come la distribuzione delle misure note usate nell'interpolazione. Nel caso di vettori, l'interpolazione è stata effettuata due volte, sulle due componenti cartesiane, e quindi è presumibile che l'effetto sia stato maggiormente enfatizzato.
Come regola generale, nelle analisi di dati geografici è importante poter confrontare i risultati di analisi con la distribuzione spaziale delle misure originarie, in maniera tale da discriminare la natura dei risultati ottenuti.


Qual'è il grado di affidabilità delle direzioni di flusso dedotte?

Dal campo vettoriale è possibile ricostruire le linee di flusso tramite il  metodo di Runge-Kutta, implementato in Paraview. La validità nella delineazione nei flussi risultanti dall'interpolazione può essere testata confrontando con le linee di flusso evidenti nelle immagini satellitari. Nel caso specifico le strutture glaciali nella zona circostante l'ostacolo subglaciale permettono di riconoscere un accettabile grado di congruenza a grande scala tra campo vettoriale ricostruito e strutture glaciali. 

Linee di flusso determinate dal campo vettoriale interpolato con il metodo di Runge-Kutta, implementato in Paraview. I flussi glaciali sono diretti da sinistra verso destra.

Analizzando nel dettaglio la zona che presenta le maggiori variazioni di direzioni, quella orientale, notiamo che la corrispondenza non è comunque ideale. Anche in questo caso una ricerca più approfondita di metodi e parametri di interpolazione ottimali potrebbero produrre risultati di migliore qualità.

Particolare della zona investigata in cui si riconoscono differenze nelle linee di flusso ricostruite dal campo interpolato (linee in rosso) e lineazioni glaciali che suggeriscono i flussi glaciali nell'immagine satellitare Landsat. Mappa in Paraview.

Bibliografia

Biscaro, 2010. Applicazioni di metodi di telerilevamento per lo studio dei ghiacciai di Terra Nova Bay e della Cook Ice Shelf, Antartide orientale (Applications of remote sensing methods to the study of Terra Nova Bay and Cook Ice Shelf glaciers, Eastern Antarctica). Ph.D. Dissertation, Siena University, Italy, 173 pp.

Mitasova, H., Mitas, L., 1993. Interpolation by Regularised Spline with Tension: I. Theory and implementation. Mathematical Geology 25, 641-655.

Mitasova, H., Mitas, L., Brown, B.M., Gerdes, D.P., Kosinovsky, I., 1995. Modeling Spatially and Temporally Distributed Phenomena: New Methods and Tools for GRASS GIS. International Journal of GIS 9, Special Issue on Integration of Environmental Modeling and GIS, 443-446.

Neteler, M., Mitasova, H., 2008. Open Source GIS: A GRASS GIS Approach. Springer, Berlin, DE, 406 pp.

Scambos, T. A., Dutkiewicz, M. J., Wilson, J. C., Bindschadler, R. A., 1992. Application of image cross-correlation to the measurement of glacier velocity using satellite image data. Remote Sensing of Environment 42, 177-186.