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