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). 

Thursday 24 November 2011

Streaming delle sessioni plenarie GFOSS del 24 e 25 novembre

Sulla lista GFOSS viene segnalato che è disponibile la visione in streaming delle sessioni plenarie di questi due giorni del convegno di Foggia.
Link:


Interessante post sulle distribuzioni Linux per i GIS

Su CadGis.it è stato pubblicato un elenco delle distribuzioni Linux che comprendono software GIS pre-inseriti, come Grass, QuantumGIS e simili.
Il link è: http://www.cadgis.it/distribuzioni-linux-per-il-gis

Sunday 13 November 2011

Un plugin Quantum GIS per il calcolo della pendenza direzionale

La pendenza direzionale è un parametro topografico che può essere utile per evidenziare l'andamento della topografia lungo direzioni specifiche. In alcuni software GIS, come p.e. Mirone, è implementato il calcolo della pendenza direzionale lungo direzioni uniformi nello spazio. Un caso più generale si ha quando le orientazioni lungo le quali calcolare la pendenza variano invece nello spazio.

In post precedenti avevo presentato uno script Python (La pendenza topografica lungo direzioni variabili) ed una prima versione di plugin per Quantum GIS (Directional slope along changing directions: a Python plug-in for Quantum GIS) che consentivano di calcolare la pendenza lungo orientazioni variabili.

In questo post presento una nuova versione con maggiori funzionalità, sempre come plugin per Quantum GIS. Le maggiori funzionalità consistono nel permettere anche il calcolo della pendenza su intervalli di orientazioni uniformi (p.e., da 0° a 360° con step di 10°), oltre alla più usuale pendenza massima.
I metodi per il calcolo della pendenza sono quelli di Zevenbergen & Thorne (1987) e di Horn (1981), i due metodi che producono le migliori stime di pendenza (Burrough & McDonnell, 1998).

Rispetto alle versioni precedenti la possibilità di usare un grid delle orientazioni variabili nello spazio è stata resa più semplice. Non è più richiesto infatti che esso abbia le stesse caratteristiche geografiche del DEM, ma può essere più esteso o anche meno esteso, avere una dimensione della cella inferiore o superiore a quella del DEM. L'algoritmo automaticamente ricampiona, per nearest neighbour, i valori di orientazione al DEM sul quale poi calcolare lo slope.

Una schermata della finestra del modulo è nella figura sottostante.



Il modulo è ulteriormente descritto e scaricabile alla pagina http://www.malg.eu/pendenzadirezionale.php



Tuesday 25 October 2011

Monday 17 October 2011

Calcolare le linee di flusso in campi vettoriali: due implementazioni open-source in Python e C++

Il tema della determinazione delle linee di flusso di dati GIS era già stato considerato nei post Metodi di Runge-Kutta per la determinazione di linee di flusso nei GIS e Generare linee di flusso nei GIS.

In questo post presento due differenti implementazioni in Python ed in C++, entrambe open-source e free, del classico algoritmo di Runge-Kutta-Fehlberg (RKF 45) per il calcolo della traiettoria di flusso in un campo di velocità variabile nello spazio ma “costante” nel tempo.  I due algoritmi calcolano il percorso di ipotetiche particelle puntuali all'interno del campo di velocità durante un intervallo di tempo prefissato. Queste implementazioni funzionano direttamente con dati nei formati classici GIS (shapefiles e raster letti da GDAL). Sono state realizzate e testate su Windows Vista e XP. Fatto salve eventuali modifiche minori, dovrebbero funzionare anche su Linux, ricompilando, nel caso dell'implementazione in C++, il codice sorgente.

Le due implementazioni, con codice sorgente, eseguibile, un dataset di esempio ed una breve descrizione teorica, sono disponibili alla pagina  www.malg.eu/streamlines.php (in inglese). 

Entrambe queste implementazioni si basano su OGR per leggere le posizioni iniziali dei punti  da uno shapefile.
Il campo di velocità è definito dalle sue componenti cartesiane vx e vy, conservate in due distinti livelli raster. Mentre lo script Python legge solo raster in formato Arc/Info ascii grid, il programma C++ legge i formati raster resi disponibili da GDAL. E' cruciale che i livelli raster delle componenti delle velocità abbiano le stesse caratteristiche geometriche e geografiche (uguale numero di righe e colonne, stessa dimensione geografica delle celle, identico posizionamento nello spazio).

Il risultato consiste in uno shapefile, sempre puntuale, i cui punti rappresentano le traiettorie delle particelle per intervalli di tempo successivi. Gli intervalli di tempo possono avere durata variabile, in quanto il metodo RKF 45 consente di scegliere dei passi adattivi per ridurre gli errori stimati al di sotto di una soglia massima.


Thursday 29 September 2011

Creare una semplice mappa web geologica con Leaflet


Leaflet è una libreria Javascript che permette di creare e visualizzare mappe interattive, sia su desktop sia in dispositivi mobili. Dal punto di vista dell'utente, un suo aspetto apprezzabile è una buona velocità di caricamento e navigazione. Per gli sviluppatori è importante la sua facilità di sviluppo di mappe.

Leaflet permette di visualizzare dati sulla base cartografica OpenStreetMap. Possiamo utilizzare dati WMS e geoJASON. GeoJSON è un formato vettoriale, testuale, che può essere letto direttamente da Javascript. File GeoJSON possono essere creati a partire per esempio da shapefile con strumenti come OGR, che è disponibile in FWTools.

In questo post consideriamo come esempio la creazione di una mappa web della geologia di una zona in Umbria (Valnerina). In questa mappa è possibile avere informazioni sulle formazioni geologiche e visualizzare due sezioni geologiche associate alle loro tracce in mappa.

Quello che descriverò sono solo alcuni aspetti della creazione della mappa. La documentazione completa e gli esempi di Leaflet sono disponibili nel sito di Leaflet: http://leaflet.cloudmade.com.
Il codice sorgente della mappa può essere visualizzato dalla URL http://www.malg.eu/webmaps/leaflet_test/ll_01.html. La mappa è  stata testata su versioni Win recenti di Chrome, Firefox, Opera, Safari e IE (9). L'immagine sottostante rappresenta un esempio di schermata .




Passi preliminari

Carichiamo la libreria Leaflet sul server nel quale vogliamo creare la nostra mappa. Registriamoci come utente Leaflet presso CloudMade e richiediamo una API key per la nostra mappa web.


Configurazioni della pagina html

La API key viene incorporata nella linea di codice html che definisce la URL di collegamento a CloudMade:

var cloudmadeUrl = 'http://{s}.tile.cloudmade.com/YOUR-API-KEY/997/256/{z}/{x}/{y}.png'

Colleghiamo la libreria Leaflet all'interno della pagina web che creiamo, nel suo head:

<link href="../../Leaflet/leaflet.css" rel="stylesheet"></link>
<!--[if lte IE 9]><link rel="stylesheet" href="../../Leaflet/leaflet.ie.css" /><![endif]-->

La seconda linea di codice è necessaria per far sì che la mappa web sia visualizzata anche in IE (io ho testato solo con la versione 9).

La struttura della pagina è ispirata a quella di http://kothic.org/js/ nella parte che riguarda la legenda.

Creazione del livello geologico

In questo esempio considereremo dati della geologia di una zona della Valnerina (Umbria). E' una mappa geologica che deriva dai dati di terreno che ho raccolto nel corso della mia tesi di dottorato in geologia strutturale. Questi dati sono in formato shapefile poligonale, in coordinate geografiche (cioè non proiettati in UTM e simili). Vengono convertiti in formato geoJSON con OGR, utilizzando FWTools con il comando ogr2ogr, per esempio:

ogr2ogr -f "GeoJSON" geologia.js geologia_geogr.shp geologia_geogr

Apriamo il file in formato JSON appena creato e aggiungiamo al suo inizio la definizione della variabile (evidenziata in giallo) che conserva i dati geologici:

var geologia = {
"type": "FeatureCollection",
"features": [
{ "type": "Feature", "properties": { "AREA": 20906.211000, "PERIMETER": 791.990000, "GEOL_CL_": 4.000000, "GEOL_CL_ID": 0.000000, "FORMATION": "Scaglia Rossa" }, "geometry": { "type": " ....

Questo file geoJSON deve essere collegato alla pagina web che crea la mappa, incorporandone il riferimento nell'head del file html:

<script src="geologia.js" language="javascript"></script>

Dopo questo passo, si può editare la pagina html, inserendo tutte le funzioni in uno script Javascript inserito o collegato nell'head della pagina.

La variabile Javascript che rappresenta i dati geologici viene richiamata nel codice Javascript che definisce la mappa, facendole creare l'oggetto geoJSON (evidenziata in giallo). Lo stile del livello viene definito attraverso la variabile formation_colors (in verde).

var geojson = new L.GeoJSON();
geojson.on('featureparse', function(e) {
// you can style features depending on their properties, etc.
var popupText = 'Formazione: ' + e.properties.FORMATION + <br \>;
if (e.layer instanceof L.Polygon) {
e.layer.setStyle({color: 'gray', weight: 1, fillColor: formation_colors[e.properties.FORMATION]});
}
e.layer.bindPopup(popupText);
});

geojson.addGeoJSON(geologia);
map.addLayer(geojson);

Nel codice sovrastante viene anche definito il contenuto di una finestra di poup che si attiva al click su ogni poligono del livello, oltre allo stile di rappresentazione  delle formazioni geologiche. Quest'ultimo si basa sulla proprietà FORMATION (e.properties.FORMATION) definita nel file geologia.js per ogni poligono, e su una 'colormap' che assegna un colore ad ogni formazione geologica:

var formation_colors = {
'Depositi fluviali':'#cFF',
'Depositi quaternari':'#cCF',
'Bisciaro':'#Fc6',
'Scaglia Cinerea':'#Fcc',
'Scaglia Variegata':'#0cc',
'Scaglia Rossa':'#3c6',
'Scaglia Bianca':'#cc0',
'Marne a Fucoidi':'#c06',
'Maiolica':'#063',
'Calcari Diasprini': '#309'
}


Definizione delle tracce delle sezioni

L'ultimo aspetto che descrivo è la creazione delle due sezioni geologiche. Il codice relativo alla prima delle due sezioni è:

var p1 = new L.LatLng(42.76554, 12.79581);
var p2 = new L.LatLng(42.83028, 12.91872);
var section_a = [p1, p2];
var sezione_a = new L.Polyline(section_a, {weight: 1, color: 'blue', opacity: 0.5});
map.addLayer(sezione_a);
sezione_a.bindPopup("<img style='width: 600px;' src='ims/sezione_a.jpg' alt='sezione A' />", {maxWidth: 650});

Vengono creati due punti con le coordinate geografiche che corrispondono agli estremi della sezione. Questi due punti vanno a creare una lista, la quale serve da base per la creazione di una Polyline, che viene aggiunta alla mappa e alla quale viene collegato un popup, con la funzione bindPopup, che apre una immagine che rappresenta la sezione specifica. E' opportuno definire la maxWidth della finestra di popup in maniera tale che contenga completamente l'immagine collegata.



Monday 29 August 2011

Esempi di visualizzazioni statistiche su dati GIS con R


R è un software di statistica open-source, con un numero di metodologie di analisi e visualizzazione ben superiore a quello tipico dei software GIS.

Un apparente svantaggio di R è l'uso di righe di comando piuttosto che di una interfaccia grafica. Ma questo consente di copiare e riutilizzare più volte in R una successione di comandi anche complessa.

Anche se è possibile usufruire di R all'interno di Python, Grass e altri, la maniera più semplice, soprattutto per chi si avvicina la prima volta ad R, è quella di utilizzare di R in maniera a se stante.

Un documento di buon livello e free su R è An Introduction to R: Software for Statistical Modelling & Computing di Kuhnert e Venables. Su temi di analisi spaziale, un testo avanzato è Applied Spatial Analysis with R di Bivand, Pebesma e Gómez-Rubio. Per un quadro generale  dei packages esistenti per l'analisi spaziale in R è utile  il sunto che si trova in CRAN Task View: Analysis of Spatial Data.

Per la lettura di dati vettoriali possono essere usati i packages maptools, shapefiles, foreign, RArcInfo, mentre per i raster i pacchetti Rgdalraster. Esistono anche dei packages specifici per Saga e Grass. Un pacchetto per la geostatistica è gstat.

Lo scopo di questo post è di presentare come sia possibile, in maniera semplice, analizzare dati GIS. Usiamo quindi i packages che risultano più immediati.


Dati vettoriali
In generale i dati vettoriali sono conservati come shapefile o possono essere facilmente convertiti in tale formato con i GIS. A quel punto è facile importarli in R, o leggendo l'intero shapefile, oppure il solo file dbf dello shapefile, nel quale sono conservati gli attributi descrittivi.
Nel primo caso è possibile anche plottare mappe ma questi utilizzi sono quasi sempre eseguibili nei GIS, che invece mancano di funzionalità statistiche avanzate.
Consideriamo quindi il più semplice package foreign, che consente di leggere un qualsiasi file dbf come quelli degli shapefile.

Useremo lo shapefile dei comuni italiani del 2001 (com2001_s.shp), rilasciato dall'ISTAT, come dataset di esempio. A questo shapefile, contenente l'informazione della popolazione, sono stati aggiunti un campo col nome della regione di appartenenza, oltre a due campi che contengono area e perimetro dei comuni.

Dopo avere installato in R il package foreign, lo si carica con il comando:

library(foreign)

Possiamo quindi leggere il nostro file dbf:

comuni <- read.dbf("com2001.dbf")

avendo se necessario l'accortezza di spostarsi nella cartella che contiene lo shapefile (File → Cambia directory).

A questo punto è utile chiamare una funzione di resoconto sui dati letti, per verificare nomi, tipo e statistiche dei campi importati:

summary(comuni)

Verrà visualizzato un elenco di elenco dei campi, con le loro statistiche.



Quelli che utilizzeremo sono i campi POP2001, AREA e COD_REG.

I codici delle regioni sono nella tabella successiva.



1 PIEMONTE
2 VALLE D'AOSTA
3 LOMBARDIA
4 TRENTINO-ALTO ADIGE
5 VENETO
6 FRIULI-VENEZIA GIULIA
7 LIGURIA
8 EMILIA-ROMAGNA
9 TOSCANA
10 UMBRIA
11 MARCHE
12 LAZIO
13 ABRUZZO
14 MOLISE
15 CAMPANIA
16 PUGLIA
17 BASILICATA
18 CALABRIA
19 SICILIA
20 SARDEGNA

I singoli campi possono essere richiamati inserendo, fra il nome del dataset e quello dello specifico campo (come compare nella lista prodotta dal comando summary), il simbolo '$'.
Per esempio:

summary(comuni$AREA)

A questo punto possiamo creare grafici di vario tipo, oltre che effettuare analisi statistiche. Seguono alcuni esempi di grafici.

Proviamo a sintetizzare con boxplot le distribuzioni dei residenti per comune, raggruppati in base alla regione (espressa dal codice Istat).


boxplot(comuni$POP2001 ~ comuni$COD_REG, main='Comuni italiani - popolazione 2001', xlab='Regioni (codice Istat)', ylab='Residenti')

Il codice (semplificato) POP2001 ~ COD_REG esprime la suddivisione dei valori di popolazione POP2001 in base a valori con significato di categoria conservati in COD_REG.



La statistica è dominata dalle grandi città (Roma, Milano, Napoli, Torino etc.), con le loro popolazioni sull'ordine del milione di abitanti. La grande maggioranza dei comuni finisce schiacciata verso il basso, per il loro numero molto ridotto di residenti.

Per rendersi conto della situazione dei comuni medio-piccoli, è quindi necessario restringere il dataset, con una semplice operazione di filtro sui comuni con meno di 100.000 abitanti ([comuni$POP2001 < 100000]) che va ripetuta per ogni campo coinvolto (due in questo caso). 


boxplot(comuni$POP2001[comuni$POP2001 < 100000] ~ comuni$COD_REG[comuni$POP2001 < 100000], col='pink', main='Comuni italiani con popolazione 2001 < 100.000 residenti', xlab='Regioni (codice Istat)', ylab='Residenti')




Dal boxplot si nota come Puglia (16) e Toscana (9) abbiano distribuzioni meno frammentate rispetto alle altre regioni.
Proviamo a comparare tramite istogrammi la densità di due differenti regioni che presentano una differente distribuzione delle popolazioni dei comuni, il Piemonte (1) e la Puglia (16), considerando l'intervallo inferiore a 200,000 residenti.
I due grafici verrano accolti in una singola una finestra ad una riga e 2 colonne usando il comando par(mfrow=c(1,2)).

pop_piemonte <- comuni$POP2001[comuni$COD_REG == 1 & comuni$POP2001 < 200000]
pop_puglia <- comuni$POP2001[comuni$COD_REG == 16 & comuni$POP2001 < 200000]
par(mfrow=c(1,2))


hist(pop_piemonte, col = 'pink', xlim=c(0,2e5), breaks=10, main='Comuni Piemonte < 200,000 res.', xlab='Residenti', ylab='Numero')

hist(pop_puglia, col = 'pink', xlim=c(0,2e5), breaks=20, main='Comuni Puglia < 200,000 res.', xlab='Residenti', ylab='Numero')



Si nota l'abbondante numero di comuni con residenti inferiori a 10,000 abitanti nel caso del Piemonte oltre alla bassa frequenza di comuni con popolazione maggiore di 20,000 abitanti. La situazione è nettamente più equilibrata in Puglia, oltre a corrispondere ad un numero ben inferiore di comuni, di circa un ordine di grandezza.

Per confronto delle aree degli stessi comuni, possiamo utilizzare i seguenti comandi:
aree_piemonte <- comuni$AREA[comuni$COD_REG == 1 & comuni$POP2001 < 200000]
aree_puglia <- comuni$AREA[comuni$COD_REG == 16 & comuni$POP2001 < 200000]
par(mfrow=c(1,2))

hist(aree_piemonte/1e6, breaks=10, xlim=c(0,6e2), col='orange', main='Comuni Piemonte < 200,000 res.', xlab='Area (kmq)', ylab='Numero')

hist(aree_puglia/1e6, breaks=30, xlim=c(0,6e2), col='orange', main='Comuni Puglia < 200,000 res.', xlab='Area (kmq)', ylab='Numero')


Si conferma anche per le aree comunali (con residenti < 200k) lo squilibrio a favore di una regione come il Piemonte rispetto all'esempio della Puglia: il Piemonte presenta una frammentazione territoriale di un ordine di grandezza superiore rispetto alla situazione della Puglia. Occorre tenere presente che le superfici di Piemonte e Puglia sono circa comparabili (19,600 kmq la Puglia, 25,400 kmq il Piemonte).


Dati raster

La struttura dei dati raster è semplice e possono essere esportati come dati tabulari, importabili in R con un comando come read.dbf o read.csv.
Se invece si preferisce importarli come livelli raster, sono utilizzabili packages come rgdal e raster.

Rgdal potrebbe richiedere l'installazione aggiuntiva di GDAL, se non è già presente perché già installata da altri software GIS.

Negli esempi che seguono useremo i dati SRTM2 a 3 secondi-arco, collegati fra loro e poi ritagliati all'estensione del Piemonte. Questi dati presentano alcune zone con dati mancanti.

library(rgdal)

e quindi si possono leggere i formati raster riconosciuti da gdal, come l' ESRI ASCII grid:
piem_elev <- readGDAL("srtm2_orig_geog.asc")

ed ottenere informazioni di sintesi con:
summary(piem_elev)

o visualizzarlo come mappa:
image(piem_elev, col = terrain.colors(12))


Per leggere ESRI ASCII grid esiste anche la funzione readAsciiGrid del package maptools.



Nel formato ESRI grid è presente una sola banda. Possiamo visualizzare un istogramma della sua distribuzione con:
hist(piem_elev$band1, col='orange')




oppure la sua empirical cumulative distribution function (ECDF) con:

plot(ecdf(piem_elev$band1), main='ECDF elevazioni Piemonte', xlab='elevazione')



Questo è solo un breve esempio delle numerose, possibili visualizzazioni con R.