I have a 2-column file containing azimuth and plunge values, representing the poles to fault planes. I would like to plot them on a Schmidt stereonet (equal-area projection), but I can’t get it to work — the points appear outside the stereonet circle.
Here is my current script:
gmt begin "STEREO_${type}" png
gmt basemap -R0/360/0/90 -JA0/0/10 -Ba20g10
gmt plot STEREO_POLE_${type}.txt -Sc0.1 -Gred
gmt end
With the first line of STEREO_POLE_${type}.txt being: 124 25
I think I could help you but I don’t remind exactly how the Schmidt stereonet works. Could you show me what is the expected output and what is the meaning of the input values.
I think this is your answer.
First, I use a -JA command to plot only the stereonet gridlines.
Then, in the second line, I switch to the -JP option to plot the polar data. For this, the azimuth is used directly, but the plunge needs to be converted to its complementary angle (i.e., 90° – plunge), which I do using awk. Please let me know if my script works for you.
#Create data file
cat > data << EOF
124 25
EOF
width=10c # Set width of the stereonet
gmt begin stereonet png
gmt basemap -Rd -JA0/0/$width -Bg15 # Plot stereonet grid
awk '{print $1,90-$2}' data | gmt plot -JP$width+a -Sc0.3c -Gred
gmt end
It works ! I haven’t think of -JP !
Thank you !
By the way, looking at the doc for -JP you can directly use -JP$width+a+fe to flip the radial direction to point inwards and it avoid the use of awk
Sorry but after few try I don’t thinks it works
Imagine a plane with a strike of 180 and a dip of 60, the pole has a trend of 90 and a plunge of 30.
So we use :
No, the Schmidt chart is Lambert with center at 0,0 (see CookBook), but it is a Lambert proj nevertheless and expects coordinates as lon,lat. Azimuth and Dip should be converted to the lon,lat that it expects.
Using -JA0/-90 might do the trick for points but I’m not sure if it’s possible to use it to plot planes
I finally managed to plot a stereonet with both fault planes and poles. It was quite tricky, and I couldn’t find clear answers anywhere (I guess I might be the only one stubborn enough to insist on doing this with GMT!).
As Joaquim pointed out, it’s necessary to convert strike and dip values to longitude and latitude. I wasn’t familiar with this kind of transformation, so it took me a while to figure it out. I don’t know if there’s a more straightforward method, but here’s how I did it:
Assume the dip is the longitude and set latitude to 0.
Convert these spherical coordinates (lat/lon) to 3D Cartesian coordinates (x/y/z).
Apply a 3D rotation matrix to rotate the points according to the strike. (Be very careful with angle conventions e.g., strike is measured clockwise from north, while Cartesian angles are typically counterclockwise from the x-axis. And angle units degree et radian)
Convert the rotated Cartesian coordinates back to spherical (lat/lon).
If there’s a simpler or more elegant way to do this, feel free to share!
If any GMT developer happens to read this a gmt stereonet command in a futur version would totally make the day of every structural geologist out there (myself included)!
Attached below is my script and the resulting plot.
gmt begin "STEREO" png
gmt basemap -R0/360/-90/90 -JA0/0/10 -Bpg2 -Bsg20 --MAP_GRID_PEN_PRIMARY=0,gray --MAP_GRID_PEN_SECONDARY=0,black
# Create the points along the fault plane and plot it
awk '
function asin(x) { return atan2(x, sqrt(1-x*x)) }
function acos(x) { return atan2(sqrt(1-x*x), x) }
function atan(x) { return atan2(x,1) }
{
pi = atan2(0,-1)
dip = -(90-$3) * (pi/180)
strike = (180+$2) * (pi/180)
print ">"
i = 0
for (val = -90; val <= 90; val++) {
lat[i] = val * (pi/180)
lon[i] = dip
i++
}
for (j = 0; j < i; j++) {
x = cos(lat[j]) * cos(lon[j])
y = cos(lat[j]) * sin(lon[j])
z = sin(lat[j])
yrot = y * cos(strike) + z * sin(strike)
zrot = -y * sin(strike) + z * cos(strike)
r = sqrt( x**2 + yrot**2 + zrot**2)
latrot = asin( zrot / r ) * (180/pi)
lonrot = atan2( yrot , x ) * (180/pi)
print lonrot, latrot
}
delete lat
delete lon
}' data | gmt plot -W0.5,red -l"fault plane"
# Calculation of the pole and plot
awk '
function asin(x) { return atan2(x, sqrt(1-x*x)) }
function acos(x) { return atan2(sqrt(1-x*x), x) }
function atan(x) { return atan2(x,1) }
{
pi = atan2(0,-1)
dip = (-$3) * (pi/180)
strike = ($2) * (pi/180)
x = cos(dip)
y = sin(dip)
z = 0
yrot = y * cos(strike) + z * sin(strike)
zrot = -y * sin(strike) + z * cos(strike)
r = sqrt( x**2 + yrot**2 + zrot**2)
lat = asin( zrot / r ) * (180/pi)
lon = atan2( yrot , x ) * (180/pi)
print lon, lat
}' data | gmt plot -Sx0.1 -l"pole"
gmt basemap -R0/360/-90/90 -JP10+a -Ba15
gmt legend -DjBR+w3+jBR -F+gwhite+p0,black
gmt end
Hey guys!
Hope you’re all doing well!
Just wondering if there’s been any progress on the attempts to create a GMT stereonet function directly in pure GMT?
