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.