用GMT绘制中国大陆及邻区地质年代图

用GMT绘制绘制地质图一直是一个难题，此前Po-Chin Tseng的博客详细介绍了制作和绘制地质图的过程，并成功绘制了东南亚的地质图。本文将介绍用相似的方法绘制中国大陆及邻区的地质图。

数据来源与格式转换#

地质图的数据来自美国USGS（Generalized Geology of the Far East (geo3al)）其分辨率为1:5,000,000。源数据和引用格式请参考这里

原始数据下载#

我们所需要的文件可以从这里下载，解压该文件后，获得.shp, .dbf, .prj 和 .shx文件。

格式转换#

GMT提供了命令ogr2ogr来进行地理矢量数据格式之间的转换，也包括这里的.shp格式。解压后的坐标不是我们希望的经纬度坐标，而是大地坐标，这里需要使用ogr2ogr中的-t_srs参数进行坐标转换。

1
ogr2ogr -f GMT geo3al.gmt geo3al.shp -t_srs EPSG:4326

NOTE

其中.dbf, .prj 和 .shx文件路径需要与该命令的运行路径为同一路径。

确保在安装GMT后所添加的环境变量中包含PROJ_LIB这个变量，请参考GMT安装说明

格式转换后geo3al.gmt文件的头段为

1
# @D21315149832.578|1546682.830|v|J|J
2
# @P
3
123.121278084975 53.3084236387908
4
123.230414901407 53.3161680969698
5
123.288482121838 53.3166334127133
6
123.328886980316 53.3062394894715
7
......

下载链接#

这里我们提供了已经完成上述步骤的地质年代文件和地质年代色标文件，可以在GMT中直接使用。

GMT脚本#

1
data=geo3al.gmt
2
cpt=geoage.cpt
3
rocksize=500
4
lengsize=0.15i
5

6
function shp2xyz(){
7
    ogr2ogr -f GMT $data geo3al.shp -t_srs EPSG:4326
8
}
9

10
function plot_age_legend(){
11
    cat > tmp << EOF
12
H 10 3 Age of rock units
13
G 1p
14
N 3
15
EOF
16
    awk '!/^($|B|F|N|#)/{print $0}' $cpt | while read label color period
17
    do
18
        echo "S 0.3c r $lengsize $color 0.3p 0.7c $period" >> tmp
19
    done
20
    gmt legend tmp -DJBR+w300p/140p+jBR+o0c/-83p+l1.3 -F+p0.7p+g255 -C3p/3p
21

22
    gmt legend -DJBR+w150p/50p+jBR+o0c/57p+l1.9 -F+p0.7p+g255 -C3p/1p <<EOF
23
H 10 3 Rock type
24

25
N 2
26
S 0.3c r 0.2i p400/28:B255 0.3p 0.7c Volcanic rocks
27
S 0.3c r 0.2i p400/29:B255 0.3p 0.7c Intrusive rocks
28
S 0.3c r 0.2i p400/44:B255 0.3p 0.7c Ultrabasic igneous rock or ophiolites
29
EOF
30
}
31

32
shp2xyz
33
gmt begin geo3al png
34
    gmt set FONT_ANNOT_PRIMARY 10
35
    gmt set MAP_FRAME_WIDTH 0.08
36
    gmt set MAP_TICK_LENGTH_PRIMARY 0.08
37

38
    gmt coast -R70/150/13/55 -JM8.6i -Baf -Df -G255 -BWsNe
39
    gmt plot $data -C$cpt -aZ="GEN_GLG" -G+z -L
40
    gmt convert $data -aI="TYPE" -S"-Iv" | gmt plot -Gp28+r500+f100+b-
41
    gmt convert $data -aI="TYPE" -S"-Ii" | gmt plot -Gp29+r500+f100+b-
42
    gmt convert $data -aI="TYPE" -S"-Iw" | gmt plot -Gp44+r500+f100+b-
43

44
    gmt coast -SCADETBLUE1
45

46
    gmt set FONT_ANNOT_PRIMARY 7p
47
    plot_age_legend
48

49
    rm tmp  tmp.cpt
50
gmt end show

特别鸣谢#

这次脚本做了大量的修改，非常感谢田冬冬博士、姚家园博士和GMT中文社区的其他小伙伴的建议和努力，在和他们的讨论中也使我受益匪浅。