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)

geoSurfDEM: a C++ console application for determining intersections between 3D geological surfaces and topography

The determination of the theoretical intersections between digital 3D geological surfaces and topography could be of potential help for studying the field attitudes of natural geological surfaces, as mapped from outcrops or from aerial and satellite images. Since geological structures have complex geometries, the analysis of the relationships between 3D surfaces and topography requires tools that can process 3D geological surfaces.

geoSurfDEM aims at determining:

a) the theoretical intersections between a 3D surface and a topography
b) the local 3D attitude of that surface at each intersection point

How does it work?

Below you see the screenshot of an application run in a Linux shell. When compiled for Windows, the procedure is identical. The total run time can be quite long, many minutes or more.

What is to note?

After the application header display, the user is asked for the name of a text file. In this example, "input_files.txt" is provided, the name of a file located in the same directory as the running application. This file provides the paths of three files:

1) DEM, in ESRI ASCII grid format
2) geological surface, in (old) VTK text format
3) output csv file storing for each row: x-y-z-dip direction-dip angle

An example of input text file is the following:

./publ_data/dem_malpi_aster_wgs84utm33.asc
./publ_data/geological_plane.vtk
./publ_data/intersections.csv


Examples of input data files (i.e., DEM ASCII grid, VTK geosurface file, CSV intersection result) are present in the publ_data subdirectory.

Afterwards, the application outputs a few informative messages about the number of found features and at the last prints out the number of found intersecting points, hopefully greater than zero. The results are stored in the text file referenced by the third path in the input text file.

Example of use

To present the application and check the validity of its results, we use a theoretical test case, i.e. a geological plane with a desired attitude 135°/35°, and with a spatial extent fitting that of the test DEM, covering the Mt. Alpi zone (Basilicata, Southern Italy), derived from global ASTER data.
You can export a DEM in ESRI ASCII grid format with Saga GIS (in addition to ArcGIS).

Creation of test geological plane

The geological plane is created and saved as a VTK text file with simSurf. With this Python 2.7 tool, it is possible to simulate geological surfaces by using analytical formulas.

simSurf is subdivide in two modules:

a) geosurface_simulation.py: creates, geolocates and saves/exports an analytical surface
b) geosurface_deformation.py: reads an analytical surface created by the previous module, deforms it and saves/exports.

Horizontal plane creation

So we start creating a horizontal plane with the Geosurface simulation tool, Analytical formula part, see figure below.

The zero in the formula section is for the horizontal plane creation. You calculate the matrix and you can see the plane in three dimensions.

Then to the geographical parameters, that have to fit the DEM extent without creating an excessively large geological plane.

We create the simulated geosurface, optionally view it in three dimensions and then have to export it in the Geo Analytical Surface (GAS) format, i.e. a jason format.

Plane rotation

We then pass to the Geosurface deformation tool, import the previously exported jason file and then apply a rotation to the plane around a N-S horizontal axis, by 35°.

Apply and then rotate by 45° around a vertical axis (plunge equal to 90°).

In this way we obtain a plane dipping 35° towards N135°.

Plane displacement to DEM extent

Now we locate the rotated plane to a geographical position that broadly fits with the DEM. I choose to use my qgSurf plugin for QGIS for quickly locating a point at the center of the used DEM, while knowing also the z value.

You see to the right the coordinates (x-y-z) of the point at the DEM center, showed within QGIS.

I copied and pasted these values in the simSurf displacement tab, so that the plane is displaced by the given delta-x, delta-y and delta-z amounts.

Done, after applying.


Save the geosurface as a VTK file and then you can see it in Paraview and use in the geoSurfDEM application. Note that the VTK file stores the plane as triangle mesh, without explicit attitude (i.e., dip direction and angle) information. So the local results calculated by the geoSurfDEM application are derived by the local geosurface triangle attitude stored in the VTK file. Using a simple plane obviously we expect the same results for all intersection points.

Input data preview

We see how are the DEM and the VTK plane data in Paraview.

You can import the DEM when in x-y-z format (could create with Saga), then applying a
Table to Point filter, while the VTK format is directly read from Paraview. Here a nadiral view.
Y axis represents the North.

And a lateral one, as seen from the South.

geoSurfDEM result

At the end, what are the results of the geoSurfDEM application?

We see them displayed in Paraview, by importing the resulting csv file and superposing on the DEM points and the plane surface. The results are symbolized by blue dots. You see them following the visual intersection between the plane with dip direction 135° and dip angle 35° and the DEM.

 

Always in Paraview we see, for a few records, that the corresponding point attitudes calculated by geoSurfDEM are as expected: 135°/35° for each point, since in this test case we were dealing with a geological plane.

----------

The code repository of geoSurfDEM is at https://gitlab.com/mauroalberti/geoSurfDEM

Executables for Linux are also availables for Linux (64 bit) at the code repository.

EDITS: 

2023-01-01: fixed broken llinks to images

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.

"Great Sand Dunes National Park - the tallest dunes in North America. Photo © copyright by Jack Brauer." From: http://www.mountainphotography.com/photo/dunes-ripples/

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:

• sin( a * a + b * b )
• a * b + 1000
• cos( a ) * 200

Note that the chosen functions are used to create Numpy arrays, so word the functions according to the Numpy nomenclature in order to use existing Numpy functions.

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

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.

The currently implemented deformation types are:

1. displacement
2. rotation
3. scaling
4. simple shear (horizontal)
5. 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:
a₁₁ = cos ω + cos² α ( 1 - cos ω )
a₁₂ = - cos γ sin ω + cos α cos β ( 1 - cos ω )
a₁₃ = cos β sin ω + cos α cos γ ( 1 - cos ω )
a₂₁ = cos γ sin ω + cos α cos β ( 1 - cos ω )
a₂₂ = cos ω + cos² β ( 1 - cos ω )
a₂₃ = - cos α sin ω + cos β cos γ ( 1 - cos ω )
a₃₁ = - cos β sin ω + cos α cos γ ( 1 - cos ω )
a₃₂ = cos α sin ω + cos β cos γ ( 1 - cos ω )
a₃₃ = cos ω + cos² γ ( 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, S_x, S_y and S_z (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 "-γ sin² α": 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.