引言
绘制震中分布图是经常需要绘制多年地震，并需要将不同年份地震进行区分，而区分的方式往往是给起填充不同的颜色。常规的绘图方法中将每一年的地震目录用awk提出后再绘制，此方法绘图时，每年的地震需要单独一条命令进行绘制，并给其一个单独的颜色。
But—C<cpt>模块为我们提供了使用CPT文件给图中元素进行填色的功能，如此一来，只要做好数据表便可一条命令进行绘图，并以ColorBar作为图例。
psxy模块下-C选项简介
在psxy模块中—C<cpt>选项用于指定CPT文件或颜色列表。该选项后跟一个CPT 文件名，也可以使用-C<color1>,<color2>,...语法在命令行上临时构建一个颜色列表，其中<color1> 对应Z 值为0 的颜色，<color2>
对应Z 值为1 的颜色，依次类推。
- 若绘制符号（即使用
-S选项），则符号的填充色由数据的第三列 Z 值决定，其
他数据列依次后移一列 - 若绘制线段或多边形（即未使用
-S选项），则需要在多段数据的头段中指定
-Z<val>，然后从cpt文件中查找<val>所对应的颜色，以控制线段或多边
形的线条颜色
下面的例子展示了-C<color1>,<color2>..用法:
绘制两条不同颜色的线段
gmt psxy -JX10c/10c -R0/10/0/10 -B1 -Cblue,red -W2p > test.ps << EOF
> -Z0
1 1
2 2
> -Z1
3 3
4 4
EOF
image
举个栗子
上例为GMT6-documentation中的例子
以下为自己的栗子：
绘制2008年以来全球7级以上地震震中分布图
bash下绘图代码如下
#!/bin/bash
# GMT v6.0.0
# 提取每一列数据
awk '{if (substr($0,2,4)>=2008) print substr($0,1,15),substr($0,16,6),substr($0,22,7),substr($0,29,4),substr($0,33,3),substr($0,36,3)}' Ms7.EQT > 2018ms7
# 将数据逐年提取为：“经度 纬度 年度 震级”格式
# 下面为for循环下的实现
for i = ((i=2008;i<=2019;i++));
do
awk '{if (substr($1,1,4)=="'$i'") print $3,$2,"'$i'",$4}' 2018ms7 >> eqdis
done
gmt set FONT 11p,4
gmt set MAP_FRAME_TYPE plain
gmt set MAP_FRAME_PEN 1p,black
# 创建cpt文件
gmt makecpt -Crainbow -T2008/2019/1 > cpt.cpt
# 绘图
gmt pscoast -JN6i -R-30/330/-90/90 -Bx60f30 -By30 -Dc -A10000 -Gwhite -Sskyblue -K > eq7dist.ps
awk '{print $1,$2,$3,$4*0.04}' eqdis |gmt psxy -J -R -Sc -W.1p,white -Ccpt.cpt -K -O >> eq7dist.ps
# 绘制colorbar
gmt psscale -Ccpt.cpt -Dx0/0+w5.6i/.2c+jBL+h+o.3c/-1c -Bx1 -By+l"Years" -K -O >> eq7dist.ps
gmt psxy -J -R -T -O >> eq7dist.ps
# 转换至位图
gmt psconvert eq7dist.ps -A -Tj -E300 -P
rm *.cpt gmt.* eqdis
结果如下
download.png
地震目录格式如下：
# txt文本
20181211102630-58.32 -26.387.00150000 南桑威奇群岛地区
20181201012927 61.33-150.057.30 40000 美国阿拉斯加
20181205121806-21.85 169.407.50 10000 洛亚蒂群岛
20181211102630-58.32 -26.387.00150000 南桑威奇群岛地区
20190222181720 -2.17 -76.877.50140000 厄瓜多尔
20190301165041-14.55 -70.107.00260000 秘鲁
20190507051936 -6.95 146.507.10130000 巴布亚新几内亚
20190514205828 -4.17 152.487.60 30000 新不列颠岛地区
因此在绘图之前需要对文本格式进行一定的处理，以满足gmt绘图需要。
- 本例子利用for循环加awk实现文本内容的提取和格式的设定，即
经度 纬度 Z值(年份) 震级(大小)
-
makecpt模块将cpt文件rainbow的Z值设置为2008～2019，并且划分间隔为1，生成新的cpt文件； - 因此绘图时使用-C选项，会直接读取第三列为填充颜色；并可以用一条命令绘制colorbar来作为图例，省去了原来繁琐的图例绘制。
绘制颜色变化的曲线
想要绘制一条颜色变化的线段，下面是演示代码：
此例引自Seisman播客
#!/bin/bash
# GMT v5.2.1
gmt makecpt -Crainbow -T-2/2/1 > lines.cpt
gmt psxy -JX15c/4c -R0/6/0/4 -B1 -Clines.cpt -W2p > test.ps << EOF
> -Z-1.5
1 2
2 2
> -Z-0.5
2 2
3 2
> -Z0.5
3 2
4 2
> -Z1.5
4 2
5 2
EOF
绘图效果如下图：
image
简单解释一下：
-
makecpt命令制作了一个 - 2 到 2 间隔为 1 的 CPT 文件 -
psxy命令中使用了-C选项，此时需要输入数据是多段表数据，且每段数据的段头记录中，需要有-Z<val>以指定每段数据的 Z 值 - 实际绘图时，对于每段数据，命令会读取该数据数据的段头记录中的
-Z<val>中的 Z 值，然后到 CPT 文件中查找 Z 值所对应的颜色，作为该段线段的颜色
