【注】部分GMT绘图示例参考于GMT中文社区及博客,并基于Linux的gmt6版本成功运行绘图脚本,用于记录学习gmt的部分启蒙&有用例子
目录
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全国地形图+震级分布
#!/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文件制作等高线+测线高程
#!/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图
#!/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