【记录】GMT6绘图示例

yep吖

已于 2024-01-03 16:37:41 修改

阅读量2.2k

点赞数 26

文章标签： linux

于 2023-12-04 23:35:13 首次发布

本文链接：https://blog.csdn.net/weixin_44259058/article/details/130757117

版权

【注】部分GMT绘图示例参考于GMT中文社区及博客，并基于Linux的gmt6版本成功运行绘图脚本，用于记录学习gmt的部分启蒙&有用例子

在这里插入图片描述

1.1折线图+特殊图案
1.2全国地形图+震级分布
1.3由grd文件制作等高线+测线高程
1.4将三维地形、卫星图片以及断层线绘制在不同高度
1.5将卫星图片和断层线直接绘制在三维地形上
1.6地震剖线+MT图
1.7南加州地震分布
1.8海岸线+比例尺+指北针
1.9GPS场
1.10中文字体
1.11板块构造
1.12构造+地形
1.13板块边界
1.14中国及邻区地质图

1.1折线图+特殊图案

#!/usr/bin/env bash
cat > 1.txt << EOF
0 6.680
1 1.18
2 0.861
3 0.533
4 0.374
EOF

cat > 2.txt << EOF
0 6.458
1 1.245
2 0.789
3 0.523
4 0.304
EOF

cat > 3.txt << EOF
0 6.627
1 3.223
2 0.650
3 0.507
4 0.289
EOF

cat > 4.txt << EOF
0 6.680
1 1.105
2 0.547
3 0.532
4 0.257
EOF

cat > 5.txt << EOF
0 6.458
1 6.248
2 0.898
3 0.370
4 0.374
EOF

gmt begin exp1 png
gmt set MAP_GRID_PEN 0.1p,DIMGRAY
gmt basemap -R0/5/0/7 -JX16c/12c -Bxa1g1+l"The Number of Iteration" -Bya1g1+l"normalized RMS" -BWS
gmt basemap -R0/5/0/7 -JX16c/12c -B0 -Ben

# 在不使用 -S 选项的情况下，默认会将所有的数据点连成线，使用 -S 选项则仅在数据点所在位置绘制符号。
# symbol space legend type length (-) thickness,color space LegendName
echo S 0.4c s 0.3c - 0.8p,DEEPSKYBLUE     >> temp.dat
echo S 0.4c c 0.2c - 0.8p,BROWN1  >> temp.dat
echo S 0.4c a 0.3c - 0.8p,PURPLE    >> temp.dat
echo S 0.4c i 0.3c - 0.8p,black    >> temp.dat
echo S 0.4c x 0.3c - 0.8p,RED1  >> temp.dat

echo G -5l >> temp.dat # 上移5行

echo S 0.4c - 0.8c - 0.8p,DEEPSKYBLUE 1.0c MT-single >> temp.dat
echo S 0.4c - 0.8c - 0.8p,BROWN1 1.0c CSEM-single  >> temp.dat
echo S 0.4c - 0.8c - 0.8p,PURPLE  1.0c Joint >> temp.dat
echo S 0.4c - 0.8c - 0.8p,black 1.0c Joint-MT >> temp.dat
echo S 0.4c - 0.8c - 0.8p,RED1 1.0c Joint-CSEM >> temp.dat
gmt legend temp.dat  -DjTR+w3.5c+o0.2c/0.2c -F+p1p

# legend1
gmt plot 1.txt -Ss0.4c -W0.8p,DEEPSKYBLUE 
gmt plot 1.txt -W0.4p,DEEPSKYBLUE

# legend2
gmt plot 2.txt -Sc0.4c -W0.8p,BROWN1
gmt plot 2.txt -W0.4p,BROWN1 

# legend3
gmt plot 3.txt -Sa0.4c -W0.8p,PURPLE 
gmt plot 3.txt -W0.4p,PURPLE

# legend4
gmt plot 4.txt -Si0.4c -W0.8p,black 
gmt plot 4.txt -W0.4p,black

# legend5
gmt plot 5.txt -Sx0.4c -W0.8p,RED1 
gmt plot 5.txt -W0.4p,RED1 

rm temp.dat
gmt end show

1.2全国地形图+震级分布

dep.dat下载🔗

#!/bin/bash

R=70/135/15/55
eqfile=eq.dat
topodata=@earth_relief_06m

gmt begin topo png A+m1c/1c/1c/1c
gmt grdcut $topodata -R$R -GcutTopo.grd
gmt grdgradient cutTopo.grd -Ne0.7 -A50 -GcutTopo_i.grd
gmt grd2cpt cutTopo.grd -Cglobe -T-10000/10000/200 -Z -D
#绘制底图
gmt set FORMAT_GEO_MAP=ddd:mm:ssF
gmt basemap -R70/135/15/55 -JM7i -Ba10f5 -BWesN
gmt grdimage cutTopo.grd -IcutTopo_i.grd -Q
gmt coast -Dh -W1/0.2p -I1/0.25p -N1/0.5p

#绘制colorbar
gmt colorbar -DjCB+w7i/0.15i+o0/-0.5i+m -Bxa2000f400+l"Elevation/m" -G-8000/8000

#分震级绘制地震
awk '$3>=5.0 && $3<6.0 {print $1,$2,$3*0.04}' $eqfile | gmt plot -Sc -Gblue
M5=`awk '$3>=5.0 && $3<6.0' $eqfile | wc -l | cut -d" " -f1`
awk '$3>=6.0 && $3<7.0 {print $1,$2,$3*0.04}' $eqfile | gmt plot -Sc -Gred
M6=`awk '$3>=6.0 && $3<7.0' $eqfile | wc -l | cut -d" " -f1`
awk '$3>=7.0 && $3<8.0 {print $1,$2,$3*0.04}' $eqfile | gmt plot -Sa -Ggreen
M7=`awk '$3>=7.0 && $3<8.0' $eqfile | wc -l | cut -d" " -f1`
awk '$3>=8.0 {print $1,$2,$3*0.06}' $eqfile | gmt plot -Sa -Gpurple -W0.4p,black
M8=`awk '$3>=8.0' $eqfile | wc -l | cut -d" " -f1`

#绘制图例
gmt legend -DjBR+w1.2i+l1.2+o0 -F+g229+p0.25p <<EOF
G 0.01i
H 8 4 MAGNITUDE
C blue
S 0.1i c 0.20 blue 0.25p,blue 0.18i 5~5.9($M5)
C red
S 0.1i c 0.24 red 0.25p,red 0.18i 6~6.9($M6)
C green
S 0.1i a 0.42 green 0.25p,black 0.18i 7~7.9($M7)
C purple
S 0.1i a 0.48 purple 0.25p,black 0.18i 8~8.9($M8)
EOF

rm cutTopo*.grd
gmt end show

1.3由grd文件制作等高线+测线高程

maunakea.grd下载🔗

#!/bin/bash
#grdcontour 模块可以根据地形网格文件绘制等高线，并进行标注。
#project 模块可以根据测线的起点与终点的经纬度坐标和步长，输出测线上采样点的坐标。
#grdtrack 模块则可以根据上述输出的采样点坐标从地形文件中得到采样点的距离和高程。

gmt begin maunakea png A+m1c/1c/1c/1c
# ==== 预设参数 ====
gmt gmtset MAP_FRAME_TYPE=plain \
    FONT_ANNOT_PRIMARY=14p,25,black \
    FONT_LABEL=14p,25,black \
    FONT_TITLE=24p,25,black \
    MAP_GRID_PEN_PRIMARY=thick,black,-- \
    MAP_TICK_PEN_PRIMARY=thick,black \
    MAP_TICK_PEN_SECONDARY=thick,black \
    FORMAT_GEO_MAP=dddmmF

# ==== 等高线图 ====
# 主要等高线 500m 间隔(-A)，次要等高线 100m 间隔(-C)，Y 把地图上移 10 厘米，为剖面图预留空间
gmt grdcontour maunakea.grd -R155:40W/155:16W/19:41N/19:57N -JM15c -Y10c -C100 -Q100 -A500+f8p,25,darkred+o

# ==== 绘制剖面的测线 ====
# 给定测线的起点(-C)和终点(-E)的经纬度，以及步长(-G)，project 输出经纬度坐标通过管道传输给 grdtrack
# grdtrack 从地形文件中取样，把获得的信息输出到 gmt.xy 文件，gmt.xy 文件四列，依次是经度、纬度、距离和高程
gmt project -C-155.65/19.75 -E-155.28/19.9 -G0.1 -Q |gawk '{ print $1,$2,$3 }' > lon-lat.txt
gmt grdtrack -Gmaunakea.grd lon-lat.txt > gmt.xy
# psxy 使用 gmt.xy 的前两列画测线
gmt plot gmt.xy -Wthick,darkred
# 测线的起点和终点标记为 A 和 B
gmt text -F+j+f14p,25,darkred << EOF
-155.65 19.75 LT A
-155.28 19.9 LT B
EOF
gmt basemap -Ba5m -B+t"Mauna Kea Volcano"

# ==== 剖面图 ====
# 使用 gmtinfo 从地形文件中获得 -R 的适当范围
# 使用 -Y 把剖面图下移 9 厘米
gmt plot gmt.xy -i2,3 $(gmt gmtinfo gmt.xy -i2,3 -I1/100) -JX15c/6c -Wthick,darkred -Y-9c

# 先画 NE 两边的外框，再画 WS 两边的外框
gmt basemap -BNE -Bxa5f1 -Bya1000f500g1000
gmt basemap -BWS -Bxa5f1+l"Distance of A-B profile (km)" -Bya1000f500+l"Height (m)"
#rm gmt.xy
gmt end show

1.4将三维地形、卫星图片以及断层线绘制在不同高度

#!/usr/bin/env bash
# 将三维地形、卫星图片以及断层线绘制在不同高度

gmt begin 3D_earth_relief_1 png
    gmt set MAP_FRAME_TYPE plain
    # 1. 绘制地形图形，图层默认位于 z 轴底部（此例中即 -1000 的高度）
    gmt grdview @earth_relief_06m -R101/108/35/42/-1000/20000 -JM10c -JZ10c -N-1000+ggray -Qi -I -Ba -BwsENZ -p60/25    
    gmt colorbar -C -Ba -DJTC+o0/1c -p

    # 2. 绘制断层线图层，-p 设置图层抬高到10000 
    gmt plot CN-faults.gmt -W2p,brown -p60/25/10000 -Ba -BWSen
    gmt basemap -TdjLT+w1.5c+l+o1c -p
    
    # 3. 绘制卫星图片图层，-p 设置图层抬高到20000
    gmt grdimage @earth_day_06m -p60/25/20000 -Ba -BWSen


gmt end show

1.5将卫星图片和断层线直接绘制在三维地形上

#!/usr/bin/env bash
# 将卫星图片和断层线直接绘制在三维地形上

gmt begin 3D_earth_relief_2 pdf,png
    gmt set MAP_FRAME_TYPE plain
    # 预处理1：生成断层的三维采样，输出数据为断层的经度、纬度以及对应的地形高程
    gmt grdtrack -R101/108/35/42 -G@earth_relief_06m CN-faults.gmt > faults.xyz

    # 预处理2：截取区域卫星图片
    # 下载 6 弧分的卫星图片数据 earth_day_06m_p 到当前目录
    gmt which -Gl @earth_day_06m_p
    # 使用 gdal_translate 截取区域的卫星图片
    # -projwin 后的四个参数指定了区域范围：左上角经度 左上角纬度 右下角经度 右下角纬度
    # gdal_translate 的详细用法见 https://gdal.org/programs/gdal_translate.html
    gdal_translate -of GTIFF -projwin 101 42 108 35 earth_day_06m_p.tif day.tif

    # 绘制卫星图片，-G 表示在 DEM 上贴卫星图片
    gmt grdview @earth_relief_06m -R101/108/35/42/-1000/5000 -JM10c -JZ4c -N-1000+ggray -Gday.tif -Qi -Ba -Bz -BwsENZ -p60/25   
    # 绘制断层线
    gmt plot3d faults.xyz -W1p,brown -p
    # 绘制色标
    gmt colorbar -C -Ba -DJTC+o0/1c -p
    # 绘制底图（含指北针）
    gmt basemap -TdjLT+w1.5c+l+o1c -p60/25/2000

    rm day.tif faults.xyz earth_day_06m_p.tif
gmt end show

1.6地震剖线+MT图

earthquake.txt下载🔗

#!/bin/bash

gmt begin seismic-profile png
gmt set MAP_FRAME_TYPE=plain
# === (a) 绘制区域震中分布图和选定的测线 ===
gmt basemap -R100/110/20/30 -JM3i -Ba5f2.5 -BWeSn -Y7i
# 绘制所有地震位置
gmt plot earthquakes.dat -Sc0.1i -Gred -W0.25p
# 沿着测线每隔0.1度生成一个测点并绘制测线
gmt project -C100/22 -E106/30 -G0.1 | gmt plot -W2p,blue
echo '(a)' | gmt text -F+cBR+f12p,1 -Dj0.2c/0.2c

# === (b) 绘制选定测线0.2度范围内地震的震中分布图 ===
gmt basemap -R100/110/20/30 -JM3i -Ba5f2.5 -BWeSn -X4i
# 筛选出测线周围0.2度范围内的地震并绘制
gmt project earthquakes.dat -C100/22 -E106/30 -Fxy -W-0.2/0.2 | gmt plot -Sc0.1i -Gred -W0.25p
# 沿着测线每隔0.1度生成一个测点并绘制测线
gmt project -C100/22 -E106/30 -G0.1 | gmt plot -W2p,blue
echo '(b)' | gmt text -F+cBR+f12p,1 -Dj0.2c/0.2c

# === (c) 绘制选定测线0.2度范围内地震的经度/深度剖面图 ===
gmt basemap -R100/106/0/70 -JX7i/-1.0i -Bxa1+lLongitude -Bya20f10+l"Depth(km)" -BWeSn --MAP_TICK_LENGTH_PRIMARY=-0.03i -X-4i -Y-1.5i
# 筛选出测线周围0.2度范围内的地震并输出地震的经度和深度
gmt project earthquakes.dat -C100/22 -E106/30 -Fxz -W-0.2/0.2 | gmt plot -Sc0.1i -Gred -W0.25p
echo '(c)' | gmt text -F+cBR+f12p,1 -Dj0.2c/0.2c

# === (d) 绘制选定测线0.2度范围内地震的纬度/深度剖面图 ===
gmt basemap -R22/30/0/70 -JX7i/-1.0i  -Bxa1+lLatitude -Bya20f10+l"Depth(km)" -BWeSn --MAP_TICK_LENGTH_PRIMARY=-0.03i -Y-1.5i
# 筛选出测线周围0.2度范围内的地震并输出地震的纬度和深度
gmt project earthquakes.dat -C100/22 -E106/30 -Fyz -W-0.2/0.2 | gmt plot -Sc0.1i -Gred -W0.25p
echo '(d)' | gmt text -F+cBR+f12p,1 -Dj0.2c/0.2c
gmt end show

1.7南加州地震分布

catalog-lt2.0.txt下载🔗
station-17.txt下载🔗
station-3.txt下载🔗

gmt begin cat1.7 png
	gmt grdcut @earth_relief_01s -R-118/-117/35.3/36.3 -Gsc.grd
	gmt grdgradient sc.grd -A100 -Ne0.8 -Gtopo.grad
	gmt makecpt -Cgray -T-7000/7000/1000 -D
	gmt grdimage sc.grd -R-118/-117/35.3/36.3 -Itopo.grad -JM10c -BWSen -E100 -Ba0.2f0.1
	gawk '{if($5>=2.0 && $5<=2.5)print $3,$2}' catalog-lt2.0.txt|gmt plot -Sc0.1i -Gred -W0.05p -l"EQ>2.0"
	gawk '{if($5>=2.5 && $5<=3.0)print $3,$2}' catalog-lt2.0.txt |gmt plot -Sc0.07i -Gtan1 -W0.05p -l"EQ>2.5"
	gawk '{if($5>=3.0)print $3,$2}' catalog-lt2.0.txt |gmt plot -Sc0.06i -Gindianred1 -W0.05p -l"EQ>3.0"
	gawk '{print $1,$2}' station-17.txt |gmt plot -St0.250c -Gblue -W0.5p,blue -l"Broad band"
    gawk '{print $1,$2,$4}' station-17.txt |gmt text -F+f5p,1,black+jTL -D-0.35c/-0.18c -C0.035/0.035+to -Gwhite -Wwhite 
	gawk '{print $1,$2}' station-3.txt |gmt plot -St0.250c -Ggreen -W0.5p,green -l"Short period"
    gawk '{print $1,$2,$4}' station-3.txt |gmt text -F+f5p,1,black+jTL -D-0.35c/-0.18c -C0.035/0.035+to -Gwhite -Wwhite
gmt end show

1.8海岸线+比例尺+指北针

gmt begin coastline png
	##绘制海岸线和湖岸线
	#gmt coast -R-130/-50/20/60 -JM15c -Baf -W0.5p,black
	##1级海岸线黑色，2级海岸线浅红色，-A控制最小湖泊面积
	#gmt coast -R-130/-50/20/60 -JM15c -Baf -W1/0.5p,black -W2/1p,lightred -A5000
	## 填充-G陆地、-S水体与-C湖泊
	#gmt coast -R-130/-50/20/60 -JM15c -Baf -A5000 -Gred -Slightblue -Clightred
	##绘制-N1国界/-N2州界
	#gmt coast -R-130/-50/20/60 -JM15c -Baf -A5000 -Gred -Slightblue -Clightred -N1/1p -N2/0.25p
	##添加比例尺，将比例尺画在北纬 25° 西经 60° 处、制纬度为北纬 25° 处的比例尺，指北针同理
	##比例尺长度为 1000 千米、比例尺的风格为黑白相间的铁轨形式、显示比例尺对应的单位
	gmt coast -R-130/-50/20/60 -JM15c -Baf -A1000 -Gred -Slightblue -Clightred -Lg-60/25+c25+w1000k+f+u -Tdg-55/57+w0.5i+jCM
	echo -55 55 N |gmt text -F+f6p,4,black -Gwhite
gmt end show

1.9GPS场

参考网页1、网页2
gps_campagin.txt下载🔗
Indian.plate下载🔗
Philippine.plate下载🔗
China_tectonic.dat下载🔗
colombia.cpt下载🔗

#!/bin/bash
gmt begin gps jpg
gmt set FONT_ANNOT_PRIMARY 10 
gmt set MAP_FRAME_WIDTH 0.08 
gmt set MAP_TICK_LENGTH_PRIMARY 0.08 
R=65/140/10/60
J=D102.5/35/36/42/7.5i
gmt basemap -R$R -J$J -Bxa10f5 -Bya10f5 -BWSne
gmt grdcut @earth_relief_01m.grd -Gtopo.grd -R$R
#gmt grdgradient topo.grd -Ne0.7 -A50 -Gtopo_i.grd
gmt grdimage topo.grd -Ccolombia.cpt
gmt coast -R$R -J$J -A5000 -Dl -W0.3p -N1

gmt plot China_tectonic.dat -W1.2p,255,-
gmt plot Philippine.plate -W1.5p,lawngreen -Sf0.5c/0.1c+r+t+o0.1
gmt plot Indian.plate -W1.5p,lawngreen -Sf0.5c/0.1c+r+t+o0.1
awk '{print $1,$2,$3,$4,$5,$6,0}' gps_campagin.txt | gmt velo -Se0.02c/0.95/0 -A0.15c+e+p0.75p,red -Gred -W0.1p,red

echo 80 23 15 1c | gmt plot -SV0.5c+e -W4p -G0 
echo 133 22 310 1c | gmt plot -SV0.5c+e -W4p -G0 
echo 80 22 ~4 cm/yr | gmt text -F+f10p+a-40 -G255
echo 134 21.8 ~8 cm/yr | gmt text -F+f10p+a40 -G255

gmt end show

1.10中文字体

#!/usr/bin/env bash
gmt begin GMT_Chinese png,pdf
# 设置中文字体配置文件 cidfmap 的目录
gmt set PS_CONVERT="C-I/home/sasa/.gmt/"

# GMT 处理中文存在一些已知BUG
# 需要设置 PS_CHAR_ENCODING 为 Standard+ 以绕过这一BUG
gmt set PS_CHAR_ENCODING Standard+
gmt set FONT_TITLE 25p,41,black
gmt set FONT_LABEL 15p,39,black

gmt text -R0/8/0/4 -JX12c/4c -Bxaf+l"X轴" -Byaf+l"Y轴" -BWSen+t"中文标题" -F+f << EOF
2 3.5 25p,39,black 中文宋体
2 2.5 25p,40,blue  中文仿宋
2 1.5 25p,41,red   中文黑体
2 0.5 25p,42,green 中文楷体
4 3.5 25p,43,black 中文宋体
5 3.5 25p,44,blue  中文仿宋
6 3.5 25p,45,red   中文黑体
7 3.5 25p,46,green 中文楷体
EOF
gmt end show

1.11板块构造

plates.gmt 下载🔗

#!/usr/bin/env bash
gmt begin plates png
    gmt basemap -JN15c -Rg -Ba60
    gmt coast -Gwhite -S167/194/223
    # 从crust_type头段中提取出所有大陆地壳板块，并填充橘黄色，设置透明度80%
    gmt convert plates.gmt -S"crust_type=continental" | gmt plot -Gorange@80 -l"continental crust plates"
    # 从plate_type头段中提取出所有变形带板块，并填充红色，设置透明度50%
    gmt convert plates.gmt -S"plate_type=deformation zone" | gmt plot -W0p -Gred@50 -L -l"deformation zone"
    # 从poly_name头段中提取出青藏高原板块、华北板块，并分别填充蓝色和淡绿色
    gmt convert plates.gmt -S"poly_name=Tibet" | gmt plot -W0.5p -L -Gblue -l"Tibet plate"
    gmt convert plates.gmt -S"poly_name=North China" | gmt plot -W0.5p -Glightgreen -l"North China plate"
gmt end show

1.12构造+地形

global_gprv.gmt下载🔗

#!/usr/bin/env bash
gmt begin global_gprv png
    gmt basemap -JD105/35/36/42/10c -R70/140/3/60 -Baf
    gmt coast -Gwhite -S167/194/223

    # 提取出所有克拉通、地盾、被动板块边缘的边界数据，并填充颜色，设置30%透明度
    gmt convert global_gprv.gmt -S"prov_type=craton" | gmt plot -G168/171/210 -t30 -l"craton"
    gmt convert global_gprv.gmt -S"prov_type=shield" | gmt plot -G217/205/228 -t30 -l"shield"
    gmt convert global_gprv.gmt -S"prov_type=passive margin" | gmt plot -G196/233/236 -t30 -l"passive margin"

    # 提取出所有和阿尔卑斯-喜马拉雅造山事件相关的地质区域边界数据，并填充红色，设置30%透明度
    gmt convert global_gprv.gmt -S"lastorogen=Alpine-Himalayan" | gmt plot -Gred@30 -l"Alpine-Himalayan"

    # 提取出南海盆地的边界数据，绘制出黑色轮廓；提取出鄂尔多斯地块的边界数据，绘制出橘黄色轮廓
    gmt convert global_gprv.gmt -S"prov_name=South China Sea Basins" | gmt plot -W1.5p -l"South China Sea Basins"
    gmt convert global_gprv.gmt -S"prov_name=Ordos Block" | gmt plot -W1.5p,orange -l"Ordos Block"
gmt end show

1.13板块边界

boundaries.gmt下载🔗

#!/usr/bin/env bash

gmt begin boundaries png
    gmt basemap -JN15c -Rg -Ba60
    gmt coast -G244/243/239 -S167/194/223
    # 提取数据头段中不同类型的边界数据并绘图
    gmt convert boundaries.gmt -S"type=spreading center" | gmt plot -W0.5p,orange -l"spreading center"
    gmt convert boundaries.gmt -S"type=extension zone" | gmt plot -W0.5p,red -l"extension zone"
    gmt convert boundaries.gmt -S"type=subduction zone" | gmt plot -W0.5p,black -l"subduction zone"
    gmt convert boundaries.gmt -S"type=collision zone" | gmt plot -W0.5p,darkgreen -l"collision zone"
    gmt convert boundaries.gmt -S"type=dextral transform" | gmt plot -W0.5p,blue -l"dextral transform"
    gmt convert boundaries.gmt -S"type=sinistral transform" | gmt plot -W0.5p,pink -l"sinistral transform"
    gmt convert boundaries.gmt -S"type=inferred" | gmt plot -W0.5p,gray -l"inferred"

    # 设置图例属性
    gmt legend -DJMR+o0.1c -F+p1p+gwhite
gmt end show

1.14中国及邻区地质图

原图网址🔗
geo3al.gmt下载🔗
geoage.cpt下载🔗

#!/usr/bin/env bash
data=geo3al.gmt
cpt=geoage.cpt
gmt begin geo3al png
    gmt coast -R70/150/13/55 -JM22c -Gwhite -Baf -BWsNe
    gmt plot $data -C$cpt -aZ="GEN_GLG" -G+z
    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