Sunday 18 July 2010

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.