20 Ottobre 2010 - Prima Giornata Mondiale della Statistica - World Statistics Day

La Prima Giornata Mondiale della Statistica. Indetta dalle Nazioni Unite.
Il link ufficiale: http://unstats.un.org/unsd/wsd/
In Italia l'Istat si occupa dell'organizzazione:  http://www.istat.it/istat/eventi/2010/giornata_mondiale_statistica/
organizzando anche il barcamp Sharing data & Statistical knowledge.
Uno dei tavoli del barcamp sarà dedicato ai Sistemi Informativi Geografici, con "facilitatore" Sandro Furieri, presidente Gfoss. 
E' previsto uno streaming: http://video2.caspur.it/barcamp/streaming.html

E' morto Benoit Mandelbrot, il padre dei frattali

Anche se nei GIS il concetto di frattale è ancora poco applicato, è stato un concetto rivoluzionario per l'economia, le scienze fisiche, la geomorfologia (reticolati idrografici e topografie frattali), oltre che anche in campo artistico.
Il padre di questa teoria, Benoit Mandelbrot, polacco emigrato in Francia e poi negli USA, è morto oggi.

http://www.nytimes.com/2010/10/17/us/17mandelbrot.html?_r=1&hp

Spostamenti nei ghiacciai

Una coppia di immagini satellitari o aeree di una zona, riprese in due periodi differenti, può raccontare molto sulle variazioni avvenute in quell'arco di tempo. Mentre per riprese a latitudini basse e intermedie, le variazioni potranno riguardare la vegetazione, i cambiamenti urbani o i movimenti del suolo, per le alte latitudini entra in gioco la dinamica della criosfera. Si potranno così monitorare gli spostamenti lungo i ghiacciai, le variazioni di estensione del ghiaccio marino o delle fronti di ice shelves (ghiaccio continentale galleggiante sul mare). Considerando la dinamica dei flussi glaciali, dalle immagini teleriprese è possibile determinare gli spostamenti lungo i ghiacciai anche nell'arco di tempo di pochi anni. Vari sono i metodi di determinazione. Per zone che presentano solo spostamenti senza significative rotazioni e deformazioni interne, è possibile usare un software open source, creato per la determinazione dei flussi glaciali. Si tratta di Imcorr, del quale trovate i riferimenti in http://nsidc.org/data/velmap/imcorr.html. L'articolo che lo presenta è Scambos et al. (1992). Bibliografia addizionale su Imcorr è presente in http://www.geo.unizh.ch/~kaeaeb/glims/glims.html#Anchor-Glacier-23240

E' un programma a linee di comando, quindi privo di interfaccia grafica, che può essere usato nei sistemi Unix e Linux o Linux-like, come Cygwin inWindows. Imcorr è basato sui linguaggi C e Fortran. C gestisce l'input e l'output dei dati e dei risultati, oltre al passaggio dei dati verso le subroutine Fortran che effettuano i calcoli. La determinazione degli spostamenti presumibili avviene con metodi di cross-correlation a loro volta basate su trasformate di Fourier. Il requisito fondamentale per i dati di input di Imcorr è che le due immagini da comparare siano fra loro già co-registrate ed abbiano le stesse dimensioni (uguale numero di righe e di colonne). La scelta dei parametri ottimali di dimensione delle sotto-immagini entro cui effettuare la ricerca è forse la parte più gravosa dell'intera procedura di elaborazione.

Un esempio di applicazione di Imcorr a dati satellitari Landsat è rappresentato dalla ricerca di I. Pino sui ghiacciai antartici della Terra Vittoria, in Antartide, reperibile all'indirizzo http://amsdottorato.cib.unibo.it/1177/

Un'altra tesi di dottorato basata su derivazione di velocità di flussi glaciali per zone dell'Antartide orientale è quella di D. Biscaro, il cui abstract è scaricabile da http://www.mna.it/italiano/Didattica/dott_polare/Tesi_abstract/Biscaro_XXI_2010.pdf

Too many elevations?

Having to manage elevation data sets with hundreds of thousands or more of records is becoming a standard, but not easy task, even with the most powerful GIS software. This can be due both to the difficulty of interpolating large data sets, and to irregular spatial distributions of data that can lead to produce DEM with artifacts, particularly evident when deriving morphological parameters (slope, aspect and so on).

As a practical example, we consider a data set of elevations for a sector of Antarctica, deriving from ICESat satellite data measured over a period from 2003 to 2008. It consists of nearly 2,800,000 data showed in the map below. Data density along satellite tracks is high (every 172 meters), with several different tracks nearly covering each other, while the spacing across track is much higher, for instance in the top-left of the map this spacing is between five and ten kilometers. In many subareas of this zone, the local topography presents periodic height variations (not visible in this map), due to the presence of glacial megadune fields. Megadunes have wavelengths of several kilometers, amplitudes of several meters and lateral continuity of tens of km. Clearly, DEM deriving from these data are much more affected by local structures in correspondence of tracks.

A possible remedy for this situation is the use of spatial filters that reduce and uniform as much as possible the data density. A straightforward methodology is the use of a grid-based filter to the values of height. For each grid cell, a filter calculates a single height value, based on some functions that process the records falling in the cell. In general algorithms calculate the average of records, whether of height or other analysed properties. Usually the position attributed to the calculated value is at the cell center. Although results calculated in this way can be considered fundamentally correct, we note three negative aspects: the use of the mean introduces a smoothing of the values; erroneous data affects the average value; positioning the resulting value at the cell center can change the variographic properties of the data set.

I present a short program in C + + command line (without GUI), which filters the data in a grid, determines the median values (when there are three or more values) and preserves the original location of the median value. Program (compiled with CodeBlocks in Windows), source code and parameter and elevation example files are available here. Session parameters are in a text file (example: 'param.txt'), and the program loads this file to read the parameters. This program processes ICESat data, but obviously it can be applied to any elevation data, provided that they share the same, basic textual input structure, i.e. with space-separated record id (integer), x, y, z, structured in a column format (see example in file 'elevations.txt'). 

The use of the median avoids the smoothing of elevations, as well as reducing the impact of inaccurate, extreme values. The median record maintains its original position, obviously falling within the borders of the cell in question. 

The algorithm creates two textual outputs. The former is a grid in the ESRI ASCII format that stores the number of records in each grid cell. Open source GIS like Saga or Quantum Gis can read this format. The latter is a file in table format that stores the median elevation values together with the original spatial location of the median value.

Running this program on the previous Antarctic data set with 1 km cell size reduces the amount of data to 70,988 from 2,783,551, i.e. with a 97.4% reduction of the data number, without substantially changing the available data structure, as the map below, representing the filtered data, evidences.

As another example, if we use a 250 meters cell size, the data number reduces to 444,532 (about 16% of the original amount). The filter also eliminates many spikes in the frequency histogram (below, blue: original data, orange: filtered data), probably corresponding to artifacts related the track oversampling of local areas, since no natural reasons for their existence.

Mappa delle offerte di lavoro GIS

The GIS Jobs Clearinghouse Map rappresenta in mappa offerte di lavoro GIS:
http://www.gjc.org/map.html.

E basandosi sulla mappa l'Europa non appare brillare...

GoogleMars

Istruzioni per visitare Marte.

Get Flash to see this player.

Operatori vettoriali nei GIS: il gradiente

Nei GIS i procedimenti e le visualizzazioni di dati vettoriali non sono molto diffuse, fatto salvo il caso del calcolo delle pendenze topografiche. I GIS calcolano implicitamente un campo vettoriale, tramite due distinti campi scalari, l'esposizione (aspect) e la pendenza (slope), che esprimono rispettivamente l'orientazione della massima inclinazione rispetto al nord della mappa ed il valore di inclinazione rispetto al piano orizzontale.

Questo campo vettoriale corrisponde al gradiente di un campo scalare, l'elevazione. Il gradiente viene calcolato applicando ad una superficie scalare ϕ(x,y,z), supposta continua e derivabile in ogni punto, l'operatore differenziale vettoriale nabla:

Nel caso di una superficie topografica la superficie è funzione di due soli variabili spaziali, x e y, per cui z = ϕ(x,y) e il suo gradiente vale:

Il risultato è un campo vettoriale bi-dimensionale su un piano orizzontale. Per ogni punto del dominio spaziale, il vettore orizzontale punta nella direzione di massima pendenza della topografia ed ha un modulo proporzionale alla pendenza di questa.

Con linguaggi di scripting e software open source è possibile calcolare e visualizzare questo campo vettoriale come entità singola, non scissa in esposizione e pendenza. In questo post vediamo un esempio che si basa su Python per il processamento e su ParaView per la visualizzazione.

Dato che nella realtà le superfici topografiche non sono rappresentate da superfici analitiche ma da misure discrete nello spazio, una implementazione del calcolo del gradiente si deve basare su dati discreti, archiviati in formati come il raster a celle quadrate. Esistono varie formulazioni per il calcolo del gradiente nei GIS. Quella di base, implementata per esempio in Surfer, calcola i due gradienti lungo le direzioni x e y come differenze discrete tra i valori delle due celle contigue alla cella centrale.

dz/dx = [ Elev(i, j + 1) - Elev(i, j - 1) ] / 2*cell_size

dz/dy = [ Elev(i - 1, j) - Elev(i + 1, j) ] / 2*cell_size

In questo script Python, modificato da quello precedentemente creato per il calcolo della pendenza lungo direzioni variabili nello spazio, questa formulazione fornisce la possibilità di definire la distanza delle due celle laterali o verticali dalla cella centrale, in maniera tale da utilizzarlo come un fattore di smoothing che regola la risoluzione effettiva per la quale il gradiente è calcolato.

dz/dx = [ Elev(i, j + n) - Elev(i, j - n) ] / 2*n*cell_size

dz/dy = [ Elev(i - n, j) - Elev(i + n, j) ] / 2*n*cell_size

Il valore rappresenta il numero di celle a sinistra e a destra per la differenza lungo le x, o sopra e sotto per le y, utilizzate nel calcolo (esempio con smoothing factor uguale a 2 in figura sottostante).

Lo script Python calcola le componenti lungo i due assi a partire da un file ascii nel formato ESRI grid e le salva, assieme ai valori di elevazione, in un file ascii nel formato vtk (legacy vtk format). Col software ParaView è possibile visualizzare ed elaborare questi dati scalari e vettoriali. ParaView non è particolarmente intuitivo come utilizzo, ma con alcuni tentativi si riesce a capire quali sono le funzioni di base utili per visualizzare dati.

Usiamo dati liberi dello Shuttle Radar Topography Mission (SRTM2), con risoluzione nominale di 90 m, del lago di Bolsena, nel Lazio.

Il modulo del campo vettoriale viene calcolato automaticamente da ParaView, a partire dalle sue componenti.

Queste stesse componenti sono visualizzabili e rendono bene l'andamento morfologico della superficie, meno evidente dalle immagini satellitari, e che consiste in caldere di varie dimensioni, fra loro coalescenti. I valori positivi delle componenti indicano salite spostandosi da ovest verso est, nel caso del gradiente lungo le x, e dall'alto verso il basso, per il gradiente lungo le y.

Ed infine, il campo vettoriale può essere rappresentato utilizzando una serie di glifi che indicano l'orientazione dei vettori ed il loro modulo, sempre a partire dalle singole componenti cartesiane. Ovviamente anche nei GIS è possibile rappresentare questa informazione, ma occorre calcolare il valore angolare dell'orientazione ed utilizzarlo come campo di rotazione del glifo. Nella figura sottostante è rappresentato il campo vettoriale per un settore limitato della zona considerata.

Lunedì 28 giugno: Live broadcast (gratuito) della GEOSTAT 2010 summer school

Domattina inizia la GEOSTAT 2010 summer school a Plasencia, Spagna. Per chi volesse seguire via web le lezioni di geostatistici come Roger Bivand, Edzer Pebesma, Gerard Heuvelink, Olaf Conrad, Markus Metz, Victor Olaya e Tom Hengl, dovrebbe essere disponibile un canale web apposito, http://geostat2010.info/Live, dal quale seguire in diretta le lezioni.
Il numero massimo di connessioni in contemporanea sarebbe di soli 60.
Slide, etc? http://geostat2010.info/programme

La pendenza topografica lungo direzioni variabili

Nota: a Novembre 2011 ho pubblicato un post su un plugin Quantum GIS per il calcolo della pendenza, direzionale e massima: http://gisoftw.blogspot.com/2011/11/un-plugin-quantum-gis-per-il-calcolo.html


La pendenza topografica che viene calcolata classicamente nei GIS è la stima della inclinazione della superficie topografica locale, approssimata da un piano. Rappresenta cioè la massima inclinazione di questo piano, corrispondente al gradiente massimo di una superficie. Questo calcolo è permesso dai vari GIS che trattano dati DEM, come Saga, Grass, QuantumGis, ArcView 3 o ArcGis con le rispettive estensioni Spatial Analyst, etc.

Per applicazioni specifiche, possiamo avere bisogno di calcolare la pendenza direzionale, cioè il gradiente della superficie topografica lungo una particolare direzione, compresa tra 0° e 360° rispetto al top della mappa. A mia conoscenza, l'unico software GIS che consente questa operazione è Mirone. Grass consente di calcolare le derivate direzionali lungo le sole direzioni x e y, cioè a 90° e a 0° rispetto al top.

In rari casi, ci può essere necessario calcolare le pendenze di una superficie, topografica o meno, lungo delle direzioni che variano nello spazio. Per questi casi, non sono a conoscenza di alcun software GIS che implementi questa funzionalità. Nel caso specifico, questo tipo di calcolo è necessario per una tesi di dottorato condotta sulle relazioni delle megadune antartiche con i fattori ambientali che le regolano. Interessa calcolare la pendenza delle superficie topografica lungo le direzioni dei venti che variano nello spazio su distanze anche relativamente brevi (ovviamente). Risultati preliminari di questa ricerca sono stati presentati recentemente al congresso IPY di Oslo [1].

Un campo di megadune antartiche, visibile dal mosaico Ramp Radarsat.

Per implementare questa funzionalità, assieme a Maja Radivojevic, la dottoranda chiamata in causa, abbiamo scelto di usare Python. Non nell'ambiente ArcGis o QuantumGis, ma svincolati da ogni particolare software, commerciale o open source che sia.

Il primo quesito per creare questo tipo di algoritmo è: quale metodo di calcolo utilizzare? Ne esistono vari, di complessità variabile. Essendo una prima versione del programma, e essendo le superfici glaciali antartiche in generale piatte, abbiamo scelto una metodologia semplice, che è usata per esempio anche in Surfer e che trovate dettagliata in http://www.malg.eu/pendenzadirezionale.php

Un secondo quesito, quando si analizzano superfici raster, come i DEM, è: quale risoluzione utilizzare? Lo slope, come l'aspect e altri parametri topografici, dipendono dalla risoluzione spaziale usata nei calcoli. Cambiando risoluzione, varia anche il risultato perché questo è relativo alla scala di analisi. La pendenza di una catena montuosa calcolata con una risoluzione di 100 m sarà molto più variabile che se calcolata con una risoluzione di 10 km. Per evitare di dovere ricampionare il DEM di base prima di ogni analisi a differente scala spaziale, è stato considerato un fattore di “smoothing”, che esprime il numero di celle a sinistra e a destra, sopra e sotto, la cella per la quale si calcola lo slope direzionale. Un alto fattore di smoothing considera celle più distanti da quella centrale per i calcoli e quindi il risultato esprime una scala spaziale più ampia rispetto a quella espressa da bassi valori di smoothing.

Esempio di calcolo della pendenza usando uno smoothing factor di 2: vengono utilizzati i valori delle celle a due pixel di distanza in orizzontale e in verticlae dalla cella per la quale si effettua il calcolo.

Chiariti questi due aspetti, la creazione dello script è abbastanza semplice, grazie all'elevato livello di astrazione del linguaggio Python.

Lo script python può essere scaricato e utilizzato in locale da Python.

I dati di partenza richiesti sono:

a) grid delle elevazioni (DEM) in formato ESRI ascii

b) grid delle orientazioni in formato ESRI ascii

c) valore dello smoothing factor

 Il risultato consiste in un grid delle pendenze direzionali espressi in gradi (da 90° a -90°, positivo verso l'alto).

Quanto segue è un esempio di analisi e confronto con i risultati “classici”, su una porzione della superficie glaciale dell'Antartide orientale. Si può notare che la pendenza direzionale secondo direzioni parallele alle orientazioni locali dei venti produce dei risultati alle medie ed alte risoluzioni che differiscono da quelle del metodo classico della massima pendenza.

DEM della zona antartica analizzata.

Orientazioni dei venti utilizzata per il calcolo delle pendenze direzionali mobili.

Mappa delle pendenze direzionali calcolate secondo le orientazioni espresse nella mappa dei venti. Valori positivi indicano pendenze in salita. I valori sono simboleggiati per classi di deviazioni standard dal valore medio.

Risultato del calcolo della pendenza direzionale variabile (valori assoluti) comparato con quello classico ottenuto dai software GIS (massima pendenza). Si nota che mentre le strutture principali non vengono modificate , le strutture di media risoluzione e di dettaglio differiscono. 

Riferimenti bibliografici

[1] Maja Radivojevic, Massimo Frezzotti, Mauro Alberti, 2010. Morphology and Constraining Factors of Antarctic Megadunes. IPY Oslo Conference.

Vettori, non vettoriale

I vettoriali nei GIS sono pane quotidiano. Ma si riferiscono al formato degli oggetti geometrici rappresentati nei GIS e non hanno diretta relazione con i vettori usati in fisica, ingegneria e nelle altre scienze. Se si vogliono trattare i vettori nei GIS, di quali funzioni si dispongono? Di solito, tranne pochi casi, nessuna....

Divergenza, rotore, angoli tra vettori, prodotti scalari e vettoriali?
Possiamo calcolarli o usando delle funzioni da noi scritte, in Python, in Avenue, per esempio, oppure migrando temporaneamente i dati in software come MATLAB. Scomponiamo i vettori in componenti cartesiane conservandoli in due o tre campi della tabella degli attributi, oppure in due o tre grid se consideriamo campi continui, e quindi possiamo applicare le classiche operazioni scalari come addizione, sottrazione, o vettoriali come il prodotto scalare o vettoriale. Se il risultato corrisponde ad un vettore, potrà essere archiviato come prima nelle sue componenti e quando necessario trasformato nelle sue componenti polari.

Maggiori informazioni in questo pdf.

Cressie e Bivand al Spatial Statistics 2011 - 23 - 25 March 2011, Olanda

Su Linkedin viene segnalata la conferenza:

Spatial Statistics 2011, 23 - 25 March 2011, The Netherlands

 

Mapping Global Change - http://www.spatialstatisticsconference.com

Fra gli speaker invitati vi saranno Noel Cressie - autore del famoso Statistics for Spatial Data, Peter Atkinson, Yong Ge, Gilberto Camara, Martin Schlater, Roger Bivand - che si interessa di Geostatistica tramite R.

Come furono calcolati i portolani? Sul Washington Post

Link
http://www.washingtonpost.com/wp-dyn/content/article/2010/05/21/AR2010052104713.html?hpid%3Dtopnews&sub=AR

Verso il Lisp?

Interessante post di Paul Graham, segnalato in mailing list di Python: i linguaggi di programmazione, dal lontano 1958, starebbe tendendo lentamente verso il Lisp, uno degli ultimi casi Python.

http://www.paulgraham.com/icad.html