Sei giorni ancora per la sottoscrizione a favore dei volontari che traducono in italiano il manuale di Quantum GIS:
http://www.kapipal.com/032645076697408c90a929dc07354bac
Wednesday, 21 May 2014
Friday, 9 May 2014
Interactive map of small particle air pollution in world cities
In the Guardian web site an interactive map, based on the data of the World Health Organization, is available:
http://www.theguardian.com/news/datablog/interactive/2014/may/08/exposure-air-pollution-city-map
North America has lower pollution levels than in European states, with south-eastern European countries having the highest pollutions for Europe. India, China, Mongolia, Pakistan, Afghanistan, Iran present the highest levels in the world, together with many Middle East cities.
http://www.theguardian.com/news/datablog/interactive/2014/may/08/exposure-air-pollution-city-map
North America has lower pollution levels than in European states, with south-eastern European countries having the highest pollutions for Europe. India, China, Mongolia, Pakistan, Afghanistan, Iran present the highest levels in the world, together with many Middle East cities.
Monday, 5 May 2014
Traduzione del manuale di QGis in italiano: vuoi contribuire con una donazione?
"QGIS è un software open source molto potente sviluppato interamente da volontari. Ma cosa sarebbe un software senza un'accurata documentazione?
Una nuova versione di QGIS viene rilasciata ogni 4 mesi e di conseguenza la documentazione cambia ed è necessario aggiornare il manuale con le nuove funzionalità in modo che tutti siano in grado di sfruttarlo al massimo.
Siamo un gruppo di volontari che nel tempo libero crea, mantiene e traduce la documentazione di QGIS. Aiutaci con un piccolo contributo!"
Una nuova versione di QGIS viene rilasciata ogni 4 mesi e di conseguenza la documentazione cambia ed è necessario aggiornare il manuale con le nuove funzionalità in modo che tutti siano in grado di sfruttarlo al massimo.
Siamo un gruppo di volontari che nel tempo libero crea, mantiene e traduce la documentazione di QGIS. Aiutaci con un piccolo contributo!"
Thursday, 17 April 2014
Interesting compilation of GIS-related open source Python modules
A listing of many Python modules for GIS tasks, from spatial analysis to interpolation, projections, geodatabases and web mapping, is available at Python GIS Resources.
Monday, 14 April 2014
Python 2.7 End of Life: 2020 or later
Life for Python 2.7 is extended to at least 2020, giving more time to developers to migrate their code to Python 3.
You can read more at: http://hg.python.org/peps/rev/76d43e52d978
You can read more at: http://hg.python.org/peps/rev/76d43e52d978
Friday, 28 March 2014
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.
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.
Monday, 10 March 2014
Please support the upcoming Vienna Code Sprint 2014
In the occasion of the upcoming Vienna Code Sprint 2014 [1], the GRASS GIS Project Steering Committee decided will join the event.
While the GRASS developers are donating their valuable time, the community of enthusiast users may contribute with donations, for further details.
As usual, all of the work done in the community sprint will be directly .
For your convenience, here our easy-to-use Paypal button:
For our alternative bank transfer option, please contact Martin Landa (landa.martin@gmail.com)
Thanks for your support!
The GRASS Developers Team
About GRASS GIS
The Geographic Resources Analysis Support System (http://grass.osgeo.org/), commonly referred to as GRASS GIS, is an Open Source Geographic Information System providing powerful raster, vector and geospatial processing capabilities in a single integrated software suite. GRASS GIS includes tools for spatial modeling, visualization of raster and vector data, management and analysis of geospatial data, and the processing of satellite and aerial imagery. It also provides the capability to produce sophisticated presentation graphics and hardcopy maps. GRASS GIS has been translated into about twenty languages and supports a huge array of data formats. It can be used either as a stand-alone application or as backend for other software packages such as QGIS and R geostatistics. It is distributed freely under the terms of the GNU General Public License (GPL). GRASS GIS is a founding member of the Open Source Geospatial Foundation (OSGeo).
Sunday, 2 March 2014
qProf supports on-the-fly reprojection and multiples profile lines
The new version, 0.2.2, of this plugin for Quantum GIS has now support for on-the-fly reprojection of layers in a project. The profile will be created with the CRS set by the user in the Quantum GIS project.
Moreover, it also possible to use multiple lines for the profile creation, that will be merged into one for the profile creation. Since it could happen that the order of the single lines is not the correct one for the profile creation, it has been added the option of reading the correct order from an integer field defined in the attribute table of the layer (see figure below).
A bug related to the saving of results as shapefiles in Linux has also been removed.
The new version can be installed by using the plugin manager of Quantum GIS, or by downloading it from the Quantum GIS plugin repository at http://plugins.qgis.org/plugins/qProf/ or at https://bitbucket.org/mauroalberti/qprof/downloads.
Moreover, it also possible to use multiple lines for the profile creation, that will be merged into one for the profile creation. Since it could happen that the order of the single lines is not the correct one for the profile creation, it has been added the option of reading the correct order from an integer field defined in the attribute table of the layer (see figure below).
 |
| Example of layer with multiple lines (4), that are digitized with an incorrect order (see in the attribute table), but with the correct order set in the field 'id' of the attribute table (also used for the labelling of the lines in the map). |
 |
| Profile created with the correct line order, set in the attribute table as in the previous figure. |
A bug related to the saving of results as shapefiles in Linux has also been removed.
The new version can be installed by using the plugin manager of Quantum GIS, or by downloading it from the Quantum GIS plugin repository at http://plugins.qgis.org/plugins/qProf/ or at https://bitbucket.org/mauroalberti/qprof/downloads.
Monday, 6 January 2014
Creating and deforming analytical surfaces in Quantum GIS: experimental tools in qgSurf plugin
Imagine you want to create a georeferenced sinusoidal surface, using a trigonometric function, and then deform it via simple shear. This surface may represent a folded geological layer, lately sheared.
This operation is possible using the presented tool, implemented in qgSurf, a plugin for Quantum GIS. Its purpose is to allow these types of operations, precisely: a) the creation of analytical surfaces in a geographical space; and b) their deformation using a few, well-known deformation matrices.
Modelling of geological surfaces
In GIS continuous surfaces are generally stored as lattices/grids (rasters) or as Triangulated Irregular Networks (TINs). Both present relative advantages and disadvantages, but more importantly they share two disadvantages: first, they are memory-less, in the sense of not storing any information on their physical or mathematical generation history, and second, they have finite spatial resolution, a limiting factor when dealing with function-derived surfaces, that have a theoretically unlimited spatial resolution. Consider for instance multiple scale folds in structural geology, with superposed folds of different orders, or, in geomorphology, ripples superposed on dunes. The finite resolution of grids and TINs does not allow to represent the spatial structure from coarse to fine scales.
On the other hand, point lattices and TINs have the advantage of being expressed by points. Transformations can by applied on points by using linear equations or equivalent matrices. Deformations can be expressed as matrices, that can be incrementally added to represent a set of successive deformations [1].
This plugin allows to save the generative and deformation history parameters of a geosurface by choosing its internal, experimental "Gas" format, a simple Jason format with the parameters written as Python dictionaries. From them, a lattice of points can be generated at convenience, with the user defined grid parameters, and the deformation matrices can be applied to the points, producing a new 3D surface via triangulations.
Module structure
In addition to the previous 'Plane geoprocessing' module, the plugin presents two new modules.
With the first module, named 'Geosurface simulation', analytical, georeferenced surfaces (here called geosurfaces) can be created using the same approach as in Saga Gis (Grid - Calculus - Function, as of Saga vers. 2.0.3), visualised and exported as VTK, Grass or the internal, experimental "Gas" (geological analytical surface) format, a Jason format used in this plugin for storing the geosurface parameters, instead of the discretised points as in the VTK and Grass formats.
With the second module, 'Geosurface deformation', geosurfaces saved in the Gas format can be loaded and transformed via displacement, rotation, scaling and simple shear.
Simulation of geosurfaces
2.5 D surfaces can be created as analytical functions of a, b coordinates: z = f( a, b ).
Ranges for a and b values are defined, as well as the number of grid columns and rows, to be used for the generation of the 3D surface.
An analytical formula is provided, in order to generate the analytical surface.
Some examples could be:
The analytical surface is created after pushing the "Calculate matrix" button, and it can be viewed in 3D by using the "View as 3D surface" button (see example below).
After the analytical surface creation, this surface can be georeferenced by using the commands in the 'Geographic parameters' widget. The following parameters have to be defined:
The figure below illustrates these concepts.
A georeferenced analytical surface is created after pressing the 'Create simulated geosurface' button and again can be viewed as a 3d surface by using the 'View as 3D surface' button.
From the 'Output' widget it is possible to save the geosurface in the VTK, Grass or 'Gas' format.
VTK and Grass formats are widely used formats that stores the parameters of the geometrical elements constituting a surface. In our case, they are the triangular faces defining the complete surfaces, by means of the coordinates of each points.
In the Gas file format, on the other hand, no geometrical information is stored as points or faces. The saved informations describe the analytical and geographical parameters as defined in the 'Analytical formula' and 'Geographic parameters', plus the deformational parameters when present (described in the following paragraph). From this information, a new geosurface can be created at will.
The deformation of analytical surfaces
Surface can be changed via displacements, rotations or strains, each one with its own matrix or vector representation [1]. Apart from displacement, that is calculated as a vector that is added to the initial point position, all the other types are expressed as matrices that are multiplied to the initial point positions in order to obtain the final one. Obviously more than a deformation type can be applied to the same original analytical surfaces, for instance a vertical simple shear followed by a displacement and then a rotation.
Analytical functions produce 2.5 D surfaces. Through subsequent deformations, such as rotation or shearing, they can become true 3D surfaces, where more than one z value is defined for a single x-y value pair (image below).
The currently implemented deformation types are:
Displacement
A geosurface can be moved in the space, without rotation or distorsion, by given offsets in the x, y and/or z directions.
The displacement is calculated as the sum of initial point and the shift vectors.
Rotation
A geosurface can be rotated around a rotation axis, characterized by given trend and plunge values, by a rotation angle ω.
The rotated position is calculated by multiplying the rotation matrix with the initial position vector (eqs. 3.11a-c in [2]).
where:
a11 = cos ω + cos2 α ( 1 - cos ω )
a12 = - cos γ sin ω + cos α cos β ( 1 - cos ω )
a13 = cos β sin ω + cos α cos γ ( 1 - cos ω )
a21 = cos γ sin ω + cos α cos β ( 1 - cos ω )
a22 = cos ω + cos2 β ( 1 - cos ω )
a23 = - cos α sin ω + cos β cos γ ( 1 - cos ω )
a31 = - cos β sin ω + cos α cos γ ( 1 - cos ω )
a32 = cos α sin ω + cos β cos γ ( 1 - cos ω )
a33 = cos ω + cos2 γ ( 1 - cos ω )
and α is the angle between the rotation axis and the x axis, β is the angle between the rotation axis and the y axis, γ is the angle between the rotation axis and the z axis, and ω is the rotation angle.
The angles between the frame axes and the rotation axis are automatically derived from the rotation axis trend and plunge values.
Scaling
The size of the geosurface is scaled along the frame axes by three scale factors, Sx, Sy and Sz (X, Y and Z in the figure below).
The transformation matrix is:
Simple shear (horizontal)
We consider a horizontal simple shear (parallel to the x-y plane) with angle ψ (psi), along a direction that makes an angle α (alpha) with the x axis, as in the figure below.
The parameters are entered in this window:
Following the matrix derivation in Ramsay and Huber (1983), p. 290, the transformation is given by:
where γ is equal to tan(ψ).
Note the minus sign in the term "-γ sin2 α": in Ramsay and Huber 1983, eq. C.14 the sign is given as positive, but it appears to be inconsistent with both the derivation and the practical application of the formula.
Simple shear (vertical)
A geosurface can be sheared in the vertical plane, by an angle ψ (psi), along a direction making an angle α (alpha) with the x axis.
The parameters are entered in this window:
The transformation is given by:
where γ is equal to tan(ψ).
Current limitations
The visualization of deformed geosurfaces may fails when using the plugin internal visualizer, based on Matplotlib. Apparently the problem is linked to some features of Matplotlib management of points, but at the moment the problem is not solved. Also in Grass the visualization of deformed 3D surfaces may be unsuccessful, even if the imported data is correctly visualised in 2D.
On the contrary, Paraview import and visualization appears to be correct and without problems, at least on the tested datasets.
Download and installation
This version is still experimental. It can be downoladed and installed autmatically from the QGis plugin manager, or downloaded from http://plugins.qgis.org/plugins/qgSurf/version/0.3.0/ and unzipped in the user Python plugin folder (e.g., C:\Users\mauro\.qgis2\python\plugins for Windows Vista).
References
[1] Ramsay, J. G., Huber, M. I., 1983. The techniques of modern structural geology. Volume 1: Strain Analysis. Academic Press, Inc. 307 pp.
[2] Allmendinger, R.W., Cardozo, N., Fisher, D. M., 2012. Structural geology algorithms. Vectors and tensors. Cambridge University Press. 289 pp.
Modelling of geological surfaces
In GIS continuous surfaces are generally stored as lattices/grids (rasters) or as Triangulated Irregular Networks (TINs). Both present relative advantages and disadvantages, but more importantly they share two disadvantages: first, they are memory-less, in the sense of not storing any information on their physical or mathematical generation history, and second, they have finite spatial resolution, a limiting factor when dealing with function-derived surfaces, that have a theoretically unlimited spatial resolution. Consider for instance multiple scale folds in structural geology, with superposed folds of different orders, or, in geomorphology, ripples superposed on dunes. The finite resolution of grids and TINs does not allow to represent the spatial structure from coarse to fine scales.
 | |
| "Great Sand Dunes National Park - the tallest dunes in North America. Photo © copyright by Jack Brauer." From: http://www.mountainphotography.com/photo/dunes-ripples/ |
This plugin allows to save the generative and deformation history parameters of a geosurface by choosing its internal, experimental "Gas" format, a simple Jason format with the parameters written as Python dictionaries. From them, a lattice of points can be generated at convenience, with the user defined grid parameters, and the deformation matrices can be applied to the points, producing a new 3D surface via triangulations.
In addition to the previous 'Plane geoprocessing' module, the plugin presents two new modules.
With the first module, named 'Geosurface simulation', analytical, georeferenced surfaces (here called geosurfaces) can be created using the same approach as in Saga Gis (Grid - Calculus - Function, as of Saga vers. 2.0.3), visualised and exported as VTK, Grass or the internal, experimental "Gas" (geological analytical surface) format, a Jason format used in this plugin for storing the geosurface parameters, instead of the discretised points as in the VTK and Grass formats.
With the second module, 'Geosurface deformation', geosurfaces saved in the Gas format can be loaded and transformed via displacement, rotation, scaling and simple shear.
Simulation of geosurfaces
2.5 D surfaces can be created as analytical functions of a, b coordinates: z = f( a, b ).
Ranges for a and b values are defined, as well as the number of grid columns and rows, to be used for the generation of the 3D surface.
An analytical formula is provided, in order to generate the analytical surface.
Some examples could be:
- sin( a * a + b * b )
- a * b + 1000
- cos( a ) * 200
 |
| Analytical surface with formula: sin(a*b) and ranges - 5 to 5 for both a and b, and grid columns and rows equal to 70. |
After the analytical surface creation, this surface can be georeferenced by using the commands in the 'Geographic parameters' widget. The following parameters have to be defined:
- the length and width of the georeferenced surface to be created,
- its rotation angle with respect to the x axis,
- the x and y values of the lower-left corner surface ('x min' and 'y min').
The figure below illustrates these concepts.
A georeferenced analytical surface is created after pressing the 'Create simulated geosurface' button and again can be viewed as a 3d surface by using the 'View as 3D surface' button.
 |
| Sinusoidal function created with the qgSurf plugin. The visualization is based on Matplotlib. |
From the 'Output' widget it is possible to save the geosurface in the VTK, Grass or 'Gas' format.
VTK and Grass formats are widely used formats that stores the parameters of the geometrical elements constituting a surface. In our case, they are the triangular faces defining the complete surfaces, by means of the coordinates of each points.
In the Gas file format, on the other hand, no geometrical information is stored as points or faces. The saved informations describe the analytical and geographical parameters as defined in the 'Analytical formula' and 'Geographic parameters', plus the deformational parameters when present (described in the following paragraph). From this information, a new geosurface can be created at will.
The deformation of analytical surfaces
Surface can be changed via displacements, rotations or strains, each one with its own matrix or vector representation [1]. Apart from displacement, that is calculated as a vector that is added to the initial point position, all the other types are expressed as matrices that are multiplied to the initial point positions in order to obtain the final one. Obviously more than a deformation type can be applied to the same original analytical surfaces, for instance a vertical simple shear followed by a displacement and then a rotation.
Analytical functions produce 2.5 D surfaces. Through subsequent deformations, such as rotation or shearing, they can become true 3D surfaces, where more than one z value is defined for a single x-y value pair (image below).
 |
| Sheared sinusoidal surface. Visualised with Paraview. |
- displacement
- rotation
- scaling
- simple shear (horizontal)
- simple shear (vertical)
Displacement
A geosurface can be moved in the space, without rotation or distorsion, by given offsets in the x, y and/or z directions.
The displacement is calculated as the sum of initial point and the shift vectors.
Rotation
A geosurface can be rotated around a rotation axis, characterized by given trend and plunge values, by a rotation angle ω.
The rotated position is calculated by multiplying the rotation matrix with the initial position vector (eqs. 3.11a-c in [2]).
where:
a11 = cos ω + cos2 α ( 1 - cos ω )
a12 = - cos γ sin ω + cos α cos β ( 1 - cos ω )
a13 = cos β sin ω + cos α cos γ ( 1 - cos ω )
a21 = cos γ sin ω + cos α cos β ( 1 - cos ω )
a22 = cos ω + cos2 β ( 1 - cos ω )
a23 = - cos α sin ω + cos β cos γ ( 1 - cos ω )
a31 = - cos β sin ω + cos α cos γ ( 1 - cos ω )
a32 = cos α sin ω + cos β cos γ ( 1 - cos ω )
a33 = cos ω + cos2 γ ( 1 - cos ω )
and α is the angle between the rotation axis and the x axis, β is the angle between the rotation axis and the y axis, γ is the angle between the rotation axis and the z axis, and ω is the rotation angle.
The angles between the frame axes and the rotation axis are automatically derived from the rotation axis trend and plunge values.
Scaling
The size of the geosurface is scaled along the frame axes by three scale factors, Sx, Sy and Sz (X, Y and Z in the figure below).
The transformation matrix is:
Simple shear (horizontal)
We consider a horizontal simple shear (parallel to the x-y plane) with angle ψ (psi), along a direction that makes an angle α (alpha) with the x axis, as in the figure below.
The parameters are entered in this window:
Following the matrix derivation in Ramsay and Huber (1983), p. 290, the transformation is given by:
where γ is equal to tan(ψ).
Note the minus sign in the term "-γ sin2 α": in Ramsay and Huber 1983, eq. C.14 the sign is given as positive, but it appears to be inconsistent with both the derivation and the practical application of the formula.
Simple shear (vertical)
A geosurface can be sheared in the vertical plane, by an angle ψ (psi), along a direction making an angle α (alpha) with the x axis.
The parameters are entered in this window:
Current limitations
The visualization of deformed geosurfaces may fails when using the plugin internal visualizer, based on Matplotlib. Apparently the problem is linked to some features of Matplotlib management of points, but at the moment the problem is not solved. Also in Grass the visualization of deformed 3D surfaces may be unsuccessful, even if the imported data is correctly visualised in 2D.
On the contrary, Paraview import and visualization appears to be correct and without problems, at least on the tested datasets.
Download and installation
This version is still experimental. It can be downoladed and installed autmatically from the QGis plugin manager, or downloaded from http://plugins.qgis.org/plugins/qgSurf/version/0.3.0/ and unzipped in the user Python plugin folder (e.g., C:\Users\mauro\.qgis2\python\plugins for Windows Vista).
References
[1] Ramsay, J. G., Huber, M. I., 1983. The techniques of modern structural geology. Volume 1: Strain Analysis. Academic Press, Inc. 307 pp.
[2] Allmendinger, R.W., Cardozo, N., Fisher, D. M., 2012. Structural geology algorithms. Vectors and tensors. Cambridge University Press. 289 pp.
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.
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).
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.
'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.
.
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.
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.
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.
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.
 |
| 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. |
 |
| 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)
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).
Subscribe to:
Posts (Atom)