参考GMT中文社区绘制区域CORS站的站点分布和速度场图（脚本）

软件：GMT5.4 操作环境：win10虚拟机下的ubuntu16.04系统
在学习绘制区域CORS站的站点图以及速度场时，参考GMT中文社区里贴出来的一些脚本文件，结合自己的数据进行了一些修改。
1.绘制站点分布图
参考脚本网址（内含各省会的经纬度）：https://gmt-china.org/example/ex004/
在绘制自己所需的区域站点分布图时，需要自己制作该区域内各个站点的经纬度文件，我的文件是cors-station.dat，用GAMIT里面的制作l文件的方法，按照链接里省会的数据格式制作而成。

#!/bin/bash
#安徽省CORS站站点分布图10.18
PS=station11.ps
R=114/121/29/35
J=M18c

#准备底图所需要的地形数据
gmt grdcut earth_relief_01m.grd -R$R -GcutTopo.grd
gmt grdgradient cutTopo.grd -Ne0.7 -A50 -GcutTopo_i.grd
gmt grd2cpt cutTopo.grd -Cglobe -S-10000/10000/200 -Z -D > colorTopo.cpt

#绘制底图
gmt set MAP_GRID_PEN_PRIMARY 0.5p,dimgray,2_2:1
gmt set FORMAT_GEO_MAP ddd:mm:ssF MAP_FRAME_WIDTH 7p
gmt set FONT_ANNOT_PRIMARY 15p
gmt psbasemap -R$R -J$J -Bf2a2 -BWESN -Xc -Yc -K > $PS
gmt grdimage cutTopo.grd -IcutTopo_i.grd -R -J -CcolorTopo.cpt -Q -O -K >> $PS
#此处使用的中国国界线有误不得用于发表
gmt pscoast -R -J -ECN+p1p,black,- -W1/0.2p -I1/0.25p -O -K >> $PS
#绘制省界线
gmt psxy CN-border-La.dat -J$J -R$R -W1.8p,black -K -O >> $PS

#绘制colorbar
#gmt psscale -DjCB+w18c/0.3c+o0/-2.5c+h -R -J -CcolorTopo.cpt -BWSEN -Bxa2000f400+l"Elevation/m" -G-8000/8000 -O -K >> $PS

#绘制各个站点 用大圆套小圆的方式标注 
gmt psxy cors-station.dat -J$J -R$R -Sc0.16c -Gblack -K -O >> $PS
gmt psxy cors-station.dat -J$J -R$R -Sc0.27c -W0.6p,red4 -K -O >> $PS

#AHCORS各个站点的名称
gmt pstext cors-station.dat -J$J -R$R -F+f18p,1,black+j -Dj0.15c/0.15c -K -O >> $PS

#绘制周围省名称
gmt pstext province.dat -J$J -R$R -Dj0.15c/0.15c  -F+f23p,17,gray28+j  -K -O >> $PS


gmt psxy -R -J -T -O >> $PS
gmt psconvert -A0.5c -Tg -P $PS
rm -rf $PS
rm gmt.* cutTopo*.grd colorTopo.cpt

2.绘制速度场
参考脚本地址（内含数据）：https://gmt-china.org/example/ex015/
在绘制这幅图的时候需要注意：下边绘制的图例需要一直做修改；速度文件需要自己制作，格式就按照链接里提供的数据，以此作为参考，我制作了gps_ahcors.txt文件作为输入文件。

#!/bin/bash
#绘制AHCORS速度场
PS=velo_feild12.ps
#全国的
#R=70/135/15/55
#AHCORS站的
R=114/121/29/35
J=M18c

#准备底图所需要的地形数据
gmt grdcut earth_relief_01m.grd -R$R -GcutTopo.grd
gmt grdgradient cutTopo.grd -Ne0.7 -A50 -GcutTopo_i.grd
gmt grd2cpt cutTopo.grd -Cglobe -S-10000/10000/200 -Z -D > colorTopo.cpt

#绘制底图
gmt set MAP_GRID_PEN_PRIMARY 0.5p,dimgray,2_2:1
gmt set FORMAT_GEO_MAP ddd:mm:ssF MAP_FRAME_WIDTH 8p
gmt set FONT_ANNOT_PRIMARY 18p
#-Bf2a2前后两个2分别表示图边的经纬度间隔 -BWESN大小写表示所绘的图上下左右是否有标注，大写表示有
gmt psbasemap -R$R -J$J -Bf2a2 -BWESN -Xc -Yc -K > $PS
gmt grdimage cutTopo.grd -IcutTopo_i.grd -R -J -CcolorTopo.cpt -Q -O -K >> $PS

#此处使用的中国国界线有误  不可用于发表
gmt pscoast -R -J -ECN+p1p,black,- -W1/0.2p -I1/0.25p -O -K >> $PS
gmt psxy CN-border-La.dat -J$J -R$R -W1.8p,black -K -O >> $PS

#绘制colorbar
#gmt psscale -DjCB+w18c/0.3c+o0/0.5c+h -R -J -CcolorTopo.cpt -BWSEN -Bxa2000f400+l"Elevation/m" -G-8000/8000 -O -K >> $PS

#绘制GPS速度场
#-Se后面三参数表示：速度值为1的矢量的长度/置信度0.95/文本的大小
#-W控制矢量以及误差椭圆的轮廓的宽度，颜色，线型
#-G矢量填充色
#-A控制矢量的属性，0.15c是矢量头的大小，+e表示在矢量尾端绘制箭头，+p0.75p矢量线段部分的宽度
#awk '{print $1,$2,$3,$4,$5,$6,0}' gps_campagin.txt | gmt psvelo -J -R -Se0.05c/0.95/0 -A0.15c+e+p0.75p -Gblue -W0.2p,blue -K -O >> $PS

#绘制矢量箭头 -A0.65c+e+p1.6p（065表示箭头大小）（1.6表示矢量线的粗细）-W1.0p,brown4 （1.0表示箭头前的误差椭圆的粗细 后面是颜色）
awk '{print $1,$2,$3,$4,$5,$6,0}' gps_ahcors.txt | gmt psvelo -J -R -Se0.05c/0.95/0 -A0.65c+e+p1.6p -Gbrown4 -W1.0p,brown4 -K -O >> $PS

#绘制各个站点 用大圆套小圆的方式标注 
gmt psxy cors-station.dat -J$J -R$R -Sc0.1c -Gblack -K -O >> $PS
#gmt psxy cors-station.dat -J$J -R$R -Sc0.27c -W0.6p,red4 -K -O >> $PS


#绘制周围省名称
gmt pstext province.dat -J$J -R$R -Dj0.15c/0.15c  -F+f23p,17,gray28+j  -K -O >> $PS

#psvelo很难画图例，以下代码展示了如何绘制图例，需要进行大量微调
#这个是通过在图上标注四个点来确定图里的位置，该区域内底为白色
gmt psxy -J -R -Gwhite -K -O >> $PS << EOF
118.5 34.4
121 34.4
121 35
118.5 35
EOF
#此处用来显示“20±1mm” 前面是坐标（图里的经纬度）
gmt pstext -J -R -F+f10p+jML -K -O >> $PS << EOF
118.6 34.8 20\2611 mm/yr
EOF

#这里用来画图例上的线段、箭头和箭头前的误差椭圆（我这里设置的和前面一样）
gmt psvelo -J -R -Se0.05c/0.95/0 -A0.65c+e+p1.6p -Gbrown4 -W1.0p,brown4 -K -O >> $PS << EOF
118.7 34.6 20 0 1 1 0
EOF

#这里是文本标注，调整得到的
gmt pstext -J -R -F+f10p+jML -M -K -O >> $PS << EOF
> 119.4 34.8 0.25 2c c
Velocity of 
EOF

gmt pstext -J -R -F+f10p+jML -M -K -O >> $PS << EOF
> 120.1 34.8 0.25 2c c
AHCORS 
EOF

gmt pstext -J -R -F+f10p+jML -M -K -O >> $PS << EOF
> 119.6 34.6 0.25 2c c
2013.0-2018.5
EOF



#gmt psvelo -J -R -Se0.05c/0.95/0 -A0.15c+e+p0.75p -Gred -W0.2p,red -K -O >> $PS << EOF
#123 16 20 0 1 1 0
#EOF

#gmt pstext -J -R -F+f6p+jML -M -K -O >> $PS << EOF
#> 120.5 27 0.25 2c c
#Continuous Station
#1998-2014
#EOF

gmt psxy -R -J -T -O >> $PS
gmt psconvert -A -P -Tg $PS
rm -rf $PS
rm gmt.* cutTopo*.grd colorTopo.cpt