diff --git a/doc/rst/source/modules-classic.rst b/doc/rst/source/modules-classic.rst index 89903753327..971b3741530 100644 --- a/doc/rst/source/modules-classic.rst +++ b/doc/rst/source/modules-classic.rst @@ -127,6 +127,7 @@ All modules are requested via a call to the :doc:`gmt` program. supplements/geodesy/earthtide supplements/geodesy/gpsgridder supplements/geodesy/psvelo + supplements/geology/psstereonet supplements/gsfml/fzanalyzer supplements/gsfml/fzblender supplements/gsfml/fzinformer @@ -298,6 +299,7 @@ Supplemental Modules - :doc:`/supplements/geodesy/earthtide` - :doc:`/supplements/geodesy/gpsgridder` - :doc:`/supplements/geodesy/psvelo` + - :doc:`/supplements/geology/psstereonet` - :doc:`/supplements/gsfml/fzanalyzer` - :doc:`/supplements/gsfml/fzblender` - :doc:`/supplements/gsfml/fzinformer` @@ -621,6 +623,13 @@ geodesy | :doc:`/supplements/geodesy/psvelo` | |psvelo_purpose| | +--------------------------------------------+----------------------+ +geology +------- + ++---------------------------------------------+-----------------------+ +| :doc:`/supplements/geology/psstereonet` | |psstereonet_purpose| | ++---------------------------------------------+-----------------------+ + GSFML ----- diff --git a/doc/rst/source/modules.rst b/doc/rst/source/modules.rst index 105054a546d..590ea767c31 100644 --- a/doc/rst/source/modules.rst +++ b/doc/rst/source/modules.rst @@ -132,6 +132,7 @@ All modules are requested via a call to the :doc:`gmt` program. supplements/geodesy/earthtide supplements/geodesy/gpsgridder supplements/geodesy/velo + supplements/geology/stereonet supplements/gshhg/gshhg supplements/img/img2grd supplements/mgd77/mgd77convert @@ -309,6 +310,7 @@ Supplemental Modules - :doc:`/supplements/geodesy/earthtide` - :doc:`/supplements/geodesy/gpsgridder` - :doc:`/supplements/geodesy/velo` + - :doc:`/supplements/geology/stereonet` - :doc:`/supplements/gsfml/fzanalyzer` - :doc:`/supplements/gsfml/fzblender` - :doc:`/supplements/gsfml/fzinformer` @@ -652,6 +654,13 @@ geodesy | :doc:`/supplements/geodesy/velo` | |velo_purpose| | +--------------------------------------------+----------------------+ +geology +------- + ++--------------------------------------------+----------------------+ +| :doc:`/supplements/geology/stereonet` | |stereonet_purpose| | ++--------------------------------------------+----------------------+ + GSFML ----- diff --git a/doc/rst/source/reference/map-projections.rst b/doc/rst/source/reference/map-projections.rst index a54b9c1eb66..2de45d7f1a8 100644 --- a/doc/rst/source/reference/map-projections.rst +++ b/doc/rst/source/reference/map-projections.rst @@ -290,6 +290,9 @@ at 0/0) is known as the *equal-area* (Schmidt) stereonet and used for plotting fold axes, fault planes, and the like. An *equal-angle* (Wulff) stereonet can be obtained by using the stereographic projection (discussed later). The stereonets produced by these two projections appear below. +The :doc:`/supplements/geology/stereonet` module builds on these two projections and +converts structural measurements (strike/dip or trend/plunge) into the coordinates they +expect. .. _GMT_stereonets: diff --git a/doc/rst/source/supplements/geology/psstereonet.rst b/doc/rst/source/supplements/geology/psstereonet.rst new file mode 100644 index 00000000000..5715103f573 --- /dev/null +++ b/doc/rst/source/supplements/geology/psstereonet.rst @@ -0,0 +1,69 @@ +.. index:: ! psstereonet +.. include:: ../module_supplements_purpose.rst_ + +*********** +psstereonet +*********** + +|psstereonet_purpose| + +Synopsis +-------- + +.. include:: ../../common_SYN_OPTs.rst_ + +**gmt psstereonet** [ *table* ] +[ |-J|\ **A**\|\ **S**\ *width* ] +[ |-A|\ [*annot*\ [/*tick*]] ] +[ |SYN_OPT-B| ] +[ |-G|\ *fill* ] +[ |-K| ] +[ |-L|\ *pen* ] +[ |-O| ] [ |-P| ] +[ |-S|\ *symbol*\ [*size*] ] +[ |-T|\ [**d**\|\ **l**\|\ **p**][**+u**] ] +[ |SYN_OPT-U| ] +[ |SYN_OPT-V| ] +[ |-W|\ *pen* ] +[ |SYN_OPT-X| ] +[ |SYN_OPT-Y| ] +[ |SYN_OPT-bi| ] +[ |SYN_OPT-di| ] +[ |SYN_OPT-e| ] +[ |SYN_OPT-f| ] +[ |SYN_OPT-g| ] +[ |SYN_OPT-h| ] +[ |SYN_OPT-i| ] +[ |SYN_OPT-p| ] +[ |SYN_OPT-qi| ] +[ |SYN_OPT-s| ] +[ |SYN_OPT-t| ] +[ |SYN_OPT-:| ] +[ |SYN_OPT--| ] + +|No-spaces| + +.. include:: stereonet.rst + :start-after: .. module_common_begins + :end-before: .. module_common_ends + +.. include:: ../../common_classic.rst_ + +Examples +-------- + +.. include:: ../../explain_example.rst_ + +To plot eight fault planes given as *strike dip* on a 12-centimeter-wide Schmidt net, +drawing the cyclographic traces in red and the poles as blue crosses, try:: + + gmt psstereonet faults.txt -JA12c -B -W1p,red -Sx0.3c -L1p,blue -P > faults.ps + +See Also +-------- + +:doc:`gmt `, :doc:`gmt.conf `, +:doc:`psbasemap `, +:doc:`psxy `, +:doc:`psrose `, +:doc:`/supplements/seis/pspolar` diff --git a/doc/rst/source/supplements/geology/stereonet.rst b/doc/rst/source/supplements/geology/stereonet.rst new file mode 100644 index 00000000000..feeb1f778e9 --- /dev/null +++ b/doc/rst/source/supplements/geology/stereonet.rst @@ -0,0 +1,273 @@ +.. index:: ! stereonet +.. include:: ../module_supplements_purpose.rst_ + +********* +stereonet +********* + +|stereonet_purpose| + +Synopsis +-------- + +.. include:: ../../common_SYN_OPTs.rst_ + +**gmt stereonet** [ *table* ] +[ |-J|\ **A**\|\ **S**\ *width* ] +[ |-A|\ [*annot*\ [/*tick*]] ] +[ |SYN_OPT-B| ] +[ |-G|\ *fill* ] +[ |-L|\ *pen* ] +[ |-S|\ *symbol*\ [*size*] ] +[ |-T|\ [**d**\|\ **l**\|\ **p**][**+u**] ] +[ |SYN_OPT-U| ] +[ |SYN_OPT-V| ] +[ |-W|\ *pen* ] +[ |SYN_OPT-X| ] +[ |SYN_OPT-Y| ] +[ |SYN_OPT-bi| ] +[ |SYN_OPT-c| ] +[ |SYN_OPT-di| ] +[ |SYN_OPT-e| ] +[ |SYN_OPT-f| ] +[ |SYN_OPT-g| ] +[ |SYN_OPT-h| ] +[ |SYN_OPT-i| ] +[ |SYN_OPT-l| ] +[ |SYN_OPT-p| ] +[ |SYN_OPT-qi| ] +[ |SYN_OPT-s| ] +[ |SYN_OPT-t| ] +[ |SYN_OPT-:| ] +[ |SYN_OPT--| ] + +|No-spaces| + +.. module_common_begins + +Description +----------- + +Reads pairs of angles in degrees from *table* [or standard +input] and plots them on a **stereonet**, the diagram that structural geologists use to +display the orientation of planes and lines. Depending on |-T| the two angles are read +as the strike and dip of a plane (using the right-hand rule, so that the plane dips to +the right of the strike direction), as the dip direction and dip of a plane, or as the +trend and plunge of a line. + +For each plane we can draw its **cyclographic trace**, i.e., the great circle where the +plane cuts the hemisphere (selected with |-W|), and the **pole** to the plane, i.e., the +point where the normal to the plane pierces the hemisphere (selected with |-S|). For +each line we only draw the point where it pierces the hemisphere. + +A stereonet is not a separate map projection: it is simply one hemisphere of a unit +sphere seen from above, centered on the nadir. We therefore use one of the two standard +GMT azimuthal projections centered on 0/0, whose default 90-degree horizon is exactly +one hemisphere: |-J|\ **A** (Lambert azimuthal equal-area) gives the equal-area net that +is also known as a **Schmidt net**, while |-J|\ **S** (stereographic) gives the +equal-angle net that is also known as a **Wulff net** (see Figure +:ref:`Stereonets `). With that setup north is up, so azimuth is measured +clockwise from the top of the plot, the center of the net is a vertical (90-degree plunge) +direction, and the perimeter is horizontal. + +Required Arguments +------------------ + +.. |Add_intables| unicode:: 0x20 .. just an invisible code +.. include:: ../../explain_intables.rst_ + +Optional Arguments +------------------ + +.. _-A: + +**-A**\ [*annot*\ [/*tick*]] + Annotate the azimuth around the perimeter of the net every *annot* degrees, with tick + marks every *tick* degrees [30/10 if **-A** is given with no argument]. If only *annot* + is given then *tick* defaults to one third of *annot*. Without **-A** no azimuth ring + is drawn at all; **-A0** is the same as omitting **-A**. + +.. _-B: + +.. include:: ../../explain_-B.rst_ + :start-after: **Syntax** + :end-before: **Description** + +The gridlines requested via |-B| *are* the net: the meridians are the cyclographic traces +of planes striking N-S whose dips step by the grid interval, and the parallels are the +small circles of constant plunge. Since the net has no meaningful longitude or latitude +annotations you will normally only ask for gridlines, and the two-level mesh that a +stereonet traditionally shows is obtained by giving both a primary and a secondary +interval, e.g., **-Bpg10 -Bsg30**. Without |-B| no frame at all is drawn, not even the +perimeter of the net; give a bare **-B** for the classic two-level mesh [-Bpg10 -Bsg30]. +A |-B| that only carries frame settings, such as a **-B+t**\ *title*, gets that same mesh, +so you can title a default net without spelling out the intervals. +Add, e.g., ``--MAP_GRID_PEN_PRIMARY=0.25p,gray`` to make the fine mesh recede into the +background. + +.. _-G: + +**-G**\ *fill* :ref:`(more ...) <-Gfill_attrib>` + Select color or pattern for filling the symbols [Default is no fill]. + +.. _-J: + +**-JA**\|\ **S**\ *width* + Select the type of net, where *width* is its diameter: + + - **A** - Lambert azimuthal equal-area, i.e., a Schmidt net [Default]. + - **S** - Stereographic (equal-angle), i.e., a Wulff net. + + Both projections use the default 90-degree horizon. The net is always centered on 0/0, + so give the width only, without a center: any other projection, or an explicit center, + is an error. If |-J| is skipped we inherit the net of an earlier **stereonet** in the + same figure, else we default to **-JA**\ 15c. + +.. _-L: + +**-L**\ *pen* :ref:`(more ...) <-Wpen_attrib>` + Set the pen used to outline the symbols selected with |-S|. + +.. _-S: + +**-S**\ *symbol*\ [*size*] + Plot the pole to each plane (or the line itself if **-Tl**) using this symbol; see + :doc:`plot ` for the available symbol codes [**-Sc**\ 0.15c]. Without |-S| no + symbols are plotted for planes. + +.. _-T: + +**-T**\ [**d**\|\ **l**\|\ **p**][**+u**] + Select what the two input angles mean: + + - **d** - Planes given as dip direction and dip. + - **l** - Lines given as trend and plunge. + - **p** - Planes given as strike and dip, with the strike following the right-hand + rule so that the plane dips to the right of the strike direction [Default]. + + Optionally, append modifier: + + - **+u** - Plot the data on the upper hemisphere [Default is the lower hemisphere, + which is the convention in structural geology]. + +.. |Add_-U| replace:: |Add_-U_links| +.. include:: ../../explain_-U.rst_ + :start-after: **Syntax** + :end-before: **Description** + +.. |Add_-V| replace:: |Add_-V_links| +.. include:: ../../explain_-V.rst_ + :start-after: **Syntax** + :end-before: **Description** + +.. _-W: + +**-W**\ *pen* :ref:`(more ...) <-Wpen_attrib>` + Set the pen used to draw the cyclographic traces of the planes. Ignored if **-Tl**, + since a line has no trace; use |-L| to outline the symbols instead. + +.. |Add_-XY| replace:: |Add_-XY_links| +.. include:: ../../explain_-XY.rst_ + :start-after: **Syntax** + :end-before: **Description** + +.. |Add_-bi| replace:: [Default is 2 input columns]. +.. include:: ../../explain_-bi.rst_ + +.. include:: ../../explain_-c.rst_ + +.. |Add_-di| unicode:: 0x20 .. just an invisible code +.. include:: ../../explain_-di.rst_ + +.. |Add_-e| unicode:: 0x20 .. just an invisible code +.. include:: ../../explain_-e.rst_ + +.. |Add_-f| unicode:: 0x20 .. just an invisible code +.. include:: ../../explain_-f.rst_ + +.. |Add_-g| unicode:: 0x20 .. just an invisible code +.. include:: ../../explain_-g.rst_ + +.. |Add_-h| unicode:: 0x20 .. just an invisible code +.. include:: ../../explain_-h.rst_ + +.. include:: ../../explain_-icols.rst_ + +.. include:: ../../explain_-l.rst_ + +.. |Add_perspective| unicode:: 0x20 .. just an invisible code +.. include:: ../../explain_perspective.rst_ + +.. include:: ../../explain_-qi.rst_ + +.. include:: ../../explain_-s.rst_ + +.. include:: ../../explain_-t.rst_ + +.. include:: ../../explain_colon.rst_ + +.. include:: ../../explain_help.rst_ + +Notes +----- + +#. Repeat **-l** to label both items: the first label goes with the cyclographic traces and + the second with the symbols. +#. To draw an empty net, give no input records, e.g., ``gmt stereonet -JA12c -B < /dev/null``. + Without **-B** no frame at all is drawn, not even the perimeter of the net. +#. Field notebooks often list a name before the two angles; use, e.g., **-i**\ 1,2 to skip + such a leading column. +#. Dips (and plunges) must lie in the 0-90 range, as that is the only physically meaningful + value. An out-of-range angle is reported as an error rather than plotted, since it would + otherwise project onto the far hemisphere and be silently clipped away. +#. There is no **-R** option: a stereonet always covers a full hemisphere, and giving one + is rejected as an error rather than silently ignored. + +.. module_common_ends + +Examples +-------- + +.. include:: ../../explain_example.rst_ + +.. include:: ../../oneliner_info.rst_ + +To plot eight fault planes given as *strike dip* on a 12-centimeter-wide Schmidt net, +drawing the cyclographic traces in red and the poles as blue crosses, try:: + + cat << EOF > faults.txt + 90 30 + 180 45 + 270 60 + 0 15 + 30 45 + 120 48 + 225 27 + 350 80 + EOF + gmt begin faults + gmt stereonet faults.txt -JA12c -B -W1p,red -Sx0.3c -L1p,blue \ + -l"Fault plane" -l"Pole" --MAP_GRID_PEN_PRIMARY=0.25p,gray + gmt end show + +To plot the poles to bedding, measured as *trend plunge*, as red circles on a 10-centimeter +Wulff net without the azimuth ring, try:: + + gmt begin bedding + gmt stereonet bedding.txt -JS10c -B -Tl -A0 -Sc0.2c -Gred + gmt end show + +References +---------- + +Lisle, R. J., and P. R. Leyshon, 2004, *Stereographic Projection Techniques for Geologists +and Civil Engineers*, 2nd edition, Cambridge University Press. + +See Also +-------- + +:doc:`gmt `, :doc:`gmt.conf `, +:doc:`basemap `, +:doc:`plot `, +:doc:`rose `, +:doc:`/supplements/seis/polar` diff --git a/doc/rst/source/supplements/module_supplements_purpose.rst_ b/doc/rst/source/supplements/module_supplements_purpose.rst_ index 630c84088b0..7c8232b869d 100644 --- a/doc/rst/source/supplements/module_supplements_purpose.rst_ +++ b/doc/rst/source/supplements/module_supplements_purpose.rst_ @@ -8,6 +8,10 @@ .. |psvelo_purpose| replace:: Plot velocity vectors, crosses, anisotropy bars and wedges +.. |stereonet_purpose| replace:: Plot structural geology data on a stereonet + +.. |psstereonet_purpose| replace:: Plot structural geology data on a stereonet + .. |fzanalyzer_purpose| replace:: Analysis of fracture zones using crossing profiles .. |fzblender_purpose| replace:: Produce a smooth blended FZ trace diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 781ba594265..cbff04ec821 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -728,7 +728,7 @@ if (BUILD_SUPPLEMENTS) # supplement directories to be built as supplement libraries # SUPPL_EXTRA_DIRS are directories for custom supplements. # See cmake/ConfigUserAdvancedTemplate.cmake for details. - set (GMT_SUPPL_DIRS geodesy gsfml gshhg img mgd77 potential segy seis spotter x2sys windbarbs nswing ${SUPPL_EXTRA_DIRS}) + set (GMT_SUPPL_DIRS geodesy geology gsfml gshhg img mgd77 potential segy seis spotter x2sys windbarbs nswing ${SUPPL_EXTRA_DIRS}) # process supplement directories and collect information for each supplement libraries set (GMT_SUPPL_LIBRARIES) # empty the library list diff --git a/src/geology/CMakeLists.txt b/src/geology/CMakeLists.txt new file mode 100644 index 00000000000..31bbd8065e5 --- /dev/null +++ b/src/geology/CMakeLists.txt @@ -0,0 +1,32 @@ +# +# +# Copyright (c) 1991-2026 by the GMT Team (https://www.generic-mapping-tools.org/team.html) +# See LICENSE.TXT file for copying and redistribution conditions. +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation; version 3 or any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# Contact info: www.generic-mapping-tools.org +#------------------------------------------------------------------------------- +# +# CMake settings for supplemental package: +# +# SUPPL_NAME: name of the supplemental package +# SUPPL_HEADERS: header files. Will be installed if BUILD_DEVELOPER is TRUE +# SUPPL_PROGS_SRCS: list of C source codes for supplemental modules +# SUPPL_LIB_SRCS: list of C source codes for supplemental library +# SUPPL_EXAMPLE_FILES: README and other example files +# SUPPL_EXAMPLE_PROGS: Example scripts +# + +set (SUPPL_NAME geology) +AUX_SOURCE_DIRECTORY (longopt SUPPL_LONG_OPT_H) +set (SUPPL_PROGS_SRCS psstereonet.c ${SUPPL_LONG_OPT_H}) +set (SUPPL_LIB_SRCS ${SUPPL_PROGS_SRCS}) +set (SUPPL_EXAMPLE_FILES README.geology) diff --git a/src/geology/README.geology b/src/geology/README.geology new file mode 100644 index 00000000000..ad50dd16138 --- /dev/null +++ b/src/geology/README.geology @@ -0,0 +1,42 @@ +GEOLOGY Supplemental Package +Distributed under the GNU Lesser Public License; see file +LICENSE.TXT in main GMT directory. + +------------------------------------------------------ +Author: Federico Esteban +Date: August 22, 2026 +Version: 6 (for GMT 6 release). +------------------------------------------------------ + +This directory contains: + + README.geology This document + psstereonet.c psstereonet C program code + longopt/psstereonet_inc.h long-option translation table + +This package contains programs that may be useful for +structural geologists. Currently there is one program: + + stereonet Plot structural geology data on a stereonet + +A stereonet (also called a stereographic projection, a Schmidt net +when equal-area or a Wulff net when equal-angle) is the standard way +to display the orientation of planes and lines measured in the field. +The module converts the strike and dip of planes, or the trend and +plunge of lines, into the longitudes and latitudes expected by the GMT +azimuthal projections, and then draws the cyclographic traces (great +circles) of the planes and/or the points where the poles and lines +pierce the hemisphere. + +Documentation via RST files produce HTML and MAN pages that +give full description of the program interfaces. + + +* Acknowledgement + +The module implements the recipe that Rom1, Joaquim Luis and Federico +Esteban worked out in the GMT forum thread + + https://forum.generic-mapping-tools.org/t/plot-a-schmidt-stereonet/6067 + +which in turn is the basis for the stereonet function in GMT.jl. diff --git a/src/geology/longopt/psstereonet_inc.h b/src/geology/longopt/psstereonet_inc.h new file mode 100644 index 00000000000..8aa69c52063 --- /dev/null +++ b/src/geology/longopt/psstereonet_inc.h @@ -0,0 +1,51 @@ +/*-------------------------------------------------------------------- + * + * Copyright (c) 1991-2026 by the GMT Team (https://www.generic-mapping-tools.org/team.html) + * See LICENSE.TXT file for copying and redistribution conditions. + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation; version 3 or any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * Contact info: www.generic-mapping-tools.org + *--------------------------------------------------------------------*/ + +#ifndef PSSTEREONET_INC_H +#define PSSTEREONET_INC_H + +/* Translation table from long to short module options, directives and modifiers */ + +static struct GMT_KEYWORD_DICTIONARY module_kw[] = { + /* separator, short_option, long_option, + short_directives, long_directives, + short_modifiers, long_modifiers, + transproc_mask */ + { 0, 'A', "azimuth", "", "", "", "", GMT_TP_STANDARD }, + GMT_C_CPT_KW, + { 0, 'G', "fill", + "p,P", "bit,bitreverse", + "b,f,r", "bg|background,fg|foreground,dpi", + GMT_TP_STANDARD }, + { 0, 'L', "symbol_pen", "", "", "", "", GMT_TP_STANDARD }, + { 0, 'M', "dump", + "c,p", "traces,poles", + "", "", + GMT_TP_STANDARD }, + { 0, 'S', "symbol|style", + "-,+,a,A,c,C,d,D,e,g,G,h,H,i,I,j,k,l,n,N,p,r,R,s,S,t,T,w,x,y", + "xdash,plus,star,star_area,circle,circle_area,diamond,diamond_area,ellipse,octagon,octagon_area,hexagon,hexagon_area,invtriangle|inverted_tri,invtriangle_area,rotrectangle|rotated_rec,custom,letter,pentagon,pentagon_area,point,rectangle,roundrectangle|roundrect,square,square_area,triangle,triangle_area,wedge,cross,ydash", + "t,f,j,s,i,a,r,p", "text,font,justify,corners,inner,arc,radial,pen", + GMT_TP_STANDARD }, + { 0, 'T', "type", + "d,l,p", "dipdir|dip_direction,line,plane", + "r,u", "rake,upper", + GMT_TP_STANDARD }, + GMT_W_PEN_KW, + { 0, '\0', "", "", "", "", "", 0 } /* End of list marked with empty option and strings */ +}; +#endif /* !PSSTEREONET_INC_H */ diff --git a/src/geology/psstereonet.c b/src/geology/psstereonet.c new file mode 100644 index 00000000000..6da97d2d6c3 --- /dev/null +++ b/src/geology/psstereonet.c @@ -0,0 +1,765 @@ +/*-------------------------------------------------------------------- + * + * Copyright (c) 1991-2026 by the GMT Team (https://www.generic-mapping-tools.org/team.html) + * See LICENSE.TXT file for copying and redistribution conditions. + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation; version 3 or any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * Contact info: www.generic-mapping-tools.org + *--------------------------------------------------------------------*/ +/* + * Brief synopsis: psstereonet reads structural geology observations given as + * pairs of angles (planes as strike/dip or dip-direction/dip, or lines as + * trend/plunge) and plots them on a stereonet, i.e., the projection of the + * lower (or upper) hemisphere of a unit sphere onto the horizontal plane. + * + * A stereonet is not a new map projection: it is simply a hemisphere seen from + * above, centered on the nadir. Hence we use the standard GMT azimuthal + * projections centered on (0,0), whose default 90-degree horizon gives exactly + * one hemisphere: + * + * -JA0/0/ Lambert azimuthal equal-area -> Schmidt net (equal-area) + * -JS0/0/ Stereographic (equal-angle) -> Wulff net (equal-angle) + * + * With that setup the north pole of the sphere (lat = +90) plots at the top of + * the net, so azimuth is measured clockwise from the top of the plot while the + * angular distance from the center of the net is 90 minus the plunge of the + * line being plotted. A line with trend T and plunge P is therefore the point + * (lon, lat) obtained from the unit vector + * + * x = cos (90-P), y = sin (90-P) * sin (T), z = sin (90-P) * cos (T) + * + * while the cyclographic trace (great circle) of a plane with strike S and dip + * D is the meridian lon = 90-D rotated by S about the x-axis, i.e., about the + * axis that points at the center of the net. A point on that trace is given + * by sweeping a parameter t from -90 to +90; since the standard geological + * rake (or pitch) R of a lineation on the plane runs from 0 at the strike + * azimuth, through 90 at the down-dip direction, to 180 at the opposite end + * of the strike, the two parameters are simply related by t = 90-R. + * + * Author: Federico Esteban + * Date: 22-AUG-2026 + * Version: 6 API + * + * Based on the recipe worked out by Rom1, Joaquim Luis and Federico Esteban in + * https://forum.generic-mapping-tools.org/t/plot-a-schmidt-stereonet/6067 + */ + +#include "gmt_dev.h" +#include "gmt_internals.h" /* For gmtlib_B_is_frame, so we can tell a -B frame setting from axis intervals */ +#include "longopt/psstereonet_inc.h" + +#define THIS_MODULE_CLASSIC_NAME "psstereonet" +#define THIS_MODULE_MODERN_NAME "stereonet" +#define THIS_MODULE_LIB "geology" +#define THIS_MODULE_PURPOSE "Plot structural geology data on a stereonet" +#define THIS_MODULE_KEYS "X},>DM,CC(" +#define THIS_MODULE_NEEDS "JR" +#define THIS_MODULE_OPTIONS "-:>BJKOPRUVXYbdefghipqstxy" + +/* What the two input angles mean and hence what we can draw */ +enum psstereonet_types { + PSSTEREONET_PLANE = 0, /* -Tp: azimuth is the strike (right-hand rule) */ + PSSTEREONET_DIPDIR, /* -Td: azimuth is the dip direction */ + PSSTEREONET_LINE}; /* -Tl: azimuth is the trend and dip is the plunge */ + +/* Which of the derived data sets -M should write */ +enum psstereonet_dumps { + PSSTEREONET_DUMP_AUTO = 0, + PSSTEREONET_DUMP_TRACE, /* -Mc: the cyclographic traces (great circles) */ + PSSTEREONET_DUMP_POINT}; /* -Mp: the poles, or the lines if -Tl */ + +#define PSSTEREONET_N_TRACE 181 /* Points used to draw one great circle (i.e., 1 degree steps) */ +#define PSSTEREONET_DEF_WIDTH 15.0 /* Default width (diameter) of the net, in cm */ +#define PSSTEREONET_DEF_ANNOT 30.0 /* Default azimuth annotation interval */ +#define PSSTEREONET_DEF_TICK 10.0 /* Default azimuth tick interval */ +#define PSSTEREONET_DEF_SYMBOL "c0.15c" /* Default symbol for the poles and lines */ +#define PSSTEREONET_DEF_PEN "default" /* Default pen for traces and symbol outlines */ +#define PSSTEREONET_DEF_FRAME 2 /* Number of default -B options we may add */ + +struct PSSTEREONET_CTRL { + struct PSSTEREONET_Out { /* -> */ + bool active; + char *file; + } Out; + struct PSSTEREONET_A { /* -A[[/]] */ + bool active; + bool draw; /* False if -A0 was given */ + double annot, tick; + } A; + struct PSSTEREONET_C { /* -C */ + bool active; + char *string; /* Since we will simply pass this on to plot */ + } C; + struct PSSTEREONET_G { /* -G */ + bool active; + char *string; /* Since we will simply pass this on to plot */ + } G; + struct PSSTEREONET_L { /* -L */ + bool active; + char *string; /* Since we will simply pass this on to plot */ + } L; + struct PSSTEREONET_M { /* -M[c|p] */ + bool active; + unsigned int mode; + } M; + struct PSSTEREONET_S { /* -S[] */ + bool active; + char *string; /* Since we will simply pass this on to plot */ + } S; + struct PSSTEREONET_T { /* -T[d|l|p][+r][+u] */ + bool active; + bool upper; /* True if +u, i.e., plot on the upper hemisphere */ + bool rake; /* True if +r, i.e., a third column gives the rake of a lineation on the plane */ + unsigned int mode; + } T; + struct PSSTEREONET_W { /* -W */ + bool active; + char *string; /* Since we will simply pass this on to plot */ + } W; + struct PSSTEREONET_l { /* -l