Saturday 26 October 2013

Plotting geological planes in profiles with Quantum GIS

Plotting the traces of geological planes in a vertical profile is a task that can be automated using a new tool of the qProf plugin for Quantum GIS. The geological measure projection on the section is determined as the intersection line between the geological plane and the section plane. The intersection position in the profile is considered to be the nearest point on the intersection line, within the geological plane and starting  from the geological measure. This technique can be used when dealing with monoclinal successions, with no or little folding, while it is not applicable to folded successions, where the projection should follow the hinge line of the folds (cf. Groshong, 3-D Structural Geology).

The qProf plugin, developed by myself and Marco Zanieri, is available for both the new Quantum GIS 2.0 version, and the previous one, 1.8. It can be directly installed from the plugin manager of Quantum GIS, or downloaded from the Quantum GIS plugin repository (vers. 0.2.1 for QGis 2.0 and vers. 0.2.0 for QGis 1.8) and then unzipped into the Python plugin folder for QGis (in my case, in Windows Vista and for QGis 2.0 is: C:\Users\mauro\.qgis2\python\plugins).

As an analysis example I consider the Mt. Alpi zone (Lucania, southern Apennines, Italy), made up of a Mesozoic carbonatic platform sequence, plus upper Miocene carbonatic-clastic sediments, presumably pertaining to the Apulian platform. The structure in the profile area is almost monoclinal, with some minor faults dissecting the sediments (not represented in the figure below). The geological stratification measures (dip direction, dip angle) derive from my undergraduate thesis, while the DEM is from Aster data.



First I select all the stratification measures in the point layer (in order to use all of them in the resulting section), then I launch the qProf plugin. In the 'Profile' tab, I calculate a profile by using the Aster data and a line layer representing the vertical profile trace. An important requirement for using the geological cross-section tool is that the profile line should consists of only two points, the start and end points, since the current version is not able to handle profiles with changing orientations. This is a requirement only for using the geological section tool, while for just creating topographic profiles there is no such restriction.



Afterwards I switch to the 'Geological cross-sections' tab, I set the source layer for the geological measures, as well as an id field (may be text or integer), and the fields for the dip direction (0-360°, clockwise from North) and the dip angle (0-90°) values. The selected points in the structural measure layer will be projected in the resulting profile.




A profile with the projected measure traces is created, with the height-length ratio scale automatically set to 1:1. You can see that the stratification is almost constantly dipping at low angles towards South-East (to the right in figure below) and that in the southern section (rightmost) the topographic profile and the stratification traces are almost parallel, as you would know as you ever climbed this 'pettata'.




The results can then be exported in csv or shapefile formats from the 'Export cross-section results' section. The result file will contain all relevant information, including the Id values that aid in finding the position of each source measure in the created profile.



Monday 5 August 2013

A method to interpolate 'cylindrical' geological surfaces

I present here a method, implemented in Python 2.7, that allows to interpolate a portion of a geological surface (e.g., horizons, faults), given its topographic trace and the inferred attitudes at the trace start and end points. The method is based on the assumption of a 'cylindrical' geometry of the analysed surface.
'Cylidrical' is in the structural geologic meaning, i.e. by citing Ramsay & Huber, 1987, vol. II, p. 311: "If at all locations on the fold surface a direction can be found which lies parallel to the fold hinge line, the fold is termed cylindrical, and the line is known as the fold axis."
Substitute 'fold' with 'surface' in the previous sentence, for here we consider a general case of geological surface (stratifications, but also faults, and so on).

If a geological surface portion can be considered cylindrical, we can imagine to rotate it so that the surface axis (the 'fold axis' in previous R&H citation) becomes vertical. We consider this rotated axis to be z'. In this way, we can project all the topographic trace points on the horizontal plane x'-y'. These projected points can then be interpolated by a 2D curve using standard 2D spline interpolation.





The interpolated 2D analytical result can be transformed into a set of 3D, vertical segments having their lower and upper extremes located at z' minimum and maximum values, respectively. These lines can be back-rotated to the original basis frame, so that they represent the surface portion to be modelled.

The implementation in Python 2.7.3, named 'gSurf.Interp', requires the following modules:
- numpy
- scipy
- PyQt4
- guiqwt
- ogr



The source data must be in a comma-separated csv format, storing x, y and z values. In this example, the input data were produced using the QuantumGIS plugin 'qProf', that extracts elevations along profiles, given one or more DEMs (see previous post, linked at the bottom). From this source file, it is possible to define which fields stores the x, y, and z values (1, 2 and 4 in the figure, and in the accompanying file, see bottom)

An important step is the choice of the spline curve for the approximation of the rotated data.
This is done through an interactive graphic (based on the guiqwt module) showing the rotated data and the interpolation spline, whose adherence to the data is controlled by a single numeric parameter, the 'spline value' in the figure below.



Interpolated data are then back-rotated by the tool, so that we dispose of the 3D representation of the surface portion. The results are saved as a 3D line shapefile (see Fig. below).

Interpolation results, represented by a series of 3D lines (in red). Fault to be modelled is in green. Mount Alpi zone (Lucania, Southern Italy). Figure created with ArcScene (ESRI).

The Python source code, together with the example data presented in this post, are available in this page .

You can read also the previous post 'Details of a fault surface', on the same subject and concerning the same data set, but at a preliminary stage of creation.


Sunday 28 July 2013

Details of a fault surface

An investigation of the 3D geometrical details of a fault trace from Mt. Alpi area (Lucania, Southern Italy) is here presented, based on some QuantumGIS plugins and an in-progress Python tool. The southern portion of a well-exposed fault, separating Mt. Alpi Unit carbonatic rocks from Frido Unit metasediments, was mapped in GIS and its possible orientations were defined using qgSurf, in combination with qProf and a still-in-development Python tool that rotates 3D fault traces. Gnuplot and R were used to visualise in 3D the fault trace points.

Analysed fault segment in Mt. Alpi area (Lucania, Southern Italy). The fault segment (red line) separates the Mt. Alpi Unit, made up of carbonatic rocks, from the metasediments of the Frido Unit. Basemap: GoogleEarth; QuantumGIS with OpenLayers plugin used for map visualization.

qgSurf, a QuantumGIS plugin, helps in determining the possible planar orientations of geological surfaces, f.i. fault traces, given a topographic DEM and satellite images and/or geological information. However, experimentations with this tool have evidenced that the calculated orientations can be sensitive to the particular examined segment, particularly when the surface is not planar. More reliable results can be obtained when fault traces are visualised and examined in 3D, in order to preliminary isolate possible planar segments making up a complex fault surfaces.

qProf, another QuantumGIS plugin, has been used to extract the elevation information of a profile from the used DEM, a 30 m Aster. Elevation info were read from qProf output and saved as an xyz text file using the in-development Python tool, so that they could be easely read by gnuplot. With gnuplot it is possible to produce interactive 2D and 3D graphics.


Gnuplot-generated 3D graphic of the xyz data of the fault trace in the Mt. Alpi sector, extracted from an Aster DEM using the qProf plugin.


Trough tentative rotations, 3D visualizations with gnuplot allowed to recognise a probable planar fault section almost 1-km long in the southern section, while the northern section is more complex (see figure below). This complexity could arise from the presence of a set of minor faults, orthogonal to the analysed fault, that are well mapped in the Mt. Alpi Unit, but with unknown prosecutions in the Frido Unit, due to the low relief and reduced exposures of this unit rocks.

.
Analysis of  Mt. Alpi fault trace, using gnuplot to visualize in 3D the topographic trace of the fault. Note that in gnuplot  graphics, the horizontal axes are not in 1:1 scale. Fault traces/points are in red in both GIS map (top) and gnoplot graphics (middle and bottom). In the middle graphics the view is top-to-down, while in the bottom graphics the view is rotated in order to evidence a southern, almost 1-km long planar segment.
Graphic created with gnuplot, gimp and Inkscape.


Using qgSurf on the previously defined fault sections, produced a value of 143.6/18.1 (dip direction/dip angle) for the sourthern planar section, and a mean value of 074.4/28.1 for the northern, more complex section.

Stereonet of the two derived fault sections. Created with Stereonet by Allmendinger.

The still in-development Python tool has the aim of interpolating 3D fault segments from a topographic trace and the local fault attitude a the two extremes of the analysed fault segment.
The tool rotates the fault traces (as xyz point values) in order to plot the trace points along a plane (x'-y') that is perpendicular to the intersection of the fault attitudes at the two extreme points (z' axis), and that contains the start point. In this way, assuming that the fault attitude changes are (almost) cilindrical, the 3D interpolation problem is reduced to a 2D one, that could be solved by user-chosen interpolation parameters.


Rotated basis frame for the fault trace, given its start and end points, coupled with inferred fault attitude at the extremes. The frame is centered on the trace start point. The z' axiz is parallel to the intersection of the two inferred local fault attitudes. The x' axiz is located in the plane containing the z' axiz and the trace end point. The y' axis is, obviously, perpendicular to both x' and z' axes. Fault traces are plotted on the x'-y' plane. 

The result for the fault trace is represented in the figure below. The southern section is well approximated by the 143.6/18.1 (dip direction/dip angle) solution (green line) while the solution for the northern section (red line) is clearly only an approximation that should be improved.


Fault trace points plotted on the x' (horizontal) - y' (vertical) plane. The green line indicates the southern section fault solution, the red one the northern fault solution. Graphic created with gnuplot, gimp and Inkscape.


It planned for the in-development tools to allow the user to interpolate the rotated fault traces in 2D by means of splines, in order to reproduce the local structures of the fault, and then expand the 2D solution to 3D assuming a "cilindrical" fault surface, and meshing it to allow the import in 3D GIS or visualization programs.





Tuesday 25 June 2013

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









Sunday 26 May 2013

Pierre Raybaut rilascia WinPython, ambiente portabile per Python 2.7 e 3.3

Raybaut è lo sviluppatore di Python(x,y) e di packages come guiqwt, guidata, spyder e userconfig, tutti rigorosamente OS.
In questi giorni ha rilasciato una nuova versione di WinPython, ambiente alternativo a Python(x,y) in quanto è portabile, ma conserva l'orientamento per l'analisi e la gestione di dati scientifici.
Presenta infatti packages come NumPy, SciPy, Matplotlib e PyQwt, PyQt, numexpr, Pandas, scikit-image, scikit-learn, sympy, tables, xlrd e xlwt ed anche Cython.
E' possibile utilizzare IPython ed IPython Notebook, interessante per testare e divulgare codice Python.

La sua portabilità consente di avere installazioni plurime di Python che non interferiscono, p.e. è possibile mantenere Python 2.7 per la produzione ed avere a disposizione  Python 3.3 per testare la nuova serie.

Link:
https://code.google.com/p/winpython/


Friday 3 May 2013

Plotting profiles from GPX files within Quantum GIS

A new version of qProf, a Quantum GIS plugin, allows the plotting of profiles by reading the track points stored inside a GPX file. This new tool complements the previous implementation of profile calculation (elevation and slope) from DEMs and line tracks.

Screenshot for a GPX profile produced by qProf plugin v. 0.1.4.
Upper window: elevations, lower w.: slopes.
Data (implicit) courtesy of S. Peduzzi.

Data in source GPX file should contain point tracks, with point coordinates in latitude and longitude, plus elevation and time. Elevations are generally referred to the WGS84 datum, not to the orthometric height.

The algorithm converts the latitude-longitude coordinates in ECEF Cartesian coordinates, avoiding the need for any planar, 2D projection. The 3D distances between consecutive track points are easily calculated, and with some trigonometry formulas, slope and 2D horizontal distances between points are then derived by assuming a vertical offset between consecutive points equal to the elevation differences derived from the GPX file.

The resulting profile can be saved as .csv files, 2D point shapefiles, or as 3D line shapefiles.

This plugin can be easily installed using the QGis plugin manager, or by downloading it from the QGis plugin repository and extracting it into your personal Python plugin folder.


Wednesday 3 April 2013

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.



Saturday 16 March 2013

qProf: un nuovo plugin QuantumGIS per la creazione di profili topografici


Questo è il primo rilascio di QProf, un plugin per la creazione di  profili di elevazione e di pendenza 2D e 3D, a partire da uno o più DEM. E' possibile visualizzare più profili contemporaneamente. I risultati possono essere salvati come dati (shapefiles puntuali e lineari, csv) o immagini. L'idea iniziale nasce da Marco Zanieri e l'implementazione è di Mauro Alberti. Il plugin è installabile direttamente come nuovo plugin in Quantum GIS.

Esempio di profili prodotti da qProf da DEM SRTM2 ed Aster di una zona della Valnerina (Umbria): in alto i profili topografici,  in basso il corrispondente profilo delle pendenze lungo la traccia: valori positivi indicano salite, negativi discese.

Tool analoghi in Quantum GIS

In Quantum GIS sono disponibili almeno altri due moduli che consentono l'elaborazione di profili topografici: Profile Tool di Boris Jurgiel, Etienne Tourigny e Patrice Verchere, e Profile from Line di Ricardo Garcia Silva. Il primo utilizza come traccia di input un percorso "dinamico", definito dall'utente al momento dell'analisi, su uno o più DEM. L'output grafico e tabellare riguarda le elevazioni ma non le pendenze lungo il profilo. Il secondo determina le elevazioni, a spaziatura costante, a partire da un livello lineare e produce come output un livello puntuale con le elevazioni determinate a partire dallo specifico DEM scelto dall'utente. Non viene creato un output grafico. Entrambi i plugin estraggono le elevazioni dai DEM per nearest neighbor.


Descrizione del plugin

Schermata del plugin.


1 - Path shapefile
La traccia del profilo è costituita da una linea singola e continua  (anche se con andamento variabile), contenuta in uno shapefile.

2 -Use DEMs
Come fonte delle elevazioni possono essere scelti i raster DEM che sono già caricati nel progetto QGIS corrente. Ovviamente in ogni raster ci deve essere una sola banda di dati; in caso contrario il modulo fornirà un messaggio di errore. Raster di notevoli dimensioni possono generare errori di memoria del modulo (MemoryError).

3 - Max sampling distance
L'utente definisce la distanza massima di campionamento lungo la traccia del profilo. Questo plugin mantiene tutti i punti di campionamento definiti nel profilo, sia l'inizio, sia quelli interni, sia quello terminale. La spaziatura risultante può essere quindi più fitta, oltre che irregolare,  della max sampling distance.
Il passo di campionamento determinerà la risoluzione spaziale del profilo, incidendo anche, ovviamente sulla distanza cumulativa 3D. Quest'ultima sarà tanto maggiore quanto minore la distanza (fino al limite di risoluzione del DEM), dato che la topografia ha carattere frattale (e.g. Turcotte, 1997).
Si suggerisce di impostare un valore circa analogo a quello della risoluzione dei DEM utilizzati: per esempio per un DEM ASTER si può inserire un valore di 30 m, per un DEM SRTM2 a 90 m un valore di 90 o superiore.

4 - Calculate profile
I dati DEM vengono letti da Python tramite il binding a GDAL. I valori di elevazione per i punti di campionamento sono interpolati per bilinear interpolation a partire dai centri delle quattro celle che racchiudono il punto da interpolare.
Quando viene modificato un parametro di analisi, come p.e. i DEM usati per il profilo, occorre ricalcolare il profilo, altrimenti verrà visualizzato il profilo coi parametri precedentemente calcolati.

5 - Plot
La visualizzazione dei profili si basa sul modulo Matplotlib, che fornisce gli usuali strumenti di zoom, pan, e salvataggio in formato grafico della visualizzazione (SVG, PDF, jpg, etc.).
Lo zoom non funziona in maniera sincronizzata fra i due grafici che rappresentano le elevazioni e le pendenze. Questo aspetto potrebbe essere risolto nelle prossime versioni.

6 - Save full results as
Tutti i dati vengono salvati come informazioni puntuali in formato csv oppure come uno shapefile 2D. Il valore di inclinazione attributo ad ogni punto di campionamento si riferisce al segmento che lo precede immediatamente. Da notare che se il profilo scelta ha il vertice iniziale al di fuori di un DEM utilizzato, la distanza cumulativa 3D non sarà calcolata per nessun punto.

7 - Export 3D line from DEM
E' possibile salvare le informazioni come dati lineari 3D per uno specifico raster. Il formato di output è sempre shapefile.


Possibili modifiche future
Alcuni spunti per possibili futuri miglioramenti del plugin come tipologie di analisi sono elencati nei link sottostanti.
Esempi di integrazione di Grass con R, per il calcolo della dimensione frattale lungo profili topografici, sono riportati in Raster profile along arbitrary line segments del California Soil Resource Lab. Esempi molto interessanti (ed avanzati) di integrazione tra profili topografici ed informazioni geologiche sono presentati in Grass Gis et R: superposition des couleurs d'un raster quelconque sur un profil topographique élaboré à partir d'un MNT di di Martin Laloux, come pure Python: utilisation des couches vectorielles et matricielles dans une perspective géologique, sans logiciel SIG, sempre di Martin Laloux.
E' prevista anche la possibilità di personalizzare la vestizione del profilo, con scelta del rapporto di scala e delle caratteristiche grafiche.


Bibliografia

Turcotte, D. L., 1997. Fractals and Chaos in Geology and Geophysics. Cambridge University Press, 398 pagine.