A dial and a circle: two ways of asking a map a question

[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. Then look at where the change happens: 151° in the last kilometre of the southern sheet and 137° in the first kilometre of the northern one, with both sides steady to within a few degrees beyond that. A 2 km window smears even a sharp geological hinge over a couple of kilometres, so a transition that abrupt cannot be a structure at all; it is a change of surveyor, of convention, or of mapping campaign. And what makes the question askable 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

Misah

Given an attitude measured at an outcrop, where must that plane cross the topography? That calculation has been rewritten several times. It is now a compiled kernel called misah, written mostly by Claude (as is this post), and as of this week it is installable:

pip install misah

What it is, and what it is not

misah is not an application. There is no window, no menu and no project file. It is a library of about twenty functions that take arrays of numbers and return arrays of numbers — written in Rust, compiled, and called from Python as if it were an ordinary module. You give it a DEM and a plane and it gives you back the trace. You give it faults and slickenlines and it gives you back a stress tensor.

The reason it is compiled is speed, and the reason speed matters is not impatience. A calculation that takes two seconds is something you run when you have decided what to ask. A calculation that takes five milliseconds is something you can put on a dial and turn, and turning a dial is a different way of thinking about a problem — you can sweep dip direction through forty degrees and watch the trace walk across the valley, rather than guessing a value, waiting, and guessing again. The last section of this post is about what that looks like.

Nothing here is a new algorithm. Most of it is a port of code I have written before and checked against for years: the plane–DEM intersection and the best-fit plane come from geoSurfDEM, the forward stress problem from a Fortran program of mine called ForwardStress.f95, the 3D density from InterpDensity3D, and the mechanism rotations from Kagan (1991), against the table printed in that paper. The tests do not check that the new code runs. They check that it agrees with the old code, number for number.

Where does that plane crop out?

This is the original question, and it is still the one I use misah for most. An unbounded plane, an outcrop point it passes through, a DEM, and the answer is the set of chords where the two surfaces meet.

 

Fig. 1. One outcrop at 1527 m on Monte Alpi, southern Apennines, and the same plane at three attitudes on an ASTER DEM at 27 m. Red is 135/35; blue is the same dip direction ten degrees steeper; green is fifteen degrees round in dip direction. All three pass through the same point, and by the northern edge of the map they are more than a kilometre apart. This is why an attitude measured at one outcrop is a hypothesis about the map and not a description of it.

Look at how far the three traces separate. That is the whole argument for doing this rather than sketching it: ten degrees of dip is well inside what two people reading the same fault plane will disagree about, and by the far side of the map it has moved the contact off one ridge and onto another. In the Valnerina case the difference between 227/58 and 219/61 was eight degrees and three degrees — and it was the difference between a trace that followed the mapped fault and one that did not.

Stress, forwards and backwards

The Wallace–Bott hypothesis says slip on a fault happens in the direction of the resolved shear stress on it. Run forwards, that predicts the slickenline a given stress tensor would leave on a given plane. Run backwards, it recovers the tensor from a population of faults with their slip.

 

Fig. 2. Left: eight fault planes of assorted strike and dip, and the slickenlines that an Andersonian normal tensor (S1 vertical, Φ = 0.5) would drive on each. Right: those faults handed back to the inversion with the tensor removed. Circles are the true axes, stars the recovered ones; they coincide. The search visited 71,280 candidate tensors and the mean misfit at the winner is 0.00°.

The search is exhaustive rather than a descent, and that is a deliberate choice with a geological reason. A fault population carrying two superposed tectonic phases has two minima by construction. A descent finds whichever one it falls into and reports it with no hint that the other exists. Walking the whole grid costs more and tells you when your data are of two minds — the function returns the runners-up as well as the winner, so a minimum standing alone can be told from one sitting on a plateau.

Two things it insists on, both of which have cost me data in the past. A slickenline that does not lie in the plane it was recorded on is refused by name rather than inverted: that is not a fault with a small error, it is two compass readings that do not belong to each other. And where the sense of movement was never determined, you say so — a fault whose sense nobody could read is scored modulo 180°, instead of being scored at 180° for being right.

Density without bins

Counting hypocentres in boxes has two problems that everyone knows about and nobody can fix: you choose the box size by hand, and the answer changes when you move the boxes. A kernel estimate removes both.

 

Fig. 3. 420 synthetic epicentres and their density at 900 m bandwidth. The field is in events per km², not in arbitrary units: summed over the grid and multiplied by the cell area it comes back to 420.0 of the 420 events that went in. That is what makes two runs at different bandwidths comparable with each other.

The bandwidth is in map units and there is one per axis, which matters more than it sounds. Depth is not interchangeable with easting even when both are in metres: a seismogenic layer is far wider than it is thick, and a single isotropic bandwidth over one will mix its top with its bottom before it mixes two neighbours. The same machinery runs in one dimension too, where a histogram of hypocentral depths becomes a smooth curve with no bin edges to argue about.

Over the same grid you can put a stress inversion at every node instead — Hardebeck and Michael's (2006) spatially varying inversion, in kernel form: each fault contributes to every node it can reach, weighted by distance, rather than the dataset being cut into bins and neighbouring solutions damped towards each other. Every node also reports how much data actually reached it, which is the number that says whether the tensor there is worth reading.

Turning the dial: gSurf

gSurf is the 2012 application, rebuilt this year around the misah kernel. It does one thing: it lays an unbounded plane on a DEM and recomputes the intersection on every frame, while you turn a dial.

Fig. 5. gSurf on Monte Alpi region, with the geological map and the mapped faults underneath. Dial for dip direction, slider for dip angle, and the red trace follows as you turn them. The status bar at the bottom reports the cost of every frame: 853 points, 5.5 ms in the kernel, 3.1 ms drawing. Note the line under the attitude — convergence +0.83°, grid 134.2°.

That last detail is worth a paragraph of its own, because it is the kind of thing that quietly ruins a result. You measure dip direction against true north. The DEM lives on a projection's grid, and grid north is not true north. In the southern Apennines the meridian convergence runs between +0.41° and +1.04°, which over five kilometres of trace is up to 91 metres of displacement — 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 it is measured in the field, and subtracts the convergence before calling the kernel.

The DEM is never loaded into memory, so it can be as large as you like: the background is a decimated overview and the kernel runs on a full-resolution window centred on the source point. You can drop the point anywhere on the map, lay the plane on a horizon that passes above today's topography, and export the trace as a shapefile when it fits.

Getting it

misah is at version 0.2.0-alpha.2, on PyPI and on crates.io, GPL-3.0-or-later, and needs Python 3.9 or later:

pip install misah

There are wheels for Linux (x86-64 and ARM64), macOS (universal2) and Windows, so there is nothing to compile unless you want to. The source is at gitlab.com/mauroalberti/misah, with three worked notebooks under docs/notebooks that generate their own data through the forward model — so what the inversion ought to return is known without trusting the inversion. gSurf is at gitlab.com/mauroalberti/gSurf.

It is an alpha and the version number means it: the function signatures can still move. What will not move is the arithmetic, which is older than any of the code it currently lives in.

 

 

Building geological cross-sections with GeoProfiler: simpler and more flexible

[Thanks to chatGPT for helping]

The GeoProfiler module in qgSurf is an attempt to make geological profile creation workflow less painful and more reproducible. GeoProfiler is built around a three-step workflow:

1.    (optional) Create (or reuse) a 3D topographic profile layer 
2.    Tell GeoProfiler which 3D profile layer you want to work on
3.    Add geological data and generate a configurable plot

Steps 1 and 3 are the “real work”. Step 2 is the glue: a small but crucial step where you just define the working 3D line layer.

What makes GeoProfiler really useful, especially from a field-geologist’s point of view, is that it allows you to think in terms of geological problems, not in terms of software tools.

Very often, when we build a cross-section, the hardest part is not drawing the line — it is deciding how the terrain should be cut in order to see something meaningful. Sometimes we want a single clean section, sometimes we want to see how the same structures evolve a few hundred meters to the left and to the right.

With GeoProfiler, this becomes natural. You can create a single topographic section, when you just need one good line across your structure — or a small set of parallel sections, when you want to compare how things change laterally. This is especially helpful for understanding the geometry of thrust systems, fold trains, or fault zones that widen or branch with distance. Instead of repeating the work five times, you simply define one baseline and let the tool generate parallel profiles around it.

Another big step forward is that profiles are no longer forced to be straight lines.
In earlier tools, the section had to be rigid: point A to point B, perfectly linear. But that is rarely how we actually think in the field. Sometimes the most meaningful section follows a valley floor, bends around a ridge, or gently curves to intersect key geological features. GeoProfiler allows those curved (but still geologically reasonable) sections, so the resulting profiles feel closer to the way we would sketch them in a notebook.

A third aspect that matters a lot in real projects is that our datasets almost never “fit” each other perfectly.
A DEM from one source, vector geology from another source, structural stations collected with yet another coordinate system. GeoProfiler is designed to work across different map projections, handling the transformations internally so that everything is brought consistently onto the profile. Instead of fighting with reprojections and fear of small misalignments, you can focus on what each dataset is telling you.

Once the profile (or profiles) exist, the geological part becomes straightforward: you can start adding what you would normally reason about in a cross-section — fault traces intersecting the terrain, lithological units cutting across the line, scattered point data such as measurements or events, or structural attitudes that define planes dipping into or out of the section. Bit by bit, the profile becomes not just a topographic cut, but a genuine geological interpretation space.

The detailed technical options are all documented in the help — the intention here is simpler: GeoProfiler tries to reduce the friction between thinking like a geologist and working inside a GIS. It gives you sections that behave more like the ones we sketch in the field, while still being reproducible, measurable, and connected to real data.

 

Modifications:

2026-01-06: added hyperlink to qgSurf plugin repository 

Copernicus Browser pre-calculated Sentinel 2 band indices useful for desertic regions

Continuing with this series of the naive photo-interpreter using Copernicus Browser, another interesting feature is the presence of many pre-calculated indices.

Considering geology and arid zones, a few useful pre-calculated indices are SWIR (Short-Wave Infrared Composite, bands B12, 8A and 4) and False Color (Urban, B12, 11 and 4). 

The shown examples come from the South-Eastern Libyan desert, almost to the boundary with Egypt.

As you can see, the natural color image does not evidence particular lithological differences within the outcropping structure surrounded by sands. 

 

 

Using the SWIR index, lithological boundaries and contacts are much more evident and can be used to draft a first preliminary geological interpretation, where there are volcanic and possibly sedimentary rocks.

 

 Also False Color (Urban), using bands B12, 11 and 4 differentiate between lithologies. 

 

 

Considering instead the different types of sands, NDSI (Normalised Difference Snow Index, bands  3, 11, 4 and 2) is particularly useful: the 'rivers' of younger and brighter sands around the older sand deposits are clearly evident.

 

 

Time series of NDVI within ROI with Copernicus Browser

 

The Copernicus Browser allows also to calculate time series for specific indices, for instance NDVI, by just drawing a ROI region with the pen or rectangle tool (see rectangle in image) and using the Statistical Info button.

 

The NDVI values are quite high in the last year, since this land plot is grassland, periodically cut for animal feed. It is possible to see the NDVI variations up to 5 years.

 

Timelapses creation within Copernicus Browser

Following this video: Introduction to the CDSE website and Copernicus Browser I learned that the Copernicus Browser allows to directly create timelapses animations from Sentinel data.

I tried with NDWI, that indicates the humidity degree of the surfaces, in a zone of the Piemonte (Italy) region, where I currently live. The Po river is the major river of this zone.

I set the parameters as in the image below: 

 and after a few minutes of tile selections, the following timelapse animation was created:

 

 

It's not clear to me why the final frames number is quite low, having set a 2016 to 2025 analysis interval with steps of 1 day. 

Apart from that, the result is interesting, since it potentially allows for temporal analysis of natural variations without any major effort.
 

Segments in Mojo

After implementing points in Mojo, with some assistance from ChatGPT and Phind, we now turn our attention to segments, which represent the next fundamental geometric structure preceding polylines (or simply lines).

A segment consists of a start point and an end point, rendering it an oriented entity. It is implemented as a Mojo structure, that implements a few methods, such as returning the segment start and end points, calculating the total or horizontal lenght, and so on.

To facilitate point copying when returning a segment's start or end point (via the start_point and end_point methods), the __copyinit__ method was introduced to the Point structure: 

  fn __copyinit__(inout self, existing: Point): self.coords = existing.coords 

 This method allows for efficient copying of point coordinates, which is crucial when working with segments. 

Moving on to the Segment structure, apart from incorporating several basic methods (such as dx, midx, length, horizontal_length, and others), two significant methods were added for calculating the azimuth and plunge of a Segment instance. These concepts are extensively utilized in geological contexts. 

The azimuth refers to the angle, within the horizontal plane, between the horizontal projection of the segment and the Y-axis. This angle is measured clockwise, commencing from the Y-axis, spanning from 0° to 360°. 

The plunge represents the dip of the segment in the vertical plane. It signifies the vertical angle between the horizontal plane and the segment. Its range extends from -90° to +90°, where positive values denote a downward dip, and negative values signify an upward dip. 

Both azimuth and plunge are expressed as Float64 values, but their values may remain undefined under specific circumstances. The azimuth becomes undefined when the segment is vertical or possesses zero length (i.e., when the start and end points coincide). Similarly, the plunge remains undefined when the segment's length is zero. 

To address these edge cases, both the azimuth() and plunge() methods return an Optional[Float64], with Optional being imported from the standard library's collections module. If the value is invalid, None is returned; otherwise, the valid value is returned. 

Notably, the valid return can also be encapsulated into an Optional — both approaches are considered valid. When presenting the result, the value must be extracted using the or_else method, where the input to this method serves as a placeholder for missing data. 

In the example code, a conventional no-data value frequently employed in Geographic Information Systems (GIS) was utilized by defining an alias: 

   alias NULL_ORIENTATION = -999999999.99999 

 This alias is subsequently applied as follows:

   print("Segment plunge:", segment.plunge().or_else(NULL_ORIENTATION)) 

In this scenario, when the result is valid, it will be printed; conversely, if the result is invalid, the no-data value defined in NULL_ORIENTATION will be displayed. 

By implementing these features, the Segment structure becomes more robust and versatile, capable of handling various geometric calculations essential in geological applications and beyond. 

 

The code is available at: https://gitlab.com/mauroalberti/kira

 

An example of defining a point structure in Mojo

Mojo is a relatively new language (introduced in 2022) designed for parallel computing. It is compatible with SIMD (Single Instruction Multiple Data) architectures and is intentionally interoperable with Python. SIMD allows for the application of vectorized operations, resulting in significant speedups in processing.

Currently, Mojo is available only for Ubuntu and macOS.

In this example, we will define a point geometry and implement it using Mojo.

In this context, using SIMD to store point coordinates is advantageous because it enables the application of vectorized operations to those coordinates.

To define a SIMD variable, it’s necessary to specify both the data type and the number of stored values. Here is an example of initializing a SIMD variable:

var vec = SIMD[DType.int8, 4](1, 2, 3, 4)

In Mojo, you cannot use Float32 or Float64 directly to define the variable type. Instead, you need to use DType.float32 or DType.float64 for floating-point values.

Another requirement for SIMD variables is that the number of elements must be a power of two. Therefore, it is not possible to use SIMD[DType.float64, 3] to store just the x, y, and z coordinates of a point. To meet the minimum size requirement of 4 elements, you can add a time variable (t) of type Float64 to the point coordinates. This avoids wasting memory space.

The point coordinates are stored as follows:

var coords: SIMD[DType.float64, 4]

Similar to Python, a class-like structure in Mojo is initialized using the __init__ method. To indicate that the instance is mutable, the self parameter is marked as inout:

fn __init__(inout self, x: Float64, y: Float64, z: Float64, t: Float64 = 0.0)

Since the values are stored in a SIMD variable, the elements are accessed by their index:

fn x(self) -> Float64:
	return self.coords[0]
    

To calculate the distance between two points, you cannot use the sum() or sqrt() functions directly on a SIMD variable. Instead, you must extract the scalar values by index and perform the scalar operations manually.

The square root function is imported from the math module like this:

from math import sqrt

Here is the complete implementation of the point structure:

from math import sqrt

struct Point :

    var coords: SIMD[DType.float64, 4]

    fn __init__(inout self, x: Float64, y: Float64, z: Float64, t: Float64 = 0.0):

        self.coords = SIMD[DType.float64, 4](x, y, z, t)

    fn x(self) -> Float64:
        return self.coords[0]

    fn y(self) -> Float64:
        return self.coords[1]

    fn z(self) -> Float64:
        return self.coords[2]

    fn t(self) -> Float64:
        return self.coords[3]

    fn distance_to(self, other: Point) -> Float64:

        var delta = self.coords - other.coords
        var sum_of_squares = delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2];
        return sqrt(sum_of_squares[0])
    
    

And here’s an example of how to use this structure in a main function:

fn main():

    var pt1 = Point(1.0, 2.0, 3.0, 1222.34)
    var pt2 = Point(2.0, 3.0, 3.0, 22.34)
    print(pt1.distance_to(pt2))
      

You can also view this example on GitLab:

Mojo Point Structure Example

geogst, a new Python module for Structural Geology

 

geogst is a new Python module for structural geology available on PyPi, making it easily installable via the classic command pip install geogst (or python -m pip install geogst).

This module was primarily developed to facilitate the creation of stereonets for geological data, determining the intersection between geological planes and topographic surfaces (expressed through Digital Elevation Models, DEMs), and generating the skeletons of geological profiles. Among its main features, the module allows for the calculation of topographic profiles, adding geological attitudes, and determining the intersections of geological traces with these profiles.

 

Additionally, geogst includes example geospatial datasets that can be used to explore its functionalities, such as creating profiles directly within Jupyter Notebooks.

Here are a few examples of Jupyter Notebooks for creating geological profiles:

• Geologic profiles map animations - the Timpa San Lorenzo structure (Southern Italy):
  Geologic profiles with animation - Timpa San Lorenzo

• Geologic profiles - the 1997 seismic sequence in Colfiorito (Central Italy):
  Seismological profiles - Colfiorito

This module will form the foundation for the modules in the qgSurf plugin for QGIS. Currently, qgSurf directly includes geogst as a submodule, so no separate installation is required. In future versions of qgSurf (initially experimental), the geogst module will be automatically installed if it is not already present.

Since QGIS does not allow the upload of compiled code in Python modules, a key advantage of using geogst as a separate module in qgSurf is the potential to incorporate compiled code in Fortran, C++, and Rust, significantly improving the speed and efficiency of these tools.

Conclusion and call to action: If you are a geologist or a developer working with geological data, we invite you to explore geogst and contribute to its development. Your experiences and feedback are crucial to improving and growing the community that uses ggSurf.

Mojo

A new AI-oriented language is on the scene: Mojo, which aims to become a superset of Python. It follows Python semantics but, thanks to recent and advanced compilation techniques, might achieve speedups as large as 100,000 times or more compared to standard Python.

Mojo was created by Chris Lattner, the creator of LLVM, Clang, Swift, Swift for TensorFlow and MLIR.

Currently, Mojo is only supported on Linux and macOS.

The source code is hosted at GitHub, from where you can clone the repository locally. Within the downloaded repository, the 'examples/notebooks' directory contains Jupyter Notebooks that you can run using Mojo (provided that Mojo and the Jupyter plugin are installed).

A very brief introduction to the Mojo language is provided in HelloMojo.ipynb. Another notebook, Matmul.ipynb, implements a matrix multiplication example in both Python and Mojo. The example showcases how optimization in Mojo can lead to a remarkable speedup, reportedly around 455,127x compared to Python. It's important to note that achieving such optimized Mojo code requires considerable effort and a solid understanding of Mojo.

That said, I tried running the code in that notebook on my old HP laptop, which lacks a GPU, and has the following specifications:

inxi -Fxxxzr

System:
  Kernel: 5.15.0-118-generic x86_64 bits: 64 compiler: gcc v: 11.4.0
  Desktop: Xfce 4.18.1 tk: Gtk 3.24.33 info: xfce4-panel wm: xfwm 4.18.0
  vt: 7 dm: LightDM 1.30.0 Distro: Linux Mint 21.3 Virginia
  base: Ubuntu 22.04 jammy
Machine:
  Type: Laptop System: HP product: HP Laptop 15-bs0xx v: Type1ProductConfigId
  serial: <superuser required> Chassis: type: 10 serial: <superuser required>
  Mobo: HP model: 832B v: 23.37 serial: <superuser required> UEFI: Insyde
  v: F.21 date: 07/04/2017
CPU:
  Info: dual core model: Intel Core i5-7200U bits: 64 type: MT MCP
  smt: enabled arch: Amber/Kaby Lake note: check rev: 9 cache: L1: 128 KiB
  L2: 512 KiB L3: 3 MiB
  Speed (MHz): avg: 1134 high: 1245 min/max: 400/3100 cores: 1: 1040
  2: 1097 3: 1245 4: 1157 bogomips: 21599

Even though the speedup I achieved was an order of magnitude lower than the original 455,127x reported in the notebook, it was still impressive: 10,032x! It’s likely that using a newer machine, possibly equipped with a GPU, I could achieve speedups on the order of 100,000x.

As reported, tools like Cython or Numba achieve speedups typically in the range of 10-100x. 

Moreover, the Mojo code ran without any issue on my old laptop running Linux Mint.

So Mojo is a very interesting language, even in its infancy.

Other AI-oriented languages to explore, as described in the insightful post by James Thomason in VentureBeat, "Mojo Rising: The resurgence of AI-first programming languages" include Bend and JAX.

 

Note: the text was checked in ChatGPT.

 

 

SAVI index and agricultural landscapes

SAVI is the acronym for Soil Adjusted Vegetation Index that, as the name suggests, it is one of the many indices for the calculation of the vegetation index.

The general formula for the SAVI index is (indexdatabase.de):

     800 nm - 670 nm
  ________________________  ( 1 + L)

   800 nm + 670 nm + L

   
where L is between -0.9 and 1.6.

In the case of Sentinel 2 data, SAVI uses two bands, B8 and B4, that have the maximum available resolution for Sentinel 2 images, i.e., 10 meters, so the results have a quite high spatial resolution.

B8 is in the NIR (842 nm) while B4 registers the red-NIR transition (665 nm).

For Sentinel 2, the index is therefore calculated as:

     B8 - B4
  _______________   ( 1 + L)

   B8 + B4 + L
   
The default value for L in indexdatabase is 0.5, and in fact testing for a few values (-0.9, -0.5, 0.0, 0.5, 1.0, 1.5) in the prescribed range, the most contrasted results are obtained for L = 0.5.

One information that is obtained from applying such an index to intensive farming landscapes is that it allows to delineate well the agricultural field limits, to compare the texture and the "spatial styles" of farming across different areas, for instances separated by rivers that may have acted as administrative or political borders.

In the following images, I present some examples of textures and styles of limits between fields for different areas in the Piedmont region, derived from a Sentinel 2B image, acquired on 2023/04/04 at 10:25. The SAVI index was calculated using the SNAP software by ESA.

The used L ("soil brightness correction factor") value is 0.5, chosen after some experimentation.

The figures below are generated from within SNAP, using a 8-4-3 composite for infrared view (upper window), while the SAVI index band is represented in the lower window. In the 843 composite, vegetation is red while rivers and roads are black.

Strong agricultural landscapes contrast between the western side (left) and eastern  side (right) of the Sesia river (central, with a North-South orientation). The western landscape is characterized by smaller sizes than in the eastern part, and is also more regular as geometries. Possibly it is the result of different political systems between the two sides in the previous centuries.

 Strong East-West differences in field sizes, orientations and textures without a clear separating physical boundary (such as the Sesia river in the previous example). It could be related to different grown agricultural products between the two zones (?rice to the East), possibly also related to distinct administrative/political situations.  

Irregular field borders probably due to the limiting rivers and their location changes with time.

 

 

Rayleigh Correction for Sentinel data using SNAP

Removing atmospheric effects from satellite images is tricky, since the theoretical foundations are not completed understood, it is difficult to model the physical conditions, and also because we generally lack many information regarding the atmospheric condition at the time and location of satellite image acquisition.

One procedure that attempts to partially remove atmospheric effects is implemented in the SNAP software application, that is released by ESA and mainly devoted to the processing of Sentinel data, both optical and radar.

This method is named "Rayleigh Correction" and attempts to mitigate gaseous and Rayleigh effects. It does however do not consider aerosol effects.

In the screenshot below you see how to open the method window.

The GUI allows to select the bands to be corrected (B2-4 in the example below), besides other parameters (e.g., sea-level pressure). For this example the output resolution was changed to 10 meters, to match that of input bands (B2-4). The other parameters were left unchanged.

The two windows below display a an example of Sentinel 2 data (S2A MSI, 2017-05-27) pre- and post- Rayleigh correction. The displayed zone is comprised between Chivasso to the West and Verolengo to the East, and the Po River to the South and the Canale Cavour to the North (Piedmont, Italy). 

In the upper window the original image is displayed as a 4-3-2 composite. In the lower window the same zone is represented using a 4-3-2 composite of Rayleigh-corrected bands.

The original bands are darker and with a low contrast, while the corrected bands present a larger chromatic variance. Note that the lightest pixels in the original band tend to saturate to white, somewhat loosing information.

Water, vegetation and urban constructions are more vivid and looks like more "natural", with more chromatic variations than in the original bands.

Quick calculation of band indices from satellite images with GDAL

Band indices of satellite images, such as the NDVI, can be routinely calculated from within applications, such as ENVI, SNAP and so on.

Sometimes it is preferable to calculate band indices using a programming language: it could be faster than using applications, we could not dispose of the required applications or we need to automate a set of processing by using a programming language.

Having GDAL installed, it is quick and easy to calculate band indices (and a lot more) directly from the prompt/shell. I was drawn to this useful solution by the answer provided by

The used command is gdal_calc.py and we can specify one, two or more input bands to perform calculations, for instance, for calculating the Normalized Difference Water Index (NDWI) from a spatial subset of a Sentinel 2 image saved from SNAP as BEAM-Dimap, we use:

gdal_calc.py -A B3.img -B B8.img --outfile=ndwi.tif --calc="((B-A)/(B+A))"

Obviously the shell working directory is the folder containing the source image files. Morevover we do not have to define the used band within the input .img file since there is only one band.
The complete description of the command is available at the page 'gdal_calc.py'.


An example of calculation of the Normalized Difference Water Index is here presented. Three Sentinel 2 SLC image of the Po river in the Piedmont area between Chivasso and Casale Monferrato (Northern Italy) were spatially subset in the range lat 45.1 -> 45.2 and long 7.8 -> 8.6 from within SNAP (Fig. 1).

 Fig. 1. Map of the studied area in the central Piedmont region (Northern Italy). Basemap: Google Earth. Created with QGIS.

The NDWI was calculated with the command line as listed above:

gdal_calc.py -A B3.img -B B8.img --outfile=ndwi.tif --calc="((B-A)/(B+A))" 

From the analysis of the histogram of the 2017 image, a value of 10 was chosen of a threshold for masking the non-water pixels:

gdal_calc.py -A ndwi.tif --outfile=ndwi_masked.tif --calc="A>10"

The result is represented in Fig. 2A-C, created with QGIS. The result quality is more than acceptable (for the 2017 image the spotted fields to the North of the Po river could possibly correspond to rice fields, but it has to be better investigated).

A - 2017-05-27

B - 2021-03-30

C- 2023-04-04

 Fig. 2. NDWI > 10 for 2017 (a), 2021 (b) and 2023 (c). Basemap: Google Earth. Created with QGIS.

By comparing the results at different times, it is possible to study the temporal evolution of the rivers (in this case, the Po river, Orco torrent and Canale Cavour near the city of Chivasso (Fig. 3).

 Fig. 2. NDWI > 10 for 2017 (red), 2021 (green) and 2023 (blue) in correspondence of the city of Chivasso. The water bodies are the Po river, the Orco torrent and the artificial Canale Cavour. Basemap: Google Earth. Created with QGIS.

 
If we want to automate the task, when we have for instance a lot of input images or we have to perform routinely the same task on new datasets, the run command can be created on-the-fly and executed from within a script created with languages such as Go, Rust, Ocaml and so on, perhaps even concurrently.