Skip to content

WIP. Add stereonet module for plotting structural geology stereonets (Schmidt/Wulff) - #9178

Open
Esteban82 wants to merge 20 commits into
masterfrom
stereonet-module
Open

Esteban82 wants to merge 20 commits into
masterfrom
stereonet-module

Conversation

@Esteban82

Copy link
Copy Markdown
Member

Opening this as a draft/WIP specifically to get feedback before going further. @rhum1-geo, your suggestions (density contouring, Fisher stats, notation auto-detection, richer -S/vector options) are noted and intentionally not in this first PR — the plan is to land plotting first, then follow up with statistics in a later PR once the plotting side is settled.

What's implemented:

-T[d|l|p] — planes as strike/dip (right-hand rule) or dip-direction/dip, or lines as trend/plunge
-W/-S/-L/-G — cyclographic traces and poles (or lines), styled like plot
-A/-B — azimuth ring and the net's own grid mesh, both opt-in
-T+u — upper hemisphere

Tested with (no automated tests yet — these are the scripts I manually ran and checked; anyone from #9156 is welcome to run them too):

0_Empty_Nets.sh — empty Schmidt and Wulff nets, with title, grid, and azimuth ring, no data
gmt set MAP_GRID_PEN_PRIMARY 0.25p,gray

gmt begin "0_Empty_Net_Schdmit_(A)" png
  echo /dev/null | gmt stereonet -JA10c -B+t"Schdmit" -B -A
gmt end

gmt begin "0_Empty_Net_Wullf_(S)" png
   echo /dev/null | gmt stereonet -JS10c -B+t"Wullf" -B -A
gmt end
0_Empty_Net_Schdmit_(A) 0_Empty_Net_Wullf_(S)
1_Plot_Data.sh — a plane + its pole from strike/dip (RHR) input, styled independently
# RHR. Rhumb Dip
cat << EOF > data.txt
270 30
180 10
EOF

gmt begin "1_Plot_Data" png
  gmt set MAP_GRID_PEN_PRIMARY 0.25p,gray
  gmt stereonet data.txt -JS10c -Sc0.5c -Gorange -L1p,blue -l"Pole" -B -A
  gmt stereonet data.txt -JS10c -W1p,red -l"Plane" 
gmt end
rm -f data.txt gmt.*
1_Plot_Data
2_Plot_Data_Inverted.sh — same data with columns swapped; -: and -i1,0 should reproduce script 1 exactly
cat << EOF > data.txt
# Dip Rhumb
30 270
10 180
EOF

gmt begin "2_Plot_Data_Inverted" png
  gmt set MAP_GRID_PEN_PRIMARY 0.25p,gray
  gmt stereonet data.txt -: -JS10c -Sc0.5c -Gorange -L1p,blue -l"Pole" -B -A
  gmt stereonet data.txt -i1,0 -JS10c -W1p,red -l"Plane"
gmt end
rm -f data.txt gmt.*
3_Plot_Planes_Conventions.sh — the same plane entered as -Tp (strike/dip, RHR) and -Td (dip-direction/dip); traces should coincide
cat << EOF > data1.txt
# Rhumb Dip
270 30
180 10
EOF

# Dip-Direction Dip
cat << EOF > data2.txt
0 30
270 10
EOF

gmt begin 3_Plot_Planes_Conventions png
	gmt set MAP_GRID_PEN_PRIMARY 0.25p,gray
	gmt stereonet data1.txt -JA12c -W2p,red    -Sc0.4c -Gred    -Tp -l"-Tp: RHR" -B -A
	gmt stereonet data2.txt        -W1p,yellow -Sc0.2c -Gyellow -Td -l"-Td: Dip direction/Dip"
gmt end
rm data*.txt
3_Plot_Planes_Conventions
4_Plot_Poles_Conventions.sh-Tp and -Td poles for the same two planes should coincide (red under blue); -Tl reuses the same numbers as a line's trend/plunge instead, which is a different quantity and lands elsewhere on purpose
# RHR: Rhumb Dip
cat << EOF > data1.txt
270 30
180 10
EOF

# Dip-Direction Dip
cat << EOF > data2.txt
0 30
270 10
EOF

# Trend / Plunge
cat << EOF > data3.txt
0 30
270 10
EOF

gmt begin 4_Plot_Poles_Conventions png
	gmt set MAP_GRID_PEN_PRIMARY 0.25p,gray
	gmt stereonet data1.txt -Sc0.4c -Gred   -Tp -l"-Tp: RHR" -B
	gmt stereonet data2.txt -Sc0.2c -Gblue  -Td -l"-Td: Dip direction/Dip"
	gmt stereonet data3.txt -Sc0.2c -Ggreen -Tl -l"Lines: Trend/Plunge"   
gmt end
rm data*.txt
4_Plot_Poles_Conventions
5_Annotations.sh — fractional azimuth interval (-A22.5) with a title and -Vi
gmt set MAP_GRID_PEN_PRIMARY 0.25p,gray
gmt begin "5_Annotations" png
  echo /dev/null | gmt stereonet -JA10c -B+t"Schdmit" -A22.5 -Vi -B
gmt end
5_Annotations

Model/effort: Claude Opus 5 for design and code review, Claude Sonnet 5 for iteration and testing.

Related to #9156.

Review fixes on the stereonet module:

- Fix a segfault: a bare -W (no pen argument) together with -Tl left
  W.string NULL while W.active was already true, so aliasing the pen onto
  -L called strdup(NULL).
- Reject dips/plunges outside 0-90 and rakes outside 0-180. Such values
  projected onto the far hemisphere, where they are clipped away, so the
  module quietly produced a figure with missing data instead of an error.
- Accept a negative rake as the usual shorthand for a rake measured from
  the opposite end of the strike, folding it into the 0-180 range (-25 is
  read as 155), matching the convention used by mplstereonet. The docs
  previously claimed negatives were rejected, which was never enforced.
- Restore test/geology/minimal.sh as a classic-mode test producing $ps so
  the test harness can compare it, and drop the minimal.png/minimal.txt
  artifacts that were committed by accident.
New modules must add themselves to four separate hardcoded name lists in
gmt_modern.c, or their modern name silently misbehaves in ways that only
show up once you specifically go looking (this is not mentioned in
devdocs/custom_supplements.rst):

- gmt_current_name(): without this, every GMT_Report message (errors,
  warnings, info) prints the classic name "psstereonet" even when the
  user typed "stereonet" and is in modern mode.
- gmtlib_get_active_name(): same classic->modern translation, used by a
  different set of callers.
- gmtlib_is_modern_name(): without this, "gmt stereonet -^" (or any bare
  invocation outside a session, e.g. querying usage) fails outright with
  "Shared GMT module not found: stereonet", because the modern-mode
  leniency that lets a bare modern name grab usage/purpose text before a
  session exists never triggers.

Found by testing "gmt stereonet ..." against the exact same scenarios for
core module ternary and supplement modules polar/coupe (same pattern,
already correctly registered) and diffing the behavior.
@Esteban82 Esteban82 added add-changelog Add PR to the changelog new module PR that implements a new module AI-assisted All (or most) of the code was written by Artificial Intelligence. labels Sep 4, 2026
Esteban82 and others added 3 commits September 7, 2026 13:56
… stereonet

Implements the density-contour suggestion from issue #9156 (Rom1's comment):
-D[e|k|s][+c[<cpt>]][+i<interval>][+p<pen>][+s<sigma>] contours the density of
plotted poles/lines per Kamb (1959) and Vollmer's (1995) modified/exponential
method, or the traditional Schmidt 1% method. Math and grid-domain derivation
checked against Vollmer (1995) and independently against mplstereonet, then
reviewed (an initial draft wrongly dropped the fabs() fold and mis-normalized
Schmidt; both fixed).

Contour lines are traced with a small marching-squares pass over an internal
grid and drawn via the existing psxy/virtual-dataset path, not grdcontour:
calling grdcontour via GMT_Call_Module under -O -K from inside an azimuthal
projection's active plot reproducibly mis-scaled the output (confirmed with a
synthetic cluster across many -R/-X/-Y combinations, real and virtual grids).
+c (CPT fill between contours) is parsed but rejected with a clear error for
now, since proper fill needs polygon extraction, not just lines.

Tested with:
  gmt psstereonet cluster.txt -JA10c -B -Sc0.08c -Gblack -De

Fixes #9156

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

add-changelog Add PR to the changelog AI-assisted All (or most) of the code was written by Artificial Intelligence. new module PR that implements a new module

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant