Showing posts with label ParaView. Show all posts
Showing posts with label ParaView. Show all posts

Thursday, 17 March 2011

Paraview, un software open source per la visualizzazione scientifica


Visualizzazione ed analisi di dati con geometrie 3D (o 4D considerando anche la variabile tempo) possono essere effettuate con software open source di visualizzazione scientifica come Paraview, VTK, VisIt e Mayavi2. Questi software consentono inoltre di processare variabili di tipo vettoriale e tensoriale, possibilità mancante nei GIS classici. In questo post accenno ad alcuni aspetti di elaborazione e visualizzazione con Paraview.

Paraview è uno dei software scientifici open source sviluppati da Kitware, assieme a VTK, ITK e CMAKE e altri. E' un'applicazione che può lavorare su un singolo processore o su più computer in parallelo. Ne è stata rilasciata recentemente la versione 3.10.0, il 9 marzo. Presenta anche un calcolatore ed una shell in Python (vers. 2.7), che consentono di estendere le funzionalità base.


Formato dati in input

Anche se può leggere decine di differenti formati dati, sono due i formati standard di import dati: VTK classico ed il più recente formato XML, creato per consentire una maggiore varietà di formati geometrici in input. Il tipo di strutturazione geometrica dei dati, come anche il loro assemblaggio nei file sono più ricchi e flessibili che nei GIS. Nei dataset GIS classici il singolo file (o set collegato di file, per esempio uno shapefile) mantiene una sola tipologia di formato geometrico (raster o vettoriale), accompagnata da informazioni descrittive sotto forma di variabili categoriali e scalari. In un file VTK o XML possono essere conservati contemporaneamente più tipi di geometrie (raster 2-3D e vettoriali) e più variabili (che siano scalari, vettoriali o tensoriali).

I dati di input possono essere creati da software GIS esportandoli in formato VTK. GRASS consente di esportare dati scalari come file VTK. Script e programmi in linguaggi come Python, C++, Fortran permettono ugualmente di convertire dati GIS anche come variabili vettoriali e tensoriali, sfruttando più livelli oppure le informazioni contenute in più campi di tabelle degli attributi.


Modalità di processamento di dati scalari e vettoriali 2D

La tecnica di elaborazione dei dati è molto simile a quella usata nei GIS. Si utilizzano filtri che possono essere applicati in cascata ad un dataset caricato. Come si vede dalla Fig. 1, sono disponibili decine di filtri, applicabili a diverse tipologie di dati, per una loro elaborazione oppure visualizzazione (glifi, linee di flusso, etc). Il risultato dell'applicazione di un filtro può essere salvato come dataset a se stante (File → Save Data, avendo selezionato il dataset da esportare).

Fig. 1. Elenco alfabetico dei filtri in Paraview 3.10.0. In nero quelli applicabili sul tipo di dati nel livello selezionato, in grigio i filtri non disponibili per lo specifico tipo di dato selezionato.

Come esempio concreto, consideriamo un dataset comprendente sia variabili scalari sia vettoriali, con zone che contengono valori nulli (zone con colore blu in Fig. 2a). A questo può essere inizialmente applicato un filtro di threshold (Filters → Common → Threshold) che limita i valori accettabili di una particolare variabile scalare, in maniera tale da scartare i valori che rappresentano misure mancanti. Nel dataset filtrato saranno quindi escluse le localizzazioni con valori della variabile fuori range (Fig. 2b). A questo dataset filtrato può essere applicato un ulteriore filtro, anche su una variabile differente da quella precedentemente considerata. Per esempio, possiamo applicare un filtro di streamline (Filters → Common → Stream Tracer)ad una variabile vettoriale che rappresenta velocità di spostamento, che produrrà le linee di flusso (streamlines) del campo vettoriale (Fig. 2c).

Fig. 2a.
Fig. 2b.
Fig. 2c.

Una volta caricati e processati una serie di dati, può essere utile salvare lo “stato” della visualizzazione (File → Save State) per un successivo riutilizzo (File → Load State), come avviene nei GIS con il salvataggio in progetti.


Esempi di visualizzazione di dati 3D

Consideriamo gli eventi sismici registrati durante l'evento di Colfiorito del 1997 (terremoto dell'Umbria-Marche), derivanti da un dataset creato dall'INGV e distribuito in occasione di un convegno per il decennale del terremoto (Roma, 2007).

I dati disponibili comprendono gli ipocentri, conservati in uno shapefile puntuale. Questo livello è stato importato in GRASS e da qui esportato come file VTK. Per poter visualizzare anche la densità degli eventi, è stato scritto un programma in C++ per il calcolo della densità di punti in 3D a partire da dati testuali x-y-z, con creazione di un raster 3D in formato VTK. L'algoritmo utilizza un kernel gaussiano per calcolare le densità di punti. Un'introduzione teorica su questo soggetto si trova in http://www.spatialanalysisonline.com/output/html/Pointdensity.html

Una visualizzazione iniziale può rappresentare congiuntamente le informazioni puntuali (localizzazioni degli ipocentri) sia raster 3D (densità degli ipocentri), come in Fig. 3.

Fig. 3. Visualizzazione degli ipocentri come punti accompagnata dalle densità più elevate degli eventi sismici. Sequenza dell'Umbria-Marche del 1997, dati INGV.

Dati “raster” 3D possono essere visualizzati con varie tecniche, fra cui le sezioni (slices) e le isosuperfici (isosurfaces).
Sezioni (slices) attraverso un raster 3D (Filters → Common → Slice, con livello selezionato) sono utili per comprendere la struttura tridimensionale di un fenomeno. In Fig. 4 uno slice orizzontale e altri verticali fra loro paralleli evidenziano come le densità massime degli eventi sismici delineano superfici che immergono verso sinistra (Ovest), con una orientazione circa Nord-Sud (solo accennata nella sezione orizzontale presentata): Queste orientazioni sono congruenti con le giaciture delle faglie normali nella zona esaminata.

Fig. 4. Slices delle densità di ipocentri ottenuta per la sequenza sismica dell'Umbria-Marche del 1997. Dati INGV.

Una isosuperficie collega le posizioni con uguale valore della variabile considerata. Più isosuperfici possono essere visualizzate contemporaneamente. In Paraview vengono create col filtro Contour (Filters → Common → Contour). In Fig. 5 si nota nella zona settentrionale un allineamento circa N-S di valori elevati di densità di ipocentri, mentre nella porzione meridionale i massimi di raggruppamento degli ipocentri sono più limitati e discontinui.

Fig. 5. Gli ipocentri rappresentati come punti bianchi, mentre i volumi in rosso rappresentano le isosuperfici di densità uguale alla metà del range dei suoi valori (13 in un range da 0 a 26).

Alcuni svantaggi e problematiche di Paraview

Paraview è un software ricco di funzionalità e non immediato da apprendere. Presenta anche una serie di svantaggi e problemi, che presumibilmente si ridurranno man mano che questo software matura:
  • presenza di tutorial e documentazione libera ridotta rispetto ad altri software, comunque in via di miglioramento;
  • l'ordine di visualizzazione dei livelli, banale da definire nei GIS, è apparentemente non gestibile in Paraview, per cui in alcuni casi occorre testare la visualizzazione in più schermate differenti;
  • episodici mancati refresh della visualizzazione quando si zooma o si naviga nella schermata; riducendo ad icona il programma e visualizzandolo nuovamente la schermata si aggiorna correttamente;
  • nel caso del filtro Stream Tracer, si è riscontrato in vari casi che nessun risultato viene prodotto, senza una evidente regola.


Link per approfondimenti



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.