GMT绘图攻略

目录
  1. 坐标轴
    1. 2D
    2. 3D
    3. 坐标单位
  2. 色标与图例
    1. 色标
    2. 图例
    3. 光照
  3. 断层滑动分布三维图
  4. 震源球的绘制
  5. 添加说明文字
  6. 绘制断层
  7. 绘制图中图
  8. 绘制彩色图
    1. 格式转化
    2. 读取/出文件操作
    3. 绘制CPT

GMT可以绘制很多复杂的图形，无论是论文还是报告中都是非常用，但同时它也很复杂，想要学号记住很难。因此，博主开设了这个栏目，一来向大家提供帮助，而来自己也可以加深理解。

坐标轴

2D

1、设置坐标轴带箭头

1
gmt gmtset MAP_FRAME_TYPE = graph

2、设置标注在内部

1
gmt gmtset MAP_FRAME_TYPE = inside

3、调整坐标轴的方向，一般来说GMT绘图是左小右大，下小上大，不过我们如果需要坐标轴反向的话，例如在表示深度时，常常需要向下为大

我们需要在-J中设置，没错，不是-B，这里很多同学可能找错地方了。

我们只需设置

1
2
3
-JX-4i/4i #represent that x axis reversal
#similarly we can set 
-JZ-4i    #represent that z axis reversal

3D

绘制三维坐标时，可以选择

1
2
3
4
5
-BZ #plot the z axis with scale and callout
-Bz #plot the z axis only with scale
-Bz1234 #by default ,z axis will only be plotted one time at left-down
        #so you can select 1234 reprensent plot the other z axises 
        #counterclock

我们在绘制三维图时，常常需要添加Z轴，但默认是绘制左下角的轴，但根据我们一般的观测视角，绘图会使得其他的轴全部反向，这时就需要绘制左上角的轴

1
-BZ4 #这样就可以了

坐标单位

在绘制坐标单位时，使用+l参数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
-Bxa3g1+ldepth
#if label contains ' ',then use "" 
-Bxa3g1+l"wave velocity(km/s)"
#如果单位里含有上下标，使用转义符
@+:开启/关闭上标
@-:开启/关闭下标
@_:开启/关闭下划线
@#:开启关闭小型大写字母
@:<size>: :切换另一字体大小
@%<font>%:切换另一字体
@;<color>;:切换另一字体颜色
@@:打印@本身
例如：
-Bxa3g1+l"density(kg/m@+3@+)"

色标与图例

色标

色标绘制使用psscale

1
gmt psscale -Dx1i/1i+e+ml+w2i/0.2i+o1c/2c -Ccpt.cpt -J -K -S -Ba1+l"bar" >ps.ps
1
2
3
4
5
6
7
8
9
-D的选项有很多，先说附加项

1、上下加一个三角形：在-D中加一个'+e' or 'eb' for bottom only and 'ef' for front only

2、把标注放到另一边：在-D中加一个'+ml'

3、调整色标的走向：+h(horizonal) +v(vertical)

4、在原点处再次偏移：+o

这里主要说一下几点，首先是色标的位置

这里使用-D说明位置，有5种分别是g|j|n|x|J

1
2
3
4
5
6
7
8
9
-Dg:use the position used in data 
#for example -Dg135/24 represent 135E 24N
-Dx:use the position used in basemap
#for example -Dx7.5i/1c 
-Dj:指定修饰物上的锚点，有9种，分别是L(left),R(rigth),T(top),B(bottom),M(middle)，C(center)的组合
#vertical use C horizonal use M
#-R and -J needed
-Dn:使用归一化坐标，X，Y均从0-1
-DJ:与-Dj呈镜像，在需要绘制在底图之外时使用

其次是修饰色标

主要有给色标去除中间黑线-S

图例

图例使用pslegend

主要使用的参数有

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
gmt pslegeng -R -J -Dx+w<width><height> -F <<EOF >ps.ps
#
D 0.1i 1p
H 16 Times-Roman LEGEND
S 0.6i - 1i - thick,blue,dashed/solid/dotted 1.2i Vp
L 14 Times-Roman C Earthquakes > Mw4
EOF
#description
#-F是绘制一个图例框，可以使用-F+g<fill>+p<pen>加上背景和框，一般使用-F+gazure3+p0.2p,black
#-D是为了调整图例的位置，参见上述，其中-Dx是以地图为基准的，而GMT默认地图偏移1i/1i,+w为必须
##表示略去开头是#的行数
#D表示绘制一条线，指明左右与图例框的距离和字号
#H绘制一个图例的标题，需要指明字号，字体，标题
#S绘制一个图形，线段等，参见-S，指明中心与左侧距离，图形，长度/尺寸，填充，pen，标号的偏移，标号
#L绘制一个文字，和H类似，不过可以用LCR表示左对齐，中对齐，右对齐

光照

-I选项是非常常用的

1
2
-I<size>
#size are from -1 to 1

在绘制彩色图时，不加光照就是原色调，加光照，值越大越白

断层滑动分布三维图

在做大地震震源分析时，常常需要给出断层的三维滑动分布，我们可以利用psxyz做到。目前关于psxyz的语法和实例不过，中文社区甚至没有翻译过来。这里博主给大家开一个路

我们的思想是利用psxyz绘制一个个矩形框，根据每个矩形框的滑移值来绘制图像。

我们以USGS给出的20181130地震断层模型文件param.txt作为源文件

首先需要把它变成能够被GMT拾取的文件格式， 即

1
2
3
4
5
6
7
8
9
> -Z<slip>
lons1 lons2 deps1
lons2 lons2 deps2
lone1 lone1 depe1
lone2 lone2 depe2
> -Z<slip>
...
...
..

程序可以自己编，不是很难，注意-Z和slip之间不能有空格，这里用一个trim和adjustl函数就可以解决了。如果还有问题，可以联系博主获取支持。

这样我们就可以键入

1
2
#!/bin/bash                                                                 R=-150.6/-149.3/61/62/0/60                                                           txt=plot2018                
gmt psxyz $txt.txt -R$R -JX6i/4i -Jz-0.2 -BSWZ2 -Bxa0.3+lLontitude -Bya0.2+lLatitude -Bza20+lDepth   -p210/35 -K -I0.7 -L -Cpsxyz.cpt -W0.5p,black >psxyz.ps

来得到断层的滑移分布图了

震源球的绘制

基本格式是

1
2
gmt psmeca -R -J -Sa$size -C1.2p -Gcyan -K <<EOF >>$nameps.ps
lon lat depth strike dip rake Mw newX newY label

注意绘制震源球的方式有很多中，比较常用的是-Sa和-Sz两种

-Sa：指定3参数和震级即可，绘制出来的图，压缩和拉伸边界有明显的黑线

-Sz：指定地震矩张量的分量，一般是gCMT提供的数据来源是这样的，绘制出来的图没有边界黑线

1
2
3
-C:指定震源位置和图中位置的连线参数
-G:指定压缩区域的颜色
-F:指定压缩轴和伸张轴的属性

添加说明文字

基本语句是：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
gmt pstext -R -J -K -F+f22p,7,blue+a70 <<EOF >>$nameps.ps
x y label
EOF
-F:指定了文字的很多信息，
	+f:指定文字的大小，字体，颜色
	+a:指定文字的角度（以逆时针旋转为正）
#个人使用来看，基本是
	海洋水域的标识使用
		7号字体，blue
	陆上构造标识（如半岛，盆地）使用
		28号字体，black
	城市名使用
		31号字体，black	
	断层名使用
		1号字体，12p，black
以上仅供参考

绘制断层

一些不知道参数的断层可以去官网查询，得到kmz或者shp文件在GoogleEarth中打开，逐个描点得到坐标。

基本语法是利用psxy

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
#绘制简单的线
gmt psxy -R -J -W<pen> -K <<EOF >>$nameps.ps
x y
x y
x y
EOF
#绘制带符号的断层
gmt psxy -R -J -W<pen> -Gblack -Sf2c/0.07i+t+l -K <<EOF >>$nameps.ps
-Sf是专门绘制断层的参数，2c表示每2cm回一个标记0.07i表示标记大小
	+t:表示绘制的标记是三角形
	+l:表示绘制的标记在左侧/右侧是+r
	+b:是矩形
	+c:是圆形
	+s:绘制左旋和右旋
		注意2c可以换成负数，则表示绘制图形的个数，不过此时无法使用-o选项了
	+o<size>:绘制时让初始符号偏移距离
-G:符号填充色

绘制图中图

在一些论文中经常需要给某一个地方的大范围说明图，我们绘制时可以

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
gmt psbasemap -R$Rbig -J -Bse -K -X<pos> -Y<pos> >>$nameps.ps
gmt pscoast   -R$Rbig -J -Bse -K -O -Dl -W1/0.5p -N1 -Ggray >>$nameps.ps
gmt psxy      -R$Rbig -J -W<pen> -L -O <<EOF >>$nameps.ps
x y
...
EOF
#先绘制底图，注意根据图的位置选择
	-Bse或是其他的轴的特性，也可以选择是否绘制标尺
#绘制海岸线时
	-D选项一般选的比较低，一共5级
		f<full>   h<high>  i<intellectual> l<low> c
	-N有四级，1就是国界了，2州界。。。
	-W绘制水域相关，先设置等级，再设画笔，否则默认画笔
	-G选择陆地颜色
#绘制大图选框
	基本就是利用psxy，注意选择
	-L使得图形围起来成一个矩形框

绘制彩色图

格式转化

我们取一个最一般的情况，我们有一个txt文件（或者dat），里面的信息是包含lon,lat,value,….分别是第1～3列，当然还包含其他列的信息。

我们首先需要转化为grd文件或者nc文件。

1
2
3
4
5
gmt surface filename.txt -R -T0.01 -I0.1/0.1 -i0,1,2 -C0.001 -Gfilename.grd
#-I表示插值密度
#-C表示收敛值
#-T表示张力系数，在0~1之间，控制曲率曲面
#-i是通用选取文件的参数，表示选取文件中的文件列数，从0起算，同时可以做一些运算

读取/出文件操作

主要是-i和-o的操作

1
2
3
4
5
6
-i0,1,2
#表示取文件的1,2,3列
-i0-2
#与上式相同
-i2+s2+o10+l
#表示取第3列，并对它先乘2，再加10，最后取log10

对于-o，与-i相同，不过没有运算

另外，一些小符号可以是的输入的x，y交换，即相当于-i1,0

1
-:

还有-h可以略去开头多少行

1
2
gmt ... -h3
#表示前三行不读取

-s则可以选择数据中还有NaN的处理方法

1
2
3
4
5
6
7
8
gmt select input.dat -s
#若第三列还有NaN，则不输出
.... -sa
#只要有NaN，就不输出
.... -sr
#只输出含有NaN的行
.... -s2,3
#抑制3,4列的NaN数据

绘制CPT

使用

1
2
3
4
5
6
7
gmt makecpt -Cblue,cyan,green,yellow,red -T-1/1-0.05 -Z -D >cpt.cpt
#-C表示色域
#-T表示值域
#-Z表示颜色连续
#-D表示越界颜色取最大值颜色
#其中-C有很多内置的色彩搭配，默认是rainbow，具体可以参考
#http://soliton.vm.bytemark.co.uk/pub/cpt-city/