上面有篇博文已经介绍了绘制中国地图的方法,这里主要介绍一下用GMT绘制一张嵌入世界地图的中国地图解决方法。权当做笔记整理。

  • OS Mac OSX 10.11.2
  • GMT version GMT 5.4.2
  • Script language bash

数据准备

全球地形起伏数据

鉴于绘制的世界地图只是作为嵌入图，因此数据分辨率不需要太高，这里选择SRTM30高程数据。 直接下载netCDF格式的网格数据(约1.7G)，该网格数据可直接用GMT进行成图。

网格数据重采样

GMT提供了grdsample命令进行网格数据重采样，可以进行数据加密，也可以进行数据降采。前面提过这里的世界地图不需要太高的分辨率，因此将原始30’'(约1km)的数据降采10km的网格数据进行全球地图成图，降采到2km的网格数据进行中国地图成图。

1
2
3
4

# -I选项确定重采样空间步长
# -R选项确定新网格数据分布，这里不想把欧洲分开，因此没用0-360
gmt grdsample -R-25/335/-60/60 -I10k topo30.grd -Gtopo_10k.grd
gmt grdsample -R70/135/15/55 -I2km topo30.grd -GChina_Topo.grd

绘制图中图

绘制中国地形图，外加世界地图作为图中图嵌入在地图左下角。绘图要点：

  • 中国地形图
  • 海岸线及行政区划
  • 世界地图图中图
 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
40
41
42

#---------Initial Standard Options
# OPTIONS FOR CHINA MAP
J1=X6id/4id
R1=70/135/15/55
B1=15/10WSEN 
# OPTIONS FOR INSIDE MAP
J2=X2.8id/1.5id
R2=-25/335/-60/60
B2=0
#
PS=map_topo-1.ps
input1=China_topo.grd
input2=topo_10k.grd
#---------DATA PREPARE
# gmt grdsample -R$R1 -I2km topo30.grd -G$input1
# gmt grdsample -R$R2 -I10k topo30.grd -G$input2

#---------OPEN GMT
gmt psxy -R$R1 -J$J1 -T -K > $PS

#---------Plot China Map
# 计算梯度文件，便于后面成图。
gmt grdgradient $input1 -Gintense.grd -A0/270 -Nt1.0 -V
# 制作cpt文件。
gmt makecpt -Cglobe -T-10000/10000/100 -Z -V >pstopo.cpt
# grd文件成图。
gmt grdimage -J$J1 -X4.5 -Y5.5 $input1 -Iintense.grd  -Cpstopo.cpt -O -K -V >>$PS
# 绘制海岸线，国界边界数据。
gmt pscoast -R$R1 -J$J1 -B1$B1 -O -K -Dh -A300  -N1/0.01p -P >>$PS

#---------Plot Inside Map
gmt grdgradient $input2 -Gintense.grd -A0/270 -Nt1.0 -V
gmt grdimage -J$J2 $input2 -Iintense.grd  -Cpstopo.cpt -O -K -V >>$PS
# 基于绘制海岸线，国界边界数据。
gmt pscoast -R$R2 -J$J2 -B$B2 -O -K -Dh -A300 -N1/0.01p -P --MAP_FRAME_TYPE=plain >>$PS
gmt pscoast -R$R2 -J$J2 -B$B2 -Df -N1 -W -A5000 -K -O --MAP_FRAME_TYPE=plain >> $PS
# 利用 psbasemap 在大区域地图内绘制小区域所在的方框
gmt psbasemap -R$R2 -J$J2 -D$R1 -F+p2p,black -K -O --MAP_FRAME_TYPE=plain >> $PS
#---------Plot STATIONS

#---------CLOSE GMT
gmt psxy -R$R -J$J -T -O >> $PS

绘制效果如下： psMAP-1

这里说明以下几点：

  1. B1，这里仍然用gmt4的语法，不过gmt5仍然兼容
  2. 注意养成良好的习惯，GMT开关。
  3. 制作cpt文件时，-T-10000/10000/100控制色标棒上下界值
  4. 生成中国地形图时,-X4.5 -Y5.5控制图片位置，方便后续修改色标棒位置或者gmt logo
  5. 绘制inside map时，语句上没有什么变化，紧紧是控制不同的标准选项而已
  6. inside map默认绘制到左小角，也可以通过-X -Y选项来手动修改
  7. --MAP_FRAME_TYPE放在句中为局部变量，并不会更改其他语句中的环境设置

绘制台站

绘制台站标识主要使用PSXY，通过更改不同option来设置不同的标识符。

1
2
3

#---------Plot STATIONS
awk '{print $1, $2}' psStation.txt | gmt psxy -J$J1 -R$R1 -B$B1 -Sc0.15c -W1p,black -Glightgray -K -O >> $PS 
awk '{print $1, $2}' psOUT.txt | gmt psxy -J$J2 -R$R2 -B$B2 -St0.2c -W1p,black -Gred -K -O --MAP_FRAME_TYPE=plain >> $PS

绘制效果如下： psMAP-2

这里需要注意：

  1. 绘制嵌入图中的台站标识时与绘制图中图一样，采用相应的STANDARD OPTIONS即可相应控制。
  2. psxy的输入为[X Y]坐标矢量，可以通过awk从台站列表文件中拾取，$1 $2分别表示第一列和第二列，注意两列之间的分隔符,(让第一列和第二列分隔打印出来)。
  3. -S选项控制标识符，c表示圆圈，t表示三角形,0.15c表示半径大小。
  4. -W选项控制字体大小或者曲线粗细及颜色。
  5. -G选项控制填充颜色

添加台站名

添加台站名使用PSTEXT,如下：

 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

#---------Plot STATIONS NAME
gmt pstext -R$R1 -J$J1 -B$B1 -K -O -N  -F+f10p,2,black >> $PS << EOF
125  31.14  Shanghai  
109 31.4   Enshi 
120.17 36.9   Qingdao  
113.7  31.40  Wuhan  
114   35.50  Zhengzhou  
84.51   46.20  Karamary 
93.28   41.7   Hami 
87.36   44.29  Urumqi   
116.6  40.4  BeiJing  
103.51  36.99  Lanzhou  
108.57  34.7  Xi'an  
103.79  25.30   Kunming 
114  38.0  Zhuozhou  
93.03 44.60  Balikun  
121 29  Hangzhou 
119 31.31  Suzhou 
121.8 32.8  Nantong  
113 27 Changsha  
110.7 29.0 Yueyang 
122 39.9 Dalian  
EOF
gmt pstext -R$R2 -J$J2 -B$B2 -K -O -N --MAP_FRAME_TYPE=plain -F+f8p,blue+jBL -Glightgray >> $PS << EOF
3  36   Delft 
3  28  Cologne 
3  20  Karlsruhe
EOF
gmt pstext -R$R2 -J$J2 -B$B2 -K -O -N --MAP_FRAME_TYPE=plain -F+f8p,red+jBL -Glightgray >> $PS << EOF
-160  16  SanFranciso 
-160  8     Boise
-160  0    Denver
-160  -8 Oklahoma  
-160  -16   New Orleans   
-160  -24 Princeton
EOF

这里需要注意：

  1. 没有使用awk直接从列表导入，是因为想手动修改台站名的相应位置，因为部分台站相距太近。当然我们也可以使用awk '{print $1+1.3, $2+0.2}, $3' psStation.txt |来批量修改。
  2. -F选项控制添加文本属性，-f控制字体，8p表示字号，red控制颜色，j控制文本排版对齐方式，可参考文本属性。

添加GMT logo及图名

这里仍然使用PSTEXT，用法与上文一致。

1
2
3

gmt pstext -R$R1 -J$J1 -B$B1 -K -O -N  -F+f20p,28,red -UBR/13.8C/0.1C/'F. Cheng' >> $PS << EOF
105 53 Where have we been?
EOF

这里需要注意：

  1. -U控制GMT时间戳，BR表示右下角对齐，/13.8c/0.1c控制相应位置，具体参考U属性。
  2. -F选项中多了个28，表示GMT内置28号字体，可参考文本属性。

PS转JPG

可以使用psconvert命令进行图像转换，亦可参考Seisman博客图片转换。

1
2
3
4

#---------CLOSE GMT
gmt psxy -R$R -J$J -T -O >> $PS
#---------PS2JPG
gmt psconvert -A -E300 -P $PS -TG

最终绘制效果如下： psMAP-3

Reference及推荐阅读