提出问题

需要画一个沿着某个经线或纬线的深度剖面图,即 X 轴表示经度或纬度,Y 轴表示深度,如下图:

分析问题

这张图比较特殊的地方在于,X 轴是地理坐标,Y 轴是笛卡尔坐标。最简单的办法是直接用线性投影 -JX 来绘制整张图,比如:

gmt psbasemap -R40/50/0/600 -JX15c/-10c -Bx2 -By100 > mix.ps

效果如下图:

显然,这张图还差了一点,X 轴的坐标都没有 “度” 符号以及 WSEN 后缀,无法体现出其地理坐标的特征。

解决问题

解决办法如下:

gmt set FORMAT_GEO_MAP +ddd:mmF
gmt psbasemap -R40/50/0/600 -JX15cd/-10c -Bx2 -By100 > mix.ps

几点说明:

  • FORMAT_GEO_MAP 用于设置地理坐标的显示方式。这里的 +ddd:mmF 表示以 度:分 的形式显示,并加上 EWSN 后缀;
  • -JX15cd/-10c 设定了线性投影方式。X 轴多了一个 d,用于显式指定 X 轴为地理坐标;Y 轴没有 d,为正常的笛卡尔坐标;这里 Y 轴的长度是负值,表示 Y 值从上到下递增,以符合常见的深度的定义;
  • GMT4 同理,但 GMT 4.5.13 似乎存在 bug,当使用 -JX15cd/-10c 时,X 轴的标注的位置会出现偏差;尚未发布的 GMT 4.5.14 中该 bug 已被修复;
  • X 轴由于是地理坐标,所以不能用 -B 选项给 X 轴加标题;如果需要加标题的话,只能使用 pstext

其它解决方法

下面列出最初对于这个问题的分析以及两种稍复杂的解决方案。虽说下面的两种方案更复杂了,但其思路与想法还是很有意思的,或许在其它问题上可以借鉴,因而将其列出。

分析问题

  1. 因为 Y 轴是线性坐标,所以必然只能选择线性投影,即 -JX
  2. 线性投影的情况下,图的主体很简单,关键在于 X 轴坐标的 “度” 符号以及后缀 E 上
  3. 尝试将 X 轴和 Y 轴都当作线性坐标,然后对于 X 轴设置其 单位 为特殊的 “度” 符号。此法看似可行,但实际上 GMT 内部设置了单位与标注之间的距离,通过单位设置的 “度” 符号明显离标注的距离较远,不太美观。
  4. 为了使 X 轴有一个 “度” 符号,X 轴必须当作地理坐标处理;而 Y 轴的范围已经超过了地理坐标的范围,所以必须当作线性坐标处理;
  5. 综上,必须使用两个命令来绘制边框,分别绘制地理坐标的 X 轴(-BxxxSN)和线性坐标的 Y 轴(-BxxxEW);
  6. 关于如何绘制地理坐标的 X 轴,有两种解决办法。

解决问题

法 1

#/bin/bash
Rx=40/50    # X 轴范围
Ry=0/600    # Y 轴范围
R=$Rx/$Ry
B=2/100     # 间隔
J=X15c/10c
PS=map.ps

psxy -R$R -J$J -T -K > $PS   # 写入 PS 文件头

psbasemap -R$R -J$J -B${B}SN -K -O --D_FORMAT='%g\260E' >> $PS    # 绘制 X 轴
psbasemap -R$R -J$J -B${B}EW -K -O >> $PS     # 绘制 Y 轴

# 这里放置其它绘图命令,不再使用 B 选项

psxy -R$R -J$J -T -O >> $PS  # 写入 PS 文件尾
rm .gmt*

这里绘制 X 轴时直接使用 --D_FORMAT=%g\260E ,使得在该命令中 D_FORMAT 的值为 %g\260E, 即设置显示浮点数时在其后加上 “度” 符号以及后缀“E”。

此法的优点在于简单,缺点在于后缀 “E” 是固定值,无法处理东西经同时存在的情况。

法 2

#/bin/bash
Rx=40/50    # X 轴范围
Ry=0/600    # Y 轴范围
Rfake=0/1   # 假轴范围
R=$Rx/$Ry
B=2/100     # 间隔
J=X20c/15c
PS=map.ps

psxy -R$R -J$J -T -K > $PS   # 写入 PS 文件头

psbasemap -Rg$Rx/$Rfake -J$J -B${B}SN -K -O --BASEMAP_TYPE=plain >> $PS    # 绘制 X 轴
psbasemap -R$R -J$J -B${B}EW -K -O >> $PS     # 绘制 Y 轴

# 这里放置其它绘图命令,不再使用 B 选项

psxy -R$R -J$J -T -O >> $PS  # 写入 PS 文件尾
rm .gmt*

其它的说明

  1. GMT 的 B 选项提供了这样一个功能,可以使用形如 -Bgxmin/xmax/ymin/ymax 的语法,其中 g 告诉命令即便使用 -JX 投影,也认为其是地理坐标。由于是地理坐标,“度”符号以及后缀 “E” 就很容易出来了
  2. 使用 -Bgxmin/xmax/ymin/ymax 存在两个问题

    1. 虽然是线性投影,但是由于使用了地理坐标,GMT 会默认将底图类型设置为 fancy。这一点可以设置 BASEMAP_TYPE 等于 plain 来解决。
    2. Y 轴被当作地理坐标,所以 ymin 和 ymax 的范围被限制在 [-90,90] 之内
  3. 在此例中在绘制 X 轴时引入了一个假的 Y 轴 0/1,以满足 -Rgxmin/xmax/ymin/ymax 形式中对 ymin 和 ymax 范围的限制。

这样,X 轴和 Y 轴就都设计好了,接下来要做的就只是保证其它命令都不使用 B 选项即可。

修订历史

  • 2014-04-26:初稿;
  • 2014-04-26:修改脚本,解决了对 Y 轴范围的限制;Thanks to Chen Zhaohui;
  • 2014-06-09:通过修改 D_FORMAT 以解决地理坐标的度符号;该方法由刘珠妹提供;
  • 2014-07-09:找到了一种非常简单的方法来解决该问题;
  • 2014-11-24:修正了 --PLOT_DEGREE_FORMAT 中的小问题;
  • 2015-06-29:GMT4.5.13 在解决该问题时有其他 bug,这里使用 GMT5;