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.