在地图上绘制波形
文章目录
将 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
先绘制地图
在考虑如何把波形画在地图上之前，先得有一个地图。用如下代码绘制一个地图：
#!/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
代码中使用了 saclst 和 awk 提取了每个波形所对应的台站经纬度，并通过管道传递给 psxy 和 pstext
命令。
地图效果如下：
第一个波形振幅图
以 ntkl.z 的波形为例，先不去想如何把波形放在台站附近，只考虑如何将波形画在地图上。
先把代码写出来，如下：
#!/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，使得将波形偏移到台站坐标处的同时不修改原底图的坐标原点，这样不会影响到后面波形的绘制。
#!/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 的输出的基础上加上额外的偏移量，波形就可以放在台站正右方了。
代码如下：
#!/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 的命令会出现问题，请勿使用。
#!/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 波到时）：
几点说明：
-Ent0中n表示绘制 Y 轴为文件序号的剖面图，因为每个 pssac 命令中只绘制一个波形，因而这个波形会被画在Y=1处，因而此处-R中 Y 轴的范围取的是 0 到 2。实际上只要 Y 轴的范围包含 1 即可，比如 Y 轴范围取0.5/1.5、0/3均可，可以尝试取不同的值查看其效果；-Ent0中t0表示以 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，使得脚本可以只需要尽可能少的修改即可适用于绘制部分波形。
文章作者 SeisMan
上次更新 2015-07-18