在地图上绘制波形

将 SAC 波形数据绘制在地图上也算是常见的需求之一。最常见的情况是,在地图上有一些台站,需要将波形画在台站的附近。

此问题的难点在于,在画地图时,-R 指定的是区域范围,即横轴是经度,纵轴是纬度,而 pssac 在绘制波形时,横轴是时间,纵轴是振幅。因而想要将波形放在特定的位置,不是一件容易的事情。

生成示例

用 SAC 自带的数据,生成几个示例数据,供接下来使用:

SAC> dg sub tele ntkl.z nykl.z onkl.z sdkl.z
SAC> w ntkl.z nykl.z onkl.z sdkl.z

saclst 查看一下这几个波形数据的相关信息:

$ saclst b e stlo stla f *.z
ntkl.z      199.965     1600.95    -114.592     62.4797
nykl.z      199.962     1600.97      -74.53     44.5483
onkl.z      199.618     1600.62    -93.7022     50.8589
sdkl.z      199.098     1600.11    -104.036     44.1204

先绘制地图

在考虑如何把波形画在地图上之前,先得有一个地图。用如下代码绘制一个地图:

1
2
3
4
5
6
7
8
9
10
11
12
13
#!/bin/bash
J=M15c
R=-120/-60/30/65
PS=map.ps
psxy -J$J -R$R -T -K > $PS
# 绘制海岸线
pscoast -J$J -R$R -B10/10 -Ggray -K -O -A1000 >> $PS
# 绘制台站位置
saclst stlo stla f *.z | awk '{print $2, $3}' | psxy -J$J -R$R -Sa0.5c -Gblack -K -O >> $PS
# 绘制台站名
saclst stlo stla f *.z | awk '{print $2, $3,"15 0 0 TR", $1}' | pstext -J$J -R$R -D-0.1c/-0.1c -K -O >> $PS
psxy -J$J -R$R -T -O >> $PS

代码中使用了 saclstawk 提取了每个波形所对应的台站经纬度,并通过管道传递给 psxypstext 命令。

地图效果如下:

第一个波形振幅图

ntkl.z 的波形为例,先不去想如何把波形放在台站附近,只考虑如何将波形画在地图上。

先把代码写出来,如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
#!/bin/bash
J=M15c
R=-120/-60/30/65
PS=map.ps
psxy -J$J -R$R -T -K > $PS
pscoast -J$J -R$R -B10/10 -Ggray -K -O -A1000 >> $PS
saclst stlo stla f *.z | awk '{print $2, $3}' | psxy -J$J -R$R -Sa0.5c -Gblack -K -O >> $PS
saclst stlo stla f *.z | awk '{print $2, $3,"15 0 0 TR", $1}' | pstext -J$J -R$R -D-0.1c/-0.1c -K -O >> $PS
# 绘制波形
pssac -JX4c/2.6c -R200/1600/0/2 -B500/1 -M1 -Ent -K -O ntkl.z >> $PS
psxy -J$J -R$R -T -O >> $PS

生成的效果如下:

解释一些 pssac 中各个选项的含义及作用:

  • -R 选项中 X 轴的范围 200/1600 是波形的时间跨度,这里要绘制整个波形,所以时间跨度由波形的 B 值和 E 值决定;当然也可以只绘制波形中的一小段,这个放后面再说;
  • -JX 中的底图长度 4c 用于控制波形在地图上的长度;
  • -M1 会对波形做振幅归一化,使得波形的最大振幅差在地图上显示的高度为 1** 英寸 **;修改这个值可以控制波形的高度;
  • -R 选择中 Y 轴的范围设置为 0/2 以及 -Ent 中的 n 是固定的,直接用就好,不需要修改;这样设计的原因会在后面解释;
  • -JX 中的 Y 轴长度 2.6c 为底图的高度;由于波形的最大振幅差已经设定为 1 英寸,即 2.54 厘米,底图的高度必须大于等于波形的高度,如果底图高度取 2.54c,可能会导致某些波形被截断,所以要底图高度要稍大一些,可以一点一点增加底图高度看具体效果;
  • 此代码中 -B500/1 为波形加了边框,这样做是为了看起来更直观,稍后会把这个边框去掉;

将波形放对位置

前面已经把波形画到了地图上,还需要将波形移动到台站所在的位置。最直观的想法是在用 pssac 绘制波形时,同时使用 -X-Y 选项移动波形的坐标原点。那么 -X-Y 的偏移量分别是多少呢?这需要借助于 GMT 的 mapproject 命令。

ntkl.z 台站的位置是 (-114.592, 62.4797)mapproject 命令可以计算出任意一点相对于当前底图左下角坐标原点的偏移量:

$ echo -114.592 62.4797 | mapproject -JM15c -R-120/-60/30/65
1.352   12.2477921602

这里的输出结果表明, (-114.6, 62.5) 这一坐标相对于底图左下角原点的距离是 (1.35c, 12.25c)

知道偏移距离之后,给代码加上 -Xa1.35c -Ya12.25c 就可以把波形偏移到台站处了。 -X-Y 选项中在偏移量之前加上了 a,使得将波形偏移到台站坐标处的同时不修改原底图的坐标原点,这样不会影响到后面波形的绘制。

1
2
3
4
5
6
7
8
9
10
11
12
13
#!/bin/bash
J=M15c
R=-120/-60/30/65
PS=map.ps
psxy -J$J -R$R -T -K > $PS
pscoast -J$J -R$R -B10/10 -Ggray -K -O -A1000 >> $PS
saclst stlo stla f *.z | awk '{print $2, $3}' | psxy -J$J -R$R -Sa0.5c -Gblack -K -O >> $PS
saclst stlo stla f *.z | awk '{print $2, $3,"15 0 0 TR", $1}' | pstext -J$J -R$R -D-0.1c/-0.1c -K -O >> $PS
pssac -JX4c/2.6c -R200/1600/0/2 -B500/1 -M1 -Ent -K -O -Xa1.35c -Ya12.25c ntkl.z >> $PS
psxy -J$J -R$R -T -O >> $PS

绘图效果如下:

上图还是有问题。因为 -X-Y 控制的是 pssac 的底图原点相当于地图底图的原点的偏移量,因而导致波形的底图边框左下角恰好位于台站位置处。

从美观等角度考虑,可能希望波形起点位于五角星的正右边一点。因而需要将 mapproject 计算出来的 X 偏移量调大一点,Y 偏移量则应减去底图高度的一半。这个例子中,可以将 -Xa1.35c -Ya12.25c 改成 -Xa1.5c -Ya11c。如果希望波形起点位于五角星的正左边一点,也是可以的,自己指定一个额外的偏移量即可。

绘图效果如下:

绘制所有波形

前面已经把一个台站的波形放对了位置,接下来就需要把所有波形都放对位置。

可以使用 mapproject 命令获取所有台站位置相对于底图左下角的偏移量,如下所示:

$ saclst stlo stla f *.z | awk '{print $2, $3, $1}' | mapproject -JM15c -R-120/-60/30/65
1.352   12.2477921602   ntkl.z
11.3675 4.57807021226   nykl.z
6.57445 6.91932494984   onkl.z
3.991   4.42902501583   sdkl.z

然后在 mapproject 的输出的基础上加上额外的偏移量,波形就可以放在台站正右方了。

代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#!/bin/bash
J=M15c
R=-120/-60/30/65
PS=map.ps
psxy -J$J -R$R -T -K > $PS
pscoast -J$J -R$R -B10/10 -Ggray -K -O -A1000 >> $PS
saclst stlo stla f *.z | awk '{print $2, $3}' | psxy -J$J -R$R -Sa0.5c -Gblack -K -O >> $PS
saclst stlo stla f *.z | awk '{print $2, $3,"15 0 0 TR", $1}' | pstext -J$J -R$R -D-0.1c/-0.1c -K -O >> $PS
pssac -JX4c/2.6c -R200/1600/0/2 -M1 -Ent -K -O -Xa1.5c -Ya11c ntkl.z >> $PS
pssac -JX4c/2.6c -R200/1600/0/2 -M1 -Ent -K -O -Xa11.5c -Ya3.25c nykl.z >> $PS
pssac -JX4c/2.6c -R200/1600/0/2 -M1 -Ent -K -O -Xa6.67c -Ya5.6c onkl.z >> $PS
pssac -JX4c/2.6c -R200/1600/0/2 -M1 -Ent -K -O -Xa4.2c -Ya3.1c sdkl.z >> $PS
psxy -J$J -R$R -T -O >> $PS

绘图效果如下:

上图还是有些问题,比如 nykl.z 台站的波形右边有点过界了, sdkl.z 似乎被限幅了,可以修改 -JX4c/2.6c 来微调。微调的事情就不再说了。

自动化

如果只有几个波形的话,直接为每个波形计算偏移距,然后手写 pssac 语句倒也不麻烦。如果台站很多的话,就需要将整个流程自动化了。一方面是因为要画多个波形的话,手写 pssac 工作量有些大;另一方面,命令中的某些参数存在关联,因而微调某个参数的同时,可能也要同时调整其他参数(比如,上面的代码中,若要微调一下地图的 -J-R ,则所有的偏移量都需要重新计算并修改)。

下面给出用 Perl 写的自动化脚本。之所以用 Perl 而不是用 Bash,一方面是因为自动化过程中涉及到一些浮点数运算,而 Bash 的浮点运算比较麻烦,另一方面是因为我不会写稍复杂的 Bash 脚本。

注意:本脚本只能使用 GMT4 的绘图命令,使用 GMT5 的命令会出现问题,请勿使用。

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
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
#!/usr/bin/env perl
use strict;
use warnings;
my $J = "M15c";
my $R = "-120/-60/30/65";
my $PS = "map.ps";
system "psxy -J$J -R$R -T -K> $PS";
# 绘制海岸线
system "pscoast -J$J -R$R -B10/10 -Ggray -K -O -A1000>> $PS";
my @sacfiles = glob "*.z";
# 绘制台站位置
open(PSXY,"| psxy -J$J -R$R -Sa0.5c -Gblack -K -O>> $PS");
foreach (@sacfiles) {
my ($fname, $stlo, $stla) = split " ", `saclst stlo stla f $_`;
print PSXY "$stlo $stla\n";
}
close(PSXY);
# 绘制台站名
open(PSTEXT,"| pstext -J$J -R$R -D-0.1c/-0.1c -K -O>> $PS");
foreach (@sacfiles) {
my ($fname, $stlo, $stla) = split " ", `saclst stlo stla f $_`;
print PSTEXT "$stlo $stla 15 0 0 TR $fname\n";
}
close(PSTEXT);
# 波形相关参数
my $width = 4; # 波形的底图宽度
my $height = 2.6; # 波形的底图高度
my $M = "-M1"; # 波形在图上的高度为 1 英寸
my $E = "-Ent"; # 波形剖面图类型
my $Rsac = "-R200/1600/0/2"; # 波形的范围
my $Jsac = "-JX${width}c/${height}c";
my $dx = 0.15; # 对 X 偏移量的校正
my $dy = -$height/2.0; # 对 Y 偏移量的校正
foreach (@sacfiles) {
my ($fname, $stlo, $stla) = split " ", `saclst stlo stla f $_`;
# 用 mapproject 计算偏移量
my ($Xshift, $Yshift) = split " ", `echo $stlo $stla | mapproject -J$J -R$R`;
# 对计算出的偏移量做微调
$Xshift += $dx;
$Yshift += $dy;
# 绘制波形
system "pssac $Jsac $Rsac $M $E -Xa${Xshift}c -Ya${Yshift}c -K -O $fname >> $PS";
}
system "psxy -J$J -R$R -T -O>> $PS";

这个 Perl 脚本的自动化程度比较高,对参数的微调基本不会影响到图的整体效果。下面稍微解释一下其中可以微调的变量:

  • $J :地图的投影方式;
  • $R :地图的区域范围;
  • $M :波形在图上的高度(inch);
  • $width :波形底图的宽度(cm),相当于波形在图上的宽度;
  • $height :波形底图的高度(cm),要求比 $M 给出的值稍大;
  • $Rsac :波形的范围,横轴范围需要自己根据数据调整,纵轴范围固定为 0/2
  • $E :固定为 -Ent,如果想要只绘制波形的一部分,可以稍做修改,参考后面一节;
  • $Jsac :波形的投影方式,不需要修改;
  • $dx :对 mapproject 计算出的 X 偏移量的微调;
  • $dy :对 mapproject 计算出的 Y 偏移量的微调;

注:该脚本中所有波形文件共用同一个 $Rsac,这是一个限制。如果波形的 B 和 E 值不同,则需要对每个波形使用不同的 $Rsac,这需要对脚本做一点微调。

只绘制一段波形

前面的代码中绘制的是整条波形,有时候会需要只绘制波形中的一段。比如 P 波的震相到时位于 SAC 文件的头段变量 T0 中,想要将数据中 T0 前 30 秒到 T0 后 30 秒的波形画在地图上。

对前面给出的 Perl 代码稍作修改即可:

  • $E 的值改成 -Ent0
  • $R 的值改成 -R-30/30/0/2 -C

最终得到的效果图如下(示例数据中本没有标记 T0,这里我随便标记了几个 T0 作为 P 波到时):

几点说明:

  • -Ent0n 表示绘制 Y 轴为文件序号的剖面图,因为每个 pssac 命令中只绘制一个波形,因而这个波形会被画在 Y=1 处,因而此处 -R 中 Y 轴的范围取的是 0 到 2。实际上只要 Y 轴的范围包含 1 即可,比如 Y 轴范围取 0.5/1.50/3 均可,可以尝试取不同的值查看其效果;
  • -Ent0t0 表示以 T0 为参考时间,此时 T0 对应的时刻相当于 0,T0 前后 30 秒就很容易用 -R-30/30/0/2 来表示了;
  • -C 选项很重要。若不使用 -C 选项, -M1 会将整个波形的最大振幅差归一化到 1 英寸,进而导致截取的部分波形最大振幅差比较小;若使用了 -C 选项,表示截取 - 30 秒到 30 秒的波形,然后将这段波形的最大振幅差归一化到 1 英寸;

额外的说明

前面的所有代码中,都使用了 -Ent-Ent0 选项,即表示 Y 轴为文件序号的剖面图。明明只有一个波形,直接画波形振幅图不就好了吗,为什么一定要弄成剖面图呢?

的确,在绘制整个波形的时候,完全可以不使用 -Ent 选项的,此时就是最简单的波形振幅图,此时波形位于 Y=0 处,因而 -R 中 Y 的范围改成 -1/1 或者其他包含 0 的范围即可。

在绘图部分波形的时候,也可以不使用 -Ent0 选项。比如 P 波到时为 600 秒,前后 30 秒也就是 570 到 630 秒,直接把 -R 选项中的 X 范围改成 570/630 即可。但不同的波形 P 波到时不同,所以 X 轴的范围也不同,这种做法每次都要先提取 T0 的值,然后加减 30 秒计算 X 轴范围,相对来说比较麻烦。

绘制部分波形时,更合理也更直观的做法,肯定是指定 T0 为参考时刻。要指定 T0 为参考时刻,就必须使用 -E 选项,而使用 -E 选项就意味着要绘制剖面图。剖面类型有 5 种,当剖面类型为 n 时,Y 值取 0/2 即可;如果使用其他剖面类型,Y 值则要依赖于波形中的头段。因而限定了剖面类型为 n

在绘制整个波形时,使用 -Ent 且 Y 轴范围为 0/2,使得脚本可以只需要尽可能少的修改即可适用于绘制部分波形。