Plot a schmidt stereonet

Q&A
11

Hi everyone,

Thanks a lot for your help!

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:

  1. Assume the dip is the longitude and set latitude to 0.
  2. Convert these spherical coordinates (lat/lon) to 3D Cartesian coordinates (x/y/z).
  3. 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)
  4. 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)! :joy:

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

STEREO_DATA
STEREO_DATA1288×1289 191 KB