画地图

昨天简单用 GMT 画了一个研究区域位置图，并不费劲，今天画地形图才算是费了大劲。首先是 pscoast 画的底图跟 DEM 对不上。我在 Google Earth 和 GRASS 之间来回对照没看出什么问题，后来打开经纬度网格才发现 gtopo30 中的青藏高原南边的一系列山峰跑到了北纬30度附近（实际应为27.5度左右），不知道是 gtopo30 这个数据集本身的 bug 还是 GRASS 转换格式时出的问题。

[caption id="" align="alignleft" width="240" caption="Topography map of Tibet"][/caption]

后来重新搜了一下，改用 UC San Diego 的 Smith & Sandwell[1. Smith, W. H. F., and D. T. Sandwell, Global seafloor topography from satellite altimetry and ship depth soundings, Science, v. 277, p. 1957-1962, 26 Sept., 1997.] 的全球地形数据集。之后就是不断地看各个命令的 man page，以及 GMT Cookbook 还有邮件列表和各种网站上的资料，并攻克了诸如调节调色板对应高度范围等种种难题。在这过程中我发现 GMT 虽然学习曲线相当陡峭，但并不是阳春白雪，许多大学的制图课都将其作为重要的内容。

绘图用的脚本如下：

#/bin/bash
# prepare the surface
# data downloaded from ftp://topex.ucsd.edu/pub/global_topo_1min/
img2grd ~/images/global_topo_1min/topo_11.1.img -T1 -S1 -R70/105/25/40 -m1 -D -Gtopo.nc -V
mv topo.nc ~/Desktop/thesis/auxillary/topo.nc
# land is a home-made cpt derived from GMT_relief.cpt
makecpt -Cland -T-1000/9000/1000 -Z -V > topo.cpt
# Relieves
grdgradient ~/Desktop/thesis/auxillary/topo.nc -Ne1 -A100 -M -Glight.nc
#basemap
grdimage ~/Desktop/thesis/auxillary/topo.nc -Ilight.nc -Ctopo.cpt -Jm87.5/32.5/0.6c -Bpa600mg600m -R70/105/25/40 -K > ~/Desktop/thesis/images/topography.ps
pscoast -Bpa600mg600mf300m -J -R -I2 -N1/thick -W -Lf77/27/33/1000 -Tf102/37/1.5c -O -K >> ~/Desktop/thesis/images/topography.ps
psxy ~/Desktop/thesis/auxillary/boundary.xy -O -J -R -Sc0.03c -Gblack -K >> ~/Desktop/thesis/images/topography.ps
# annotations
pstext -O -J -R -K >> ~/Desktop/thesis/images/topography.ps <75 29 18 0 1 LM @#India@%%
80 39 16 0 21 LM XUAR
84 33 16 0 21 LM TAR
93 37 16 0 21 LM Qinghai
100 30 16 0 21 LM Sichuan
END
pstext -O -J -R -Wwhite,Othinnest,white -K >> ~/Desktop/thesis/images/topography.ps <88 37 14 0 5 LM [1]
92 33 14 0 5 LM [2]
89 30 14 0 5 LM [3]
100 29 14 0 5 LM [4]
END
# [1]: Kunlung Mt.
# [2]: Tangula Mt.
# [3]: Niangqen-Tanglha Mt.
# [4]: Hengduan Mt.
# making and drawing contours
grdcontour ~/Desktop/thesis/auxillary/topo.nc -Jm87.5/32.5/0.6c -Q10000 -L-1000/9000 -C500 -O -K >> ~/Desktop/thesis/images/topography.ps
psscale -D20c/4.5c/9c/1c -Ctopo.cpt -Bp2000:@Topography:/:m: -O >> ~/Desktop/thesis/images/topography.ps

# convert format
ps2raster ../images/topography.ps -Tg -P -A

# clean up
rm topo.cpt
rm light.nc