A Linux tool for calculating local best-fit plane attitudes from geological traces

The topographic traces of geological surfaces store information on the 3D attitudes of the geological surfaces. Knowing the coordinates of the intersection points between a topographic surface and a geological surface, it is therefore possible to estimate the local attitudes of a geological surface.

 

The tools

 

A short description of geoSurfDEM: it is composed of two tools, IntersectDEM and the new BestFitGeoplanes.
The former allows to calculate the intersection points between a geological surfaces (stored in the VTK format) and a DEM. The latter, that will be described in the current post, allows to estimate the local attitudes given a set of 3D points, all deriving from the intersection between a single geological surface and a topographic surface. BestFitGeoplanes is developed using C++ and Fortran. The algorithm uses Singular Value Decomposition (SVD) in order to invert the local traces into local best-fit planes. 

The input data is constituted by a set of points (their x, y and z coordinates), all related to a single, continuous geological surface. How to process them in order to derive the local attitudes? Since the source points are unconnected points (not lines), all deriving from a single, continuous surface, the 2D space is discretized into a raster, using a user-defined cell size.
Based on the points falling into a cell, we have three possible cases:

1. no points falling into the cell;
2. one or two points falling into the cell;
3. three or more points in the cell.

In the first and second case, no attitude inversion via SVD is possible.
In the third case, when the points are not all collinear, a solution is provided by the SVD method. This solution is attributed to the grid cell.
The algorithm output will therefore consist in a gridded set of points for which the local attitudes have been inverted.
 
What is the difference with deriving the attitude via a spatial interpolation using for instance kriging? These interpolations implicitly assume a 2.5D surface and do not allow 3D surfaces. Geological surfaces, on the other hand, due to folding, can be 3D surfaces, i.e. with more then one point for each x-y position. Local inversion via SVD does not constrain the geological surfaces to be 2.5D,  but that may be locally modeled via a planar surface.

 

Compilation

 

This application is developed in Linux and is available at: https://github.com/mauroalberti/geoSurfDEM.
A makefile is available for compiling this tool in a Linux environment.
The compilation sequence is:
cd path/to/source/files
make
make clean

Lapack and BLAS libraries must be available. The makefile assumes that Lapack is available in usr/lib/lapack (with name lapack) and BLAS in usr/lib/libblas (name blas). Modify the makefile accordingly, to adapt to your settings. Alternatively, an example of commands to build it (always in a Linux environment and with same libraries settings) is in https://github.com/mauroalberti/geoSurfDEM/blob/master/BestFitGeoplanes/compile

 

Use

 

Having compiled the application, it is possible to run it as console application (Fig. 1). The only user interaction after launching the application is providing a parameter file name ("param.txt" in Fig. 1 example).

Fig. 1. Example of a run of the application.

 

The input parameter file is a text file that lists six pieces of information:

1. the path of text file storing the 3D coordinates (x, y and z) of the intersection points of a single, continuous, geological surface;
2. the number of header lines in the file referenced in point 1;
3. the path  of the text file in which the georeferenced results will be stored;
4. the path to the analysis report file;
5. the path of the output grid (in ESRI ASCII grid format) that will list the number of intersection points for each grid cell;
6. the output grid cell size;

 An example is available in https://github.com/mauroalberti/geoSurfDEM/tree/master/BestFitGeoplanes, in the param.txt file:

/home/mauro/Documents/Ricerca/Codice/Geostrutturale/geoSurfDEM/test_data/BestFitGeoplanes/inters_malpi_135_35.csv
1
/home/mauro/Documents/Ricerca/Codice/Geostrutturale/geoSurfDEM/test_data/BestFitGeoplanes/bfg_malpi_13535_100.txt
/home/mauro/Documents/Ricerca/Codice/Geostrutturale/geoSurfDEM/test_data/BestFitGeoplanes/rep_malpi_13535_100.txt
/home/mauro/Documents/Ricerca/Codice/Geostrutturale/geoSurfDEM/test_data/BestFitGeoplanes/pts_malpi_13535_100.asc
100 

The application output is represented by the three files listed in points 3-5 of the previous list. The more important is the inverted result file (point # 3), that consists in a csv file listing a few fields (see Fig. 2):

1. x: x coordinate of the cell grid center;
2. y: y coordinate of the cell grid center;
3. pt_num: number of points used for the inversion. Required minimum is 3.
4. dip_dir: inverted dip direction for geoplane points in cell;
5. dip_ang: inverted dip angle for geoplane points in cell;
6. x_range: spatial range along the x-direction (E-W) for the inverted points in the considered cell;
7. y_range: spatial range along the y-direction (N-S) for the inverted points in the considered cell;
8. z_range: spatial range along the z-direction (vertical) for the inverted points in the considered cell;
9. pseudo-volume: "pseudo"-volume defined by the inverted points in the cell, given by the product of the three previous ranges (i.e., x_range, y_range and z_range) - possibly removed in successive releases.

Fig. 2. Tabular view of the content of the inversion result file, as viewed by importing the file in QGIS.

 

A theoretical case study 

 

This post presents a practical assessment of the result, by using a theoretical surface for which the expected local attitudes are known (currently, as of February 2017, all example data are available at https://github.com/mauroalberti/geoSurfDEM/tree/master/test_data).
A theoretical plane has been generated using simSurf, with dip direction equal to 135° and dip angle of 35°.
This plane, saved in VTK format from within simSurf, has been used as input plane to be intersecated with a natural topographic surface, of the Mt. Alpi zone (Lucania, Southern Italy), using the IntersectDEM of geoSurfDEM.
The resulting intersection points are saved as a csv file, that can be used as input for the BestFitGeoplanes. The intersection points, to be used for the best-fit-plane local inversions, are represented in Fig. 3.

Fig. 3. Input points (yellow) representing the theoretical intersection between a geoplane oriented 135°/35° (dip dir. and dip angle) and a DEM topography (Mt. Alpi zone, Lucania, Southern Italy). Visualization using QGIS.

 

Using BestFitGeoplanes with the previously described param.txt file, we obtain a result that, for the georeferenced point part, imported in QGIS is shown in Fig. 4.

Fig. 4. The same as in Fig. 3, plus superposed the gridded points with inverted results (orange dots). Visualization with QGIS.

 

How much the inverted results conform to the theoretical input source, that as said is a plane with a dip direction of 135° and a dip angle of 35°?
A stereoplot, created with the geocouche plugin for QGIS, illustrates the degree of  concordance between the source attitude (135°/35°, blue great circle in Fig. 5) and the inverted local attitudes (semi-transparent orange great circles in Fig. 5).
The majority of inverted data conform closely to the expected result, while a few inversions show a minor deviation from the expected result.

 

Fig. 5. Stereonet representing the inferred local plane attitudes as semi-opaque orange great circles, and the source geological plane attitude (135°/35°) as blue great circle. Created with geocouche.

The deviations of the inverted results from the expected value were calculated using the "Geological angles" of geocouche and the statistics calculated with QGIS (see Fig. 6). The maximum is 9.3° and the minimum almost zero, while the median and mean deviation values are lower than 0.5°. The standard deviation is about 1.2°. So in general we can be quite confident in the generated results. More in-depth analyses of the deviations of expected-versus-inferred results could be the subject of a still to-be-written paper.

Fig. 6. Statistics for angular deviations of the calculated results from the theoretical test case.

Edits

2022-12-29: improved paragraph styling; modified input parameters description (items 5 and 6)

qgSurf now supports QGis on-the-fly projection

qgSurf is a plugin for QGis, that allows to calculate the best-fit-plane given a set of point with DEM-derived heights, or the intersection between a geological plane and a DEM. Recently it has added some basic tools for the creation of 3D analytical surface and their deformation, even if they are still at an experimental stage.
The new version, 0.3.1, adds support for on-the-fly projection, so that, for instance, it is easier to use satellite data providers as Bing or Google for geological calculations.

When using the Best-Fit-Plane tool, DEMs can be in lat-long, but in that case the project must be set to a projected, planar CRS. DEM elevations and projected planar distances must be in the same measure unit (e.g., meters).
On the other hand, the DEM-plane intersection tool does not work properly with DEM in lat-long, since it requires that the horizontal distance measure unit to be the same as the vertical one, impossible for DEMs with data in lat-long. Also mixing meters for horizontal distances with feets for elevations will produce erroneous results. It is safe however to have a DEM in UTM 32 (with x, y and z values in meters) and the project CRS in Lambert Conformal Conic.

The plugin can be installed from the QGis plugin manager, or downloaded from http://plugins.qgis.org/plugins/qgSurf/ or https://bitbucket.org/mauroalberti/qgsurf/downloads.
It works for QGis >= 2.0, and has been tested in Windows Vista and Ubuntu LTS 10.2.
Errors in Linux or Windows for the 3D geosurface simulation tools could possibly be related to an old installed Matplotlib version, and solved by updating this last.
 

Best-fit geological plane inferred from DEM and georeferenced points along traces

We want to calculate the plane attitude that best fits a set of points located on a DEM. These points could follow a geologic trace, for instance a stratification or a fault that we or others have mapped in the field. In other cases, we could clearly see lineaments in a satellite or aerial image, and we want to extract the geological information exposed by these lineaments.
qgSurf, a plugin for QuantumGIS, an open-source and free GIS software, adds this functionality in its new version (0.2.1, installable from QGIS plugin manager or downloadable from QGis plugin repository). Previously its aim was to calculate the intersections between a geological plane and a DEM. The calculation of the best-fit plane given points on a DEM is a somewhat inverse procedure so that they complement each other.
To see a short help for the plugin, look at: https://bitbucket.org/mauroalberti/qgsurf/wiki/Help

The basis of the algorithm is the application of singular value decomposition to derive the eigenvectors of a set of measures. See for instance the discussion in Best fit plane algorithms why different results?
Obviously we have to dispose of a DEM, e.g., Aster or SRTM or others.
A very interesting QGis plugin, OpenLayers, allows to load GoogleEarth or Bing maps in a project, in order to use them as backgrounds for our analysis. However, this requires the on-the-fly reprojection of layers while qgSurf does not support it for DEMs, so we need to set as the project projection that of the DEM: UTM, Lambert, or whatsoever.

The processing sequence is:

a) load in the QGis project the required DEM(s) layers and whatsoever vector or image layers needed for your analysis

b) use "Get current raster layers" in qgSurf plugin: this will allow the plugin to know which raster layers are currently loaded

c) from the "Use DEM" combo box, choose the required DEM and make sure that the QGis project and the DEM have the same projection

d) from the "Best-fit-plane calculation", press "Define points in map": this will allow you to define in the canvas at least three, and possibly more, points, whose coordinates will be listed in the plugin widget.

e) with at least three points defined, you can calculate the best-fit plane by pressing "Calculate best-fit-plane": a message box will report the dip direction and dip angle of the calculated plane

f) you can add even more points and again calculate the best-fit plane; otherwise, if you want to start a new analysis on the same DEM, go to e), or if you want to use another DEM, go to c) if it is already loaded in the project, or load it in the project and then go to b)

Example

Mt. Raparo, Basilicata, limestones of the Panormide Complex.

Source DEM is Aster, projected in WGS84-UTM33N.  Base map is Bing, loaded in QGis via openLayers plugin. Project CRS is set to the same of the source DEM, i.e. WGS84-UTM33N, otherwise the plugin wouldn't work.

Two level are recognizable in the aerial image, and a few points are drawn to derive their attitude using the new plugin functionality. The results are 336.2 / 8.7 and 337.8 / 13.0 (dip direction nd dip angle).

The results are checked with field measurement made by myself and M.C. Lapenta during field work in 1991: 320 / 5, 290 / 15 and 345 / 15. The geological plane attitudes are congruent both with themselves and with field measurements (note that while the field measurements refer to magnetic North, GIS-derived values refer to map top for the current projection).

Ricuperare giaciture geologiche da immagini satellitari e DEM: nuova versione del plugin qgSurf per QGis

Presento una nuova versione sperimentale di qgSurf, la 0.2.0, che ha come intento quello di permettere, all'interno di Quantum Gis, di utilizzare immagini satellitari, come quelle liberamente disponibili in GoogleEarth, accompagnandole a DEM, per ricavare in maniera interattiva, per trial-and-error, la giacitura di piani geologici .

Schermata del plugin, con un esempio di analisi su M.te Alpi (Basilicata), da immagine Google, e dati DEM Aster a 30 m di risoluzione.

Come detto la versione è sperimentale (ogni segnalazione di baco è benvenuta -> alberti.m65@gmail.com) ed è stata testata fondamentalmente sulla versione Windows di Quantum GIS 1.8.0. A differenza della versione precedente, l'interazione con la mappa avviene all'interno del canvas di Quantum GIS, il che permette di utilizzare come base cartografica qualsiasi livello vettoriale o raster, come appunto le immagini satellitari disponibili dai servizi di GoogleEarth, Bing, etc, caricabili in Quantum GIS attraverso il plugin OpenLayers.

Finestra del plugin qgSurf v. 0.2.0, tab 'Geographic data'.

Brevemente, la successione di utilizzo di qgSurgf v. 0.2.0 è la seguente:
0) occorre settare come sistema proiettivo del progetto QGis quello del DEM utilizzato per la determinazione;
1) dal tab 'Geographic data' si definisce il DEM sul quale viene calcolata l'intersezione del piano. Questo permette quindi di definire nella mappa il punto sorgente (source point) che contiene il piano geologico;
2) dal tab 'Geological data' si settano i parametri (immersione ed inclinazione) del piano, e col bottone 'Calculate intersection' vengono calcolati i punti di intersezione teorici tra DEM e piano geologico. Questi parametri possono essere variati sino a trovare un 'best-fit' con un lineamento (stratificazione, faglia, etc.) visibile in immagine satellitare.
3) nel tab 'Output' troviamo gli strumenti per salvare le intersezioni, come punti o linee, in shapefile.

Ovviamente se si dispone di misurazioni strutturali di terreno è possibile testare la loro corrispondenza con gli andamenti visibili in immagini satellitari/aeree o in cartografie geologiche georiferite.

Il plugin può essere scaricato direttamente all'interno di Quantum GIS, oppure da http://plugins.qgis.org/plugins/qgSurf/ per essere poi decompresso nella cartella dei plugin Python di Quantum GIS.

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 (x₀, y₀, z₀) e perpendicolare ad una retta con coseni direttori (l, m, n) è (Zwirner, 1983, p. 198):
l (x – x₀) + m ( y – y₀ ) + n ( z – z₀ ) = 0
sviluppabile come:
z = - ( l / n ) (x – x₀ ) – ( m / n ) ( y – y₀ ) + z₀

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 – x₀ ) + b ( y – y₀ ) + z₀
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 ( x_n, z_n ) e ( x_n+1, z_n+1 ).

Se i segmenti relativi al DEM e al piano lungo l'asse delle x hanno coefficienti angolari ed intercette ( m₁, q₁ ) e ( m₂, q₂ ), l'intersezione ha valore:
x_inters = ( q₂ – q₁ ) / ( m₁ – m₂ )

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.

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.

Mirone e la visualizzazione di DEM

sito: http://w3.ualg.pt/~jluis/mirone/index.htm
creatore: Joaquim Luis

Mirone è un software free (coffee-right, per la precisione) che presenta delle funzioni molto interessanti per il processamento di DEM (e non solo, come vedremo) accoppiato ad una notevole facilità d'uso. Questo grazie, oltre alle capacità del suo programmatore, Joacquim Luis, anche alla scelta di basarlo su un ottimo software scientifico commerciale, MatLab. Ma non c'è da preoccuparsi di dover reperire MatLab, se non disponiamo di una regolare licenza. Mirone infatti è disponibile anche in una versione stand-alone, che incorpora già legalmente le librerie necessarie di MatLab.

Mirone è un programma eclettico, sui-generis, a cavallo tra la geofisica e il gis: per darvi un'idea, presenta funzioni che vanno dal calcolo del tempo di percorso degli tsunami (mai provato) al calcolatore del movimento delle placche (dato solo un'occhiatina veloce), a funzioni gis avanzate per i raster sino a fornire delle stupende e ricchissime colormap per i raster. Al confronto, sia per le colormap sia per determinate funzioni di elaborazione dei raster, un noto software commerciale zoppica molto. Un software che invece presenta funzioni analoghe, GMT, è nettamente più complicato come uso. Come aspetto negativo, le capacità di Mirone nel gestire i dati vettoriali sono rudimentali: si limitano all'import e loro visualizzazione e pochissimo più. Anche capacità di gestione dei sistemi proiettivi sono marginali nella filosofia del programma, anche se forse verranno ampliati nelle successive versioni.

In questo post tratteremo di come Mirone permette di visualizzare DEM. Useremo la versione 1.4 di questo software (la più recente è la 1.5).

Innanzitutto, quali sono i tipi di dati leggibili e in quale sistema? Partiamo dal sistema proiettivo. Geografico o cartesiano [File - Preferences]. Non è una ampia scelta, perché il cartesiano a questo punto raccoglie tutti i dati proiettati. Ma infatti non è questo il punto di forza di Mirone.
Quali tipi di dato possiamo importare? Tutti i più diffusi tipi grid e raster, dai grid di Arc/Info a quelli di Surfer [File - Open Grid/Image]

Una volta importati i dati, siamo pronti per lavorare....
Vedremo brevi esempi di tre tipi di usi: la scelta delle colormap, lo shading, i profili dinamici e la creazione di mappe delle derivate direzionali.

Le colormap di Mirone: esempi di visualizzazione usando dati SRTM a 90 m di risoluzione

Gli stessi dati, che si riferiscono alla zona ad Est di Foligno-Spoleto, sono rappresentati con differenti colormap.
Come si nota, a seconda della scelta della colormap influenza forle strutture che vengono enfatizzate

ml - hsv



gmt -split

gmt - dem_screen

gimp horizon_1

Applicazione dello shading

Scelta del tipo di shading per una zona dell'Umbria - dati SRTM 90 m

In alto ombreggiato, in basso per confronto senza ombreggiatura.