使用GMT绘制地球化学稀土元素配分图

最近绘制地球化学图解，发现以前上学时用的Geokit已经收费咯...
其它相关的软件，如基于Excel VBA的Geoplot、GCDplot，基于R的GCDkit，基于CorelDraw的CGDK，以及最近基于matlab的matplot等等，都没咋用过，也没动力去研究了。
还是使用GMT吧，一招鲜。而且以后再绘制类似的图，稍微改一下就可以一键成图了。

# 稀土元素配分图解
#!/bin/bash
cat << EOF > xannots.txt
1 a La
2 a Ce
3 a Pr
4 a Nd
5 a Sm
6 a Eu
7 a Gd
8 a Tb
9 a Dy
10 a Ho
11 a Er
12 a Tm
13 a Yb
14 a Lu
15 a Y
EOF

# 球粒陨石标准值(S.-s. Sun and W. F. McDonough,1989)
La=0.237; Ce=0.612; Pr=0.095; Nd=0.467; Sm=0.153; Eu=0.058; Gd=0.2055;Tb=0.0374; Dy=0.254;
Ho=0.0566; Er=0.1655; Tm=0.0255; Yb=0.17; Lu=0.0254; Y=1.57 
# 上地壳标准化值，Rudnick and Gao(2003)
# La=31; Ce=63; Pr=7.1; Nd=27; Sm=4.7; Eu=1.0; Gd=4.0; Tb=0.7; Dy=3.9; Ho=0.83;
# Er=2.3; Tm=0.30; Yb=2.0; Lu=0.31; Y=21

# REE.csv数据格式
# La,Ce,Pr,Nd,Sm,Eu,Gd,Tb,Dy,Ho,Er,Tm,Yb,Lu,Y
# 53.4,101,12.9,46.2,7.36,1.73,5.49,0.54,1.62,0.22,0.53,0.05,0.23,0.079,6.07
data=REE.csv

# 开始绘图
gmt gmtset FONT_LABEL = 8p,39 # 坐标轴标签
gmt gmtset FONT_ANNOT_PRIMARY = 7p,4 # 刻度标注
gmt gmtset MAP_LABEL_OFFSET = 1p # 坐标轴标签到坐标轴的距离
gmt gmtset MAP_FRAME_PEN = 0.6p # 图框宽度
gmt begin 
gmt figure REE pdf,jpg A0.2c,E800
gmt basemap -R0/16/1/1000 -JX8c/5cl -Bxcxannots.txt -Bya1f3p+l"岩石/球粒陨石" -BWStr

# 数据处理
awk -F"[, ]+" 'NR>1{print $1/"'$La'", $2/"'$Ce'", $3/"'$Pr'", $4/"'$Nd'", $5/"'$Sm'", \
    $6/"'$Eu'", $7/"'$Gd'", $8/"'$Tb'", $9/"'$Dy'", $10/"'$Ho'", $11/"'$Er'", $12/"'$Tm'", \
    $13/"'$Yb'", $14/"'$Lu'", $15/"'$Y'"}' $data > tmp
# awk中NR为行数，NF为一行中的字段数
# γδPt1
awk 'NR>=1 && NR <=5{for(i=1;i<=NF;i++){print i, $i} print ">"}' tmp | gmt plot -W0.2p,red  -l@~'\147\144'@~Pt@-1@-+jRT+p0.2p
awk 'NR>=1 && NR <=5{for(i=1;i<=NF;i++){print i, $i} print ">"}' tmp | gmt plot -Sc0.1c -W0.2p,red -Gred  
# ηγβPt1
awk 'NR>=6 && NR <=9{for(i=1;i<=NF;i++){print i, $i} print ">"}' tmp | gmt plot -W0.2p,blue  -l@~'\150\147\142'@~Pt@-1@-
awk 'NR>=6 && NR <=9{for(i=1;i<=NF;i++){print i, $i} print ">"}' tmp | gmt plot -St0.1c -W0.2p,blue -Gblue
# ηγβmPt1
awk 'NR>=10 && NR <=12{for(i=1;i<=NF;i++){print i, $i} print ">"}' tmp | gmt plot  -W0.2p,green -l@~'\150\147\142'@~mPt@-1@-
awk 'NR>=10 && NR <=12{for(i=1;i<=NF;i++){print i, $i} print ">"}' tmp | gmt plot  -Ss0.1c -W0.2p,green -Ggreen
# ηγPt1
awk 'NR>=13 && NR <=15{for(i=1;i<=NF;i++){print i, $i} print ">"}' tmp | gmt plot  -W0.2p,cyan -l@~'\150\147'@~Pt@-1@-
awk 'NR>=13 && NR <=15{for(i=1;i<=NF;i++){print i, $i} print ">"}' tmp | gmt plot  -Ss0.1c -W0.2p,cyan -Gcyan
# δηοPt1
awk 'NR>=16 && NR <=18{for(i=1;i<=NF;i++){print i, $i} print ">"}' tmp | gmt plot  -W0.2p,brown -l@~'\144\150\157'@~Pt@-1@-
awk 'NR>=16 && NR <=18{for(i=1;i<=NF;i++){print i, $i} print ">"}' tmp | gmt plot  -Sa0.1c -W0.2p,brown -Gbrown

gmt end
rm xannots.txt gmt.* tmp