不同比例尺曲线的画法

本文由 Joe Wang(cjmuqiao@163.com)投稿。

目标效果图

问题的提出

在地磁场中 X、Y、Z 三分量幅值两两之间通常相差一两个数量级。如果直接将其画到相同比例尺的一个 坐标系下,幅值小的通常就成了一条直线。此时需要将大幅值适当缩小,或者小幅值适当放大,即在一张图 中绘制不同比例尺的曲线。在丁鉴海、卢振业、余素荣编著的《地震地磁学概论》中经常会看到不同比例 尺的图。

使用 GMT 的 -Jx 选项,绘图时 X 坐标或Y坐标就会根据数据的范围进行缩放,但这种放大或 缩小的显示方式,会带来一个刻度值显示的问题。即:

  1. 若将数据点用线连起来,由于数据点连线的时候不支持 -N 选项,则 -R 的范围必须包括 整个 Y 轴的取值范围,Y 轴的刻度会完全显示,这样有点啰嗦;
  2. 若使用 -S 选项,将数据点用某种符号显示,此时再加上 -N 选项,则可以只显示部分比例尺, 同时又可以把数据完整地表示出来,但是这样的散点图看起来断断续续的,无法连接起来;

解决思路

  1. 使用 psxy 将数据点用线连起来,此时不绘制刻度;
  2. 选择要绘制的刻度的范围,并使用 psbasemap 绘制部分数据的 Y 轴刻度,在此之前需要 仔细计算坐标原点在 Y 轴的偏移量。

脚本

  • 软件环境:GMT5 windows 版,gawk
  • 使用命令:psxy, psbasemap, gmtmath, gawk
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
set sta=CHL
set day=20120101
set ps=%sta%%day%.ps
set Rhour=0/24/0/5
set Rsec=0/86400/0/5
set Jhour=x0.72/0.1
set Jsec=x0.0002/0.1
REM preparing the input data
REM 1. calculate mean value
REM 2. subtract mean value from original data
REM 3. output line number and data
gmt gmtmath CHL2_20120101Hx.dat CHL2_20120101Hx.dat MEAN SUB = | gawk "{print NR, $1}" > x.dat
gmt gmtmath CHL2_20120101Hy.dat CHL2_20120101Hy.dat MEAN SUB = | gawk "{print NR, $1}" > y.dat
gmt gmtmath CHL2_20120101Z.dat CHL2_20120101Z.dat MEAN SUB = | gawk "{print NR, $1}" > z.dat
gmt psxy -R%Rhour% -J%Jhour% -T -K -P >%ps%
REM plot x-axis
gmt psbasemap -R%Rhour% -J%Jhour% -Bxa2f1+lTime/hour -Bya5 -BS -K -O >>%ps%
gmt psbasemap -R%Rsec% -J%Jsec% -Bxa7200f3600+lTime/s -Bya5 -BS -Y2c -O -K >>%ps%
REM plot Z
gmt psxy z.dat -R0/86400/-20/15 -Jx0.0002/0.1 -O -K>>%ps%
REM plot scale of z data
gmt psbasemap -R0/86400/0/5 -J -Y2c -Bya5 -BW+t"@%%12%%\104@%%6%%Z" -K -O >>%ps%
REM plot Hy
gmt psxy y.dat -R0/86400/-400/400 -Jx0.0002/0.005 -Wblue -O -K >>%ps%
gmt psbasemap -R0/86400/100/200 -J -Y2.5c -Bya100+l"nT" -BW+t"@;blue;@%%12%%\104@%%6%%Y" --FONT=blue --MAP_DEFAULT_PEN=+blue --MAP_TITLE_OFFSET=0p -O -K >>%ps%
REM plot Hx
gmt psxy x.dat -R0/86400/-55/30 -Jx0.0002/0.05 -Wred -O -K >>%ps%
gmt psbasemap -R0/86400/10/20 -J -Y3.25c -Bya10 -BW+t"@;red;@%%12%%\104@%%6%%X" --MAP_TITLE_OFFSET=0p --FONT=red --MAP_DEFAULT_PEN=+blue -O -K >>%ps%
gmt psxy -R%Rhour% -J%Jhour% -T -O >>%ps%
del x.dat y.dat z.dat
del gmt.history

一些说明:

  • gmtmath 命令实现是一个逆波兰式的计算器,此处其作用在于先计算数据的均值,然后从原始数据中减去,并将计算得到的结果通过管道传送给 gawk
  • 绘制每个数据时,分别使用了两个命令,其中 psxy 用于绘制数据,而 psbasemap 用于绘制比例尺;
  • 脚本中的参数的选择是有技巧,也是有规律的,下面会说明;
  • 脚本中大部分数值都是可以从数据中自动计算得到的,因而将脚本改一改就可以实现完全的自动化;本文的脚本中没有对此进一步实现,如数据范围,移动距离等都是硬编码在脚本中;
  • 命令中使用了一些额外的选项,来达到更好的绘图效果;

该脚本中某些参数的选取是很有技巧的,以 Z 分量数据为例:

  • psxy 命令中 -R0/86400/-20/15 是 Z 数据的范围,其值可以通过 gmtinfo 命令自动获得;
  • psxy 命令中 -Jx0.0002/0.1,即对于 Y 轴而言,数据的 1,在图上是 0.1 厘米,则 Z 数据的 35 的范围,对应于图上 3.5 厘米;
  • 为了在 Y 轴的 0 到 5 处绘制一个比例尺,需要使用 psbasemap 命令,该命令有两个需要注意的地方:
    • -R0/86400/0/5 设定了 Y 轴范围为 0 到 5;
    • 为了在真正的 0 处绘制比例尺,需要移动坐标系原点,此处的移动距离为 \((0-(-20))*0.1=2cm\),即 -Y2c 选项;

修订历史

  • 2015-02-01:Joe Wang 投稿;
  • 2015-02-01:SeisMan 整理语句、优化脚本;
  • 2015-02-02:使用 FONT 和 MAP_DEFAULT_PEN 代替其他参数;