檀香山是夏威夷主岛火奴鲁鲁的旧称,因孙文为代表的革命党人在这里的革命活动而为华人熟知。

现在我们想绘制檀香山的地形。我们希望地形的精度尽可能的高,以得到最清晰的地图。 目前,免费开放的最精确的地形文件的精度是 1 弧秒,即 30m 一个采样点,但只包含陆地。包含了海洋的地形的精度最高是 15s。 这样,我们就需要用 1 弧秒的数据来绘制地形,15 弧秒的数据来绘制海底。

思路:首先用 15 弧秒的文件绘制檀香山及周围海域(既有陆地,又有海洋)的地形。 然后,用 1 弧秒的数据去绘制陆地部分的地形。因为 1s 的数据后绘制,就把前面不太准确的 15 弧秒陆地地形覆盖了。 要实现这个思路的关键在于绘制 1 弧秒数据的时候,怎样才能只覆盖陆地部分,而保留海洋部分。 实现方法是:gmt 的 pscoast 模块提供了一个 clip 的功能,实现用海岸线进行裁剪。具体请看代码中的注释:

#!/bin/bash
PS=hawaii.ps
landdata=hawaii_land_1s.grd
oceandata=hawaii_ocean_15s.grd
cpt=hawaii.cpt

# 设置光照
gmt grdgradient $landdata -Ne0.7 -A315 -G$landdata.light
gmt grdgradient $oceandata -Ne0.7 -A315 -G$oceandata.light
# 底图
gmt psbasemap -R-156.5/-154.5/18.5/20.5 -JM15c -Baf -BWSEN -K -P > $PS
# 地形
#  海洋部分
gmt grdimage $oceandata -R -J -I$oceandata.light -C$cpt -O -K >> $PS
#  陆地部分
gmt pscoast -R -J -Da -Gc -O -K >> $PS # 用最高精度(-Da)的海岸线划出裁剪范围(-Gc)
gmt grdimage $landdata -R -J -I$landdata.light -C$cpt -O -K >> $PS # 绘制陆地部分地形
gmt pscoast -O -K -Q >> $PS # 关闭裁剪
# 文件尾
gmt psxy -R -J -T -O >> $PS
# 输出
gmt psconvert $PS -A -Tg
rm $PS *.light gmt.history

绘图效果如下:

夏威夷檀香山