「GMT」How to plot fault geometry with coseismic displacement by GMT

Posted on | In | | Visitors:
The very important thing about our researches are fault slip and GPS data. This webpage is aiming at showing how to plot these two things by GMT.

Click here to download example files

Mainly commands:

  • psxy - to plot fault geometry with fault slip
  • pscoast - to plot basemap
  • psvelo - to plot GPS velocity

Material needed：

  1. Fault_geo1.dat: Fault geometry 1 with coseismic slip
  2. Fault_geo2.dat: Fault geometry 2 with coseismic slip
  3. GPS_coseis.dat: GPS data corresponding to the coseismic slip.
  4. Plot_fault.ksh: Plotting script.

Material Format

The format for Fault_geo*.dat is

1
2
3
4
5
6
-Z ${slip magnitute}
${lon1} ${lat1}
${lon2} ${lat2}
${lon3} ${lat3}
${lon4} ${lat4}
...

The format for GPS_coseis.dat is

1
${lon}  ${lat} ${east_component} ${north_component} ${east_significant} ${noth_significant}

Plotting script

Next we’ll mainly focus on our plotting script.

1. Give name to some initial parameters

1
2
3
4
5
6
7
8
#!/bin/bash
ps=plot.ps
fault_geo1=fault_geo1.dat
fault_geo2=fault_geo2.dat
gps_name=GPS_coseis.dat
cpt=fault_slip.cpt
R=-157/-144/55/62
J=M5i

2. Plot basic picture

1
2
gmt makecpt -Chot.cpt -G0.4/1 -Ic -T0/22 > $cpt
gmt pscoast -R$R -J$J -Bx4 -By2 -Da -G229/195/148 -S156/233/253 -K -P > $ps

Using makecpt to define the colormap: Choose hot colormaphen, cut the range to 0.41 with -G, invert colormap with -Ic, choose 022m’s slip range with -G, finally output to $cpt.(see makecpt)
Using pscoast to plot the land and ocean with different color.(see pscoast)

3. Plot fault geometry

1
2
3
4
5
6
gmt psxy $fault_geo1 -R$R -J$J  -W0.1p -L -C$cpt -O -K >> $ps
gmt psxy $fault_geo2 -R$R -J$J  -W0.1p -L -C$cpt -O -K >> $ps
gmt psscale -D3.5i/0.5i+w2.5i/0.25i+h+jCT+m -C$cpt \
        -Bx5+l"Coseismic slip (m)" -O -K -P >> $ps
gmt pscoast -R$R -J$J -Bx4 -By2 -Da -W229/195/148 -O -K -P >> $ps

Here we use psxy to plot the fault geometry, giving different color to fault patches according to slip maganitude and assigned colormap.(see psxy)
And we plot the coast line again to avoid everything being covered by fault geometry.

4. Plot GPS vectors

1
2
3
4
5
6
7
8
    awk '{print $1,$2,$3,$4,$5,$6,0.0}' ${gps_name} |\
    gmt psvelo    -R${R}  -W0.25p,black  -Gblack  -L  -Se0.3c/0.95/25  -J${J} -A0.3c+p1.5p+e -P -V -O -K >> $ps

    gmt psvelo <<END  -h2  -R${R}  -W0.25p,black  -Gblack  -L  -Se0.3c/0.95/12 -J${J} -A0.3c+p1.5p+e -P -V -O  >> $ps
         Long.   Lat.   Evel   Nvel   Esig   Nsig  CorEN SITE
          (deg)  (deg)    (m)        (m)
      -155.5      61.5     2.0    0      0     0     0     2m
END

We use psvelo to plot GPS displacement with the first two lines.(see psvelo)
The left 5 lines are used to plot displacement’s legend.

5. Convert to png

1
2
gmt psconvert $ps -A -Tg
rm $ps gmt.* $cpt

Finally, we can convert the .ps file to .png file for better Visualization and remove files we don’t need anymore.

The output picture should looks like this.

And this is just a simple example for a quick start, you can change and add whatever you like once you get used to the GMT software. Good luck!