geo3al: 中国及邻区地质图数据

贡献者: 徐弥坚, 陈箫翰, 刘珠妹, 田冬冬, 姚家园

geo3al 数据是由 U.S. Geological Survey (USGS) 提供的中国及邻区地质图数据，是“世界地质地图”的一部分，数据分辨率为 1:5,000,000。

数据文件

中文社区提供了可供 GMT 直接使用的数据文件：

geo3al.gmt：中国及邻区地质图数据

数据的具体下载地址和使用方式请见：中国地理空间数据集。

数据中包含了如下非空间元信息（可以参考 geo3al.shp.xml 以获取元信息的更多细节）：

AREA：地层多边形区域的面积（单位 $m^2$）
PERIMETER：地层多边形区域的周长（单位 m）
TYPE：岩性
GLG：原始地图上显示地表露头的地质年代
GEN_GLG：编者从 GLG 中计算得到的地质年代，在单元年龄跨越一个以上地质年龄的情况下用“代”表示

示例图

下面的示例图使用 geo3al.gmt 数据绘制了中国大陆及邻区地质图：

地质年代从数据的 GEN_GLG 属性中提取，用不同的颜色表示
岩性从数据的 TYPE 属性中提取，用不同的填充图案表示

示例代码中使用的地质年代 CPT 色标文件 geoage.cpt 修改自 GTS2012_epochs。

#!/usr/bin/env bash
#
# Geological map of China and adjacent regions
#

data=geo3al.gmt
cpt=geoage.cpt

gmt begin geo3al pdf,png
    # Paint the land using white color
    gmt coast -R70/150/13/55 -JM22c -Gwhite -Baf -BWsNe

    # Plot the geological map, colored by geologic ages
    # -aZ="GEN_GLG": use the "GEN_GLG" metadata as the Z value
    # -G+z -C$cpt: the fill color is determined by the Z value and the CPT file
    gmt plot $data -C$cpt -aZ="GEN_GLG" -G+z

    # Plot rock types using diffenrent patterns (28, 29 and 44)
    # -aI="TYPE": using the "TYPE" metadata as ID
    # -S"-Iv": only keep data segments that matches "-Iv" (v means "Volcanic rocks")
    # -S"-Ii": only keep data segments that matches "-Ii" (i means "Intrusive rocks")
    # -S"-Iw": only keep data segments that matches "-Iw" (w means "Ultrabasic igneous rock or ophiolites")
    # -Gp28+r500+f100+b-: fill with pattern 28, dpi=500, foreground color=100 (gray), background color=- (transparent)
    gmt convert $data -aI="TYPE" -S"-Iv" | gmt plot -Gp28+r500+f100+b-
    gmt convert $data -aI="TYPE" -S"-Ii" | gmt plot -Gp29+r500+f100+b-
    gmt convert $data -aI="TYPE" -S"-Iw" | gmt plot -Gp44+r500+f100+b-

    # Paint the ocean and lake using color "cadetblue1"
    gmt coast -SCADETBLUE1

    # Plot geologic age legend
    gmt set FONT_ANNOT_PRIMARY 7p  # font for legend labels
    cat > age_legend.txt << EOF
H 10 3 Age of rock units
N 3
EOF
    # Output non-empty lines that don't start with "B", "F", "N" or "#"
    gawk '!/^($|B|F|N|#)/' $cpt | while read label color period
    do
        echo "S 0.3c r 0.38c $color 0.3p 0.7c $period" >> age_legend.txt
    done
    gmt legend age_legend.txt -DJBR+w10.5c+jTR+o0c/-2c+l1.3 -F+p0.7p+g255 -C3p/3p
    rm age_legend.txt

    # Plot rock type legend
    gmt legend -DJBR+w5.0c+jBR+o0c/2c+l1.9 -F+p0.7p+g255 -C3p/1p <<EOF
H 10 3 Rock type
N 2
S 0.3c r 0.5c p28+r500+f100+b255 0.3p 0.7c Volcanic rocks
S 0.3c r 0.5c p29+r500+f100+b255 0.3p 0.7c Intrusive rocks
S 0.3c r 0.5c p44+r500+f100+b255 0.3p 0.7c Ultrabasic igneous rock or ophiolites
EOF
gmt end show

https://user-images.githubusercontent.com/3974108/123166472-7376e080-d443-11eb-9289-d8913afed442.png

数据来源与处理

以下内容仅供数据维护者参考，一般用户无需关心。

从 “World Geologic Maps”下载 “Generalized Geology of the Far East” 对应的 Shapefile 压缩文件
使用 ogr2ogr 将 Shapefile 格式转换为 GMT 可识别的 OGR/GMT 格式:
```
$ ogr2ogr -t_srs EPSG:4326 -f OGR_GMT geo3al.gmt geo3al.shp
```
- 原始数据使用的是以米为单位的投影坐标，而非以经纬度表示的地理坐标。因而在使用 ogr2ogr 时需要加上 -t_srs EPSG:4326 参数进行坐标转换
- 转换时 .dbf、.prj、.shx 等文件必须和 shp 文件放在同一目录下

数据转换与处理参考了如下博文：

Mijian Xu 的博客：用 GMT 绘制中国大陆及邻区地质年代图
Po-Chin Tseng 的 GMT5 教程：地质图

扩展阅读

World Geologic Maps 提供了几乎所有区域的地质图数据。本文只提供了中国及邻区地质图（即 Generalized Geology of the Far East）。需要绘制其他地区地质图的用户可以参考本文的数据处理方式和绘图脚本，但需要注意：

部分原始数据使用的是经纬度坐标，因此格式转换时不需要使用 -t_srs EPSG:4326
每个数据提供的属性可能不同，因而用户应查阅从源文件中解压得到的 .xml 文件以了解数据提供的元信息，并相应地修改绘图脚本。

数据引用

Steinshouer, D.W., Qiang, J., McCabe, P.J., and Ryder, R.T, 1999, Maps showing geology, oil and gas fields, and geologic provinces of the Asia Pacific region: U.S. Geological Survey Open-File Report 97-470-F, 16 p., https://doi.org/10.3133/ofr97470F.