GSEASRX

For Earthquake Science

0%

在psxy中使用cpt填色

发表于 更新于 分类于 阅读次数： Disqus：

  绘制震中分布图时经常需要绘制多年地震，并需要将不同年份地震进行区分，而区分的方式通常是给其填充不同的颜色。常规的绘图方法中将每一年的地震目录用awk提出后再绘制，此方对于每年的地震都需要单独一条命令进行绘制，并给其一个单独的颜色。如此一来同样的绘图语句就要重复写好几遍，个人觉得还是有点不方便。
  好在gmt的—C<cpt>模块为我们提供了使用cpt文件给图中元素进行填色的功能，这样的话，只需要做好数据表便可一条命令绘制所有地震，并且使用cpt进行渲染，做到按照发生年份进行区分，同时方便ColorBar形式的时间图例的绘制。

psxy模块下-C选项简介

在psxy模块中—C<cpt>选项用于指定CPT文件或颜色列表。该选项后跟一个CPT 文件名，也可以使用-C<color1>,<color2>,... 语法在
命令行上临时构建一个颜色列表，其中 对应Z 值为0 的颜色，
对应Z 值为1 的颜色，依次类推。

  1. 若绘制符号（即使用-S选项），则符号的填充色由数据的第三列 Z 值决定，其
    他数据列依次后移一列
  2. 若绘制线段或多边形（即未使用-S选项），则需要在多段数据的头段中指定
    -Z<val>，然后从cpt文件中查找<val>所对应的颜色，以控制线段或多边
    形的线条颜色

下面的例子展示了-C<color1>,<color2>..用法:
绘制两条不同颜色的线段

1
2
3
4
5
6
7
8
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

图1

举个栗子

上例为GMT6-documentation中的例子

下面以绘制2008年以来全球7级以上地震震中分布图为例

bash下绘图代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
#!/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

结果:

图2

地震目录格式如下：

1
2
3
4
5
6
7
8
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绘图需要。

  1. 本例子利用for循环加awk实现文本内容的提取和格式的设定，即

     经度 纬度 Z值(年份) 震级(大小)

  2. makecpt模块将文件rainbow.cpt的Z值设置为2008～2019，并且划分间隔为1，生成新的cpt文件；

  3. 因此绘图时使用-C选项，会直接读取第三列为填充颜色；并可以绘制colorbar作为图例，省去了繁琐的图例绘制。

绘制颜色变化的曲线

想要绘制一条颜色变化的线段，下面是演示代码：

此例引自Seisman播客

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
#!/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

结果：

图3

👉简单解释一下：

  1. makecpt 命令制作了一个 - 2 到 2 间隔为 1 的 CPT 文件
  2. psxy 命令中使用了 -C 选项，此时需要输入数据是多段表数据，且每段数据的段头记录中，需要有 -Z<val> 以指定每段数据的 Z 值
  3. 实际绘图时，对于每段数据，命令会读取该数据数据的段头记录中的 -Z<val> 中的 Z(val)值，然后到 CPT 文件中查找 Z 值所对应的颜色，作为该段线段的颜色。