生成测线上的高程并绘制地形

 1 w=90
 2 e=130
 3 n=30
 4 s=28
 5 start=($w $n)
 6 end=($e $s)
 7 #采样点为200个，精度为15弧分，生成数据有两列：distance(km)和elevation(m)
 8 echo ${end[*]} | awk ‘{
 9     r[NR]=$1;p[NR]=$2}END{
10         for (i=0; i<201; i++)
11         print r[1]+i*(r[2]-r[1])/200,p[1]+i*(p[2]-p[1])/200
12     }‘ | gmt grdtrack -G@earth_relief_15m_p | gmt project -Fpz -Q -S -C$w/$n -E$e/$s > topo_dh.tmp
13 #画图时为了方便需要将数据区域闭合，生成数据可以将海拔正负分段，每段开头为"> >"
14 cat topo_dh.tmp | awk ‘{l=$1;h=$2}{
15     if (NR==1) print "> >\n"$1,0"\n"$0
16     else if (hh*h<0) print ll,0"\n> >\n"l,0"\n"l,h
17     else if (ll*l>0) print l,h
18 }{ll=l;hh=h}END{print $1,0}‘ > topo_plot.tmp
19 #画图
20 x=$(sort -n -k1 topo_dh.tmp)
21 length=$(echo "scale=1; ${length}*5/8" | bc)
22 gmt plot topo_plot.tmp -JX${x}c/5c -R0/$x/-4000/4000 -BWS -Byg5000+l"Elevation (m)" -Bxa500+l"Distance (km)" -W0.1p,black -Ggray85 -L -png topo_AB