[Written by Claude, reviewed by myself]
You are standing on an outcrop with a compass in your hand, and the bedding reads 135/35 — a plane dipping thirty-five degrees towards the south-east. Written in a notebook, that is a statement about a square metre of rock. Drawn on a map it becomes a much larger claim: that the same surface will be found there, a kilometre away, on the far side of the valley, under the scree where nobody can go and check.
The distance between those two statements is the whole subject of this post. I put it on this blog on 30 December 2011 as a problem rather than a program, and the answer had two halves which have not changed since. Where the trace of a surface is not mapped, the computed intersection is a hypothesis about where to go and look for it. Where it is mapped, the computation is a test: vary the attitude, and see whether some orientation brings the calculated trace into agreement with the drawn one. Three weeks later there was a program, gSurf, and three weeks after that a QGIS plug-in, qgSurf.
gSurf slept for most of a decade and was rebuilt this year, thanks to the availability of LLM agents, i.e., Claude. It now has two tools, and they answer the same question from opposite ends: one takes a single measurement and extends it across the topography, the other takes a few hundred measurements and asks whether they belong to one folded surface. This post is about what each is for in the field, with the southern Apennines as the test ground throughout. I write attitudes as dip direction/dip, so 135/35 is the plane above, and fold axes as trend/plunge.
One window, one question
gSurf is a small desktop application: one map, one panel of controls, no project file, nothing to configure before you can ask something. It is not a GIS and is not trying to become one. Its single organising idea is that the answer is recomputed on every frame rather than behind a button marked Calculate, so a parameter becomes something you sweep through instead of something you guess, wait for, and guess again.
That sounds like a detail of interface design. It is actually a change in what the tool can tell you, and I will come back to it twice — once for each tool — because in both cases the useful output is not the number on the screen but how much that number can move before the picture stops fitting.
Three projects, and who does what
Behind the application there are two libraries, and the division of labour between them is worth explaining because it is decided by geology rather than by programming taste.
misah is a kernel: about twenty functions written in Rust, compiled, and called from Python. It does the arithmetic that has to run over millions of grid cells — the intersection of a plane with a DEM, among others. I wrote about it last week. geogst is a Python library with no user interface: geometries, orientations, the statistics of a set of measured directions, stereonets, profiles, coordinate systems. gSurf is the hand on the control, and almost nothing else — it opens files, draws a map, and keeps the frame rate up.
Which of the two libraries a tool leans on tells you what kind of question it is asking. Laying a plane on a DEM means visiting every cell of a grid: a million of them per frame, and in pure Python that is 1.7 seconds — fine for one answer, hopeless for a dial. Compiled, the same work takes 22 milliseconds, a factor of about seventy-seven, and that factor is the entire reason the dial exists. Finding a fold axis, by contrast, means taking the statistics of a few dozen directions: 0.4 milliseconds for a window of twenty measurements, 4.2 for three hundred. There is no problem there to solve, so that tool calls geogst directly and no compiled code is involved at all. Writing a faster version of mathematics that already costs less than a millisecond would have bought nothing, and cost a second copy of a formula to keep correct.
Why a second application, next to a plug-in
qgSurf is the one with users. It is on the QGIS plugin repository, it has grown well past the original question — best-fit planes, distances to a plane, stereonets, and since 2022 the cross-section tools of qProf — and it runs inside the GIS where your data are already loaded. If you want to do this on your own project, that is the one to install.
gSurf is where an idea gets tried before it is worth putting in front of those people. The difference is not ambition, it is what each is allowed to risk. gSurf can depend on an alpha kernel whose function signatures still move; it can be rewritten in a weekend; it can throw its whole interface away, because nobody's Tuesday depends on it. A plug-in with users cannot do any of those things, and should not.
What comes back the other way is the reason the arrangement earns its keep. Two examples from this summer, both discovered in gSurf because there the cost of a frame is on the screen while you work. Drawing thousands of short line segments as one broken polyline instead of thousands of separate objects is four to five times cheaper — that lesson went straight into the plug-in, which had been creating one marker per computed point and then walking over all of them at every pan and zoom. And a convenient little helper that reprojects a coordinate pair turned out to be rebuilding the whole projection machinery on each call: on twenty thousand points, 215 milliseconds down to 32. Neither of those is a new algorithm. They are the sort of thing you find when the same calculation has to finish sixty times a second, and not otherwise.
The dial: a plane on a DEM
The first tool lays an unbounded plane — a bedding surface, a fault, a contact, extended forever in its own plane — through a point you choose on the map, and draws the line where it meets the topography. Turn the dial and the trace moves.
Fig. 1. Timpa San Lorenzo, Calabria, on a 5 m DEM with the geological map and the mapped faults underneath. The plane is 135/35 through the yellow point at 758 m; the red line is where it crops out. The dashed rectangle is the window the calculation actually reads. The status bar gives the cost of this frame — 853 points, 5.5 ms of arithmetic, 3.1 ms of drawing — and the grey line under the attitude reads convergence +0.83°, grid 134.2°, which is the subject of a paragraph further down.
There are four things a field geologist can do with that, and they get progressively less obvious.
Predict where to look. You have one good exposure of a contact and a day to find more of it. The trace tells you which spurs it should cross and at what elevation, which is a walking plan.
Test whether a measurement is representative. This is the example I used to introduce the plug-in in February 2012, on a normal fault in the Valnerina, Umbria, from my own mapping. The measured attitude was 227/58, and the trace it produced did not follow the fault as drawn on the map. Turning it to 219/61 did. Neither number is wrong: the first describes the outcrop, the second describes the surface. Eight degrees of dip direction and three of dip separated them — which is about what two people reading the same fault plane will disagree by, and exactly the amount that moves a contact off one ridge and onto the next. I then checked the fitted value the way it was done before any of this software existed, by measuring the separation of the computed trace between contours: 450 m of relief over 245 m of horizontal distance gives a dip of 61.4° against the 61 the program had been given, and a dip direction of about 220 against 219.
See how well constrained the answer is. Here is what the button used to hide. In 2012 I found 219/61 by typing numbers and pressing Calculate; what I could not see was whether 217 or 223 would have done just as well. Sweeping the dial answers that in one gesture: the range over which the trace still follows the mapped contact is the uncertainty of the estimate, read straight off the map. On strongly varied topography that range is a couple of degrees. On a planar hillside it is enormous — every plane through that slope produces much the same line — and watching it happen is a better lesson in ill-conditioned problems than any amount of being told.
Project a horizon that is not at the surface. The three coordinates of the source point are independent, so you can lift it off the ground and keep it there while you drag it around the map. That is how you lay a plane on a marker horizon that passes a hundred metres above today's topography, or below it.
Two traps, both geological rather than computational. The first is that the two norths do not agree. You measure dip direction against true north, once declination is corrected; the DEM lives on a projection's grid, and grid north points somewhere else. In the southern Apennines on EPSG:25833 the meridian convergence runs between +0.41° and +1.04°, which over five kilometres of trace is up to 91 metres — twenty DEM cells, small enough to look like a mapping error and large enough to make you doubt a perfectly good measurement. gSurf reads and writes true azimuth, as the compass gives it, subtracts the convergence before computing, and shows you both numbers. The second is that the DEM is never loaded into memory: the background is a coarse overview and the calculation reads one full-resolution window at a time. A 234-megapixel mosaic therefore costs the same per frame as a small crop, where loading it in the obvious way would have been 2.5 GB.
The cost of a frame, on that mosaic, dragging the dial:
| window | ms per frame | frames per second |
|---|---|---|
| 500² cells (2.5 km) | 10.9 | 92 |
| 1000² cells (5 km) | 26.6 | 38 |
| 2000² cells (10 km) | 95.7 | 11 |
Real time holds to five kilometres of window and no further. Beyond that the arithmetic alone eats the frame, and the honest thing is to say so rather than to let the dial feel sticky.
The circle: fold axes under a moving window
The second tool starts from the other end: not one measurement, but every bedding attitude on a map sheet. Drag a circle across the map, and it reports what the measurements inside that circle have in common.
The reasoning is the classical π-diagram, and it is worth stating carefully because everything the tool does follows from it. Plot the pole of each bedding plane — the direction perpendicular to it — on a stereonet. If the beds inside the circle are folded about a single line, then every bed contains that line, so every pole is perpendicular to it, and the poles fall on a great circle (a girdle). The pole of that girdle is the fold axis. If instead the beds all dip much the same way — a homocline — the poles fall in a tight cluster, and there is no fold axis to report at all.
Fig. 2. The circle at 590000/4500000 on the Potenza–Irsina sheet, radius 2 km, with 38 bedding attitudes inside it (red) out of those on the map (black). The net follows the circle as you drag it — poles as dots, the best-fit girdle dashed, and its pole, the fold axis, as the red diamond: 312/13, with K = 0.59. It normally floats over the map and is docked here to fit in one picture. The status bar times the frame you are looking at: 0.46 ms to find the attitudes inside the circle, 1.3 ms for the statistics, 5.4 ms to draw both pictures.
Telling those two cases apart is not a matter of judgement, and it is the one piece of statistics in this post. Woodcock's K compares how elongate the set of poles is with how flattened it is: below 1 the poles form a girdle and there is a fold; above 1 they cluster and there is not. C says how strongly oriented they are at all, as against scattered. gSurf refuses to call anything a fold axis unless K is at most 1, C is at least 1, and there are at least ten measurements — and you can move all three thresholds while you work.
The refusal matters more than it sounds. In a cluster, the direction the tool would report as an axis is the direction with the fewest poles, which is the worst-determined direction in the data rather than a structure. It will happily produce a number. On the 1757 usable attitudes of the Potenza–Irsina sheet, three windows in four fail that test, and taking the sheet as a whole gives K = 1.28 — that is, over most of that map the bedding is homoclinal, and an axis computed there means nothing. A refused axis is drawn in grey rather than hidden, because hiding it would answer "is this a fold?" with a blank screen, which reads the same as no data; greyed, you watch it turn colour as the circle crosses a real hinge.
The radius is the interesting control. Read the same place at three sizes and you see what the window is for.
Fig. 3. One locality, three radii, the same ground each time. At 1 km there are ten attitudes, they cluster (K = 1.16), and the axis is refused — drawn grey on the net rather than hidden. At 2 km, 38 attitudes give 312/13 with K = 0.59. At 5 km, 162 attitudes give 312/19. An axis that survives a change of scale is a structure; one that does not is a coincidence.
The radius is also a decision about the data, not only about the geology. On Monte Alpi the survey is thin — 0.69 measurements per square kilometre against 1.10 on the Val d'Agri sheet to the north — and the median window at 2 km holds seven attitudes where at 3 km it holds eighteen. Below ten the tool refuses, so on that massif the choice is 3 km or nothing, and you should know that you are making it. In the Pollino massif the question cannot be asked at all: sheet 534 Castrovillari was never published, and Monte Pollino has no bedding measurement within ten kilometres of it. The tool has nothing to say there, and says so.
Asking everywhere at once
Dragging a circle by hand answers one question at a time. The same window can be laid down at every node of a square grid instead, which turns a sheet of scattered measurements into a map of fold axes — drawn only where the test is passed, and coloured by plunge, because a tick mark can carry the trend in the direction it points but not the plunge in its length, and the plunge is half of what an axis is.
Fig. 4. Left: 1835 fold axes over two adjacent 1:25,000 sheets — one tick for every 500 m cell whose window passed the test, out of 6270 cells that held data, coloured by plunge. The grey dots are the bedding stations; the blank ground between the patches is where the poles clustered instead of forming a girdle. Right: those same axes as trend against northing (grey), with their mean per kilometre of northing in red. The mean turns by fourteen degrees across the dashed line — and the dashed line is not a structure. It is where one map sheet ends and the next begins.
This is where a tool has to be honest about what it is doing, because a picture of a thousand ticks is extremely persuasive. Three things are worth knowing before believing one.
A step in the data is not a step in the rocks. Averaged over the six kilometres on either side, the axes of the two sheets above trend 157° and 138° — nineteen degrees apart, which one would happily interpret as a structural boundary. Before interpreting it as anything, it is worth asking what the same measurement does across a line that means nothing: cut the field at parallels every 250 m inside a single sheet, and compare three kilometres above each cut with three below. Eighty-four such lines give a median step of 3° and a largest of 10°. The sheet edge gives 21°, and the six-kilometre and one-kilometre versions of the test agree. So the step is real, and it falls on an administrative line — but what put it there is still open, and the obvious answer is the wrong one. A change of recording convention can be ruled out: both sheets give azimuth to the degree and dip to five, in much the same proportions, and rewriting a sheet's attitudes to the coarser convention shifts the computed trend by 0.7° in the median. Whatever the two campaigns did differently, they did not do it with their notebooks. What is left is which outcrops were visited and which surfaces were called bedding — and neither is recoverable from the delivered map, which is the uncomfortable part. What makes the question askable at all is merging the two sheets into one layer: asked inside either one, no test can see it.
Overlapping windows are not independent observations. A circle of radius R laid down every S metres covers πR²/S² cells, which at R = 2 km and S = 500 m is about fifty. Every attitude is therefore counted in fifty windows, and neighbouring ticks share nine tenths of their data. A field of a thousand axes looks like a thousand observations and carries about fifty windows' worth of information. The tool prints both numbers under the spacing box — each attitude in ~50 cells; ~50 windows would tile the area — because the ratio between the two controls is never something you set on purpose, and the count of axes is what gets quoted. A finer spacing buys resolution in the picture and no further information underneath it.
The threshold is a choice, and it shows. Because the field is kept as numbers rather than as a picture, moving the K threshold re-decides the entire map in three to seven milliseconds without recomputing anything. Sweeping it from 0.5 to 2.5 takes the same sheet from 369 axes to 2049, out of 3556 cells that hold data. How much of a map depends on where you put a threshold is not a question you can answer by running the calculation four times over a coffee break — but it is one you can answer by dragging a slider, and then it is not a question you can avoid. The exported grid carries the thresholds along with the verdict for every cell, because an answer recorded without the choice that produced it cannot be checked afterwards.
Building the field is the slowest thing either tool does, and even that is not slow: 4140 cells in 2.6 seconds, 16289 in 9.6. Dragging the circle by hand costs 4.5 milliseconds a frame, of which the search for the attitudes inside it is 0.1 and the statistics 0.8; the rest is drawing two pictures.
What the field of axes does not tell you
One result from this summer's work in the southern Apennines belongs here, because it is the sort of thing a tool should be used to find out about itself.
On two sheets in northern Calabria there are 94 fold axes measured directly at outcrop — someone stood at the hinge and read the line. Computed axes are available at the same places from the bedding attitudes on the same sheets, so the two can be compared, and I expected the measured ones to serve as a yardstick. They do not agree: the median difference in trend is 34° to 44°. What does agree is the plunge — 15.0° against 15.2° — and the computed field is internally consistent, in that neighbouring nodes differ from each other by about four degrees. So this is not noise, and it is not a coding error either.
The most likely reading is that the two are not measuring the same thing. A 2 km window pools bedding across whatever is inside it and returns the geometry that best describes the pool, which in an area folded more than once is not necessarily the generation whose hinge you were standing on. The conclusion I draw for the tool is narrow and useful: in that area the absolute computed axes should not be trusted, while the way the field changes from place to place still can be. A window is a scale of observation, and an answer returned at one scale is not a claim about all of them.
What the two tools have in common
They look different — a dial and a circle, one measurement and four hundred — and they are the same idea twice. Both take something whose validity depends on a scale, and put that scale in your hand. Both are built so that the control can be moved continuously, because the answer that matters is the behaviour of the number and not the number. And both are willing to say nothing: a trace that fits over a forty-degree range of attitudes, an axis drawn in grey because the poles do not form a girdle. A tool that always produces an answer is not being helpful.
Getting it
gSurf is at gitlab.com/mauroalberti/gSurf, GPL-3, Python 3.9 or later:
pip install misah numpy rasterio PyQt6 matplotlib pyproj pip install geogst mplstereonet geopandas shapely # for the fold axes python -m gsurf
It opens on the two tools and asks for nothing until you have picked one; then it asks only for what that tool needs — the plane wants a DEM, the fold axes want a layer of attitudes and the two columns its angles are in, and neither is asked for the other's.
It is alpha and the word is meant. Two things from the older version are worth having back and are not in it yet: topographic profiles with attitudes projected onto them, and the fault-and-slickenline stereonet that reads a rake and a sense of movement. If you want these analyses on data you already have loaded, the plug-in — qgSurf — is the place to start, and it is where whatever survives its trial here will eventually turn up.
The trail on this blog
- 2011-12-30 — the problem, before any code
- 2012-01-18 — gSurf 0.1.0
- 2012-02-10 — the plug-in, and the Valnerina test
- 2024-08-26 — geogst
- 2026-09-09 — misah, the compiled kernel
Modified: 2026-09-18 H 14:32
