GMT的海岸线空间数据分析功能

下面介绍：使用GMT和海岸线数据绘任意距离的海岸延伸线（如12海里线等）；使用海岸线距离进行站点数据的筛选（如提取沿海50km以内的GNSS点位）。

根据离岸的距离对GNSS站点进行提取

在利用GNSS水汽对卫星高度计进行湿延迟评估时，需要提取距离海岸线小于50km的站点做比对。这里我们使用GMT快速分析提取数据并绘图展示。

绘图采用的数据主要是中国地震局陆态网、IGS和自然资源部沿海业务化GNSS系统，卫星数据主要是我国的海洋二号以及国外的Jason，绘图数据见文末链接。

中国陆态网数据和国际IGS数据为共享开放资料。由于政策原因，沿海业务化GNSS资料暂未公开。海洋二号高度计资料和国外高度计资料已全面开放共享。

绘图任务主要步骤总结为：

• grdmath的LDISTG计算网格点到GSHHG岸线的距离，生成网格文件；
• grdlandmask分离陆海；
• grdcontour提取某值的等值线；
• select依据点位到某等值线的距离筛选站点。

Code

下面是windows下的代码：

REM             GMT EXAMPLE xxxx
REM             2019-04-14 leiyang@fio.org.cn
REM Purpose:    Using grdmath to show the distance from GSHHG  coastline data
REM
gmt gmtset FORMAT_GEO_MAP = dddF MAP_FRAME_WIDTH=2p
gmt gmtset FONT_ANNOT_PRIMARY 7p,Helvetica,black FONT_LABEL 7p,Helvetica,black 

set ps=wet.ps

gmt grdmath -R100/140/15/45 -A1000000/0/4 -Dc -I10m LDISTG = dist_to_gshhg.nc
gmt grdlandmask -R -Dl -I10m -N1/-1 -Gland_mask.nc -V
gmt grdmath dist_to_gshhg.nc land_mask.nc MUL = file.nc
gmt grdsample dist_to_gshhg.nc -R100/140/15/45 -I1m -Gfile2.nc
REM Outut 50km coordinates.
gmt grdcontour file2.nc -JM4.5i -R100/140/15/45 -C50, -D > out.d50
gmt makecpt -Ccopper -T-0/50/1> t.cpt

gmt grdimage file2.nc -JM4.5i -R100/130/15/45 -P -K -Ct.cpt -t80> %ps%
gmt psxy -R -J out.d50 -W0.4p,green -O -K >>%ps%
gmt pscoast -R -JM4.5i -Dl -A10000/0/1 -B5g5 --FONT_TITLE=10p -K -Swhite -W0.1p --MAP_ANNOT_OBLIQUE=45 -O -N1>> %ps%

REM Begin the sites selection.
gmt pscoast -R -M -W -Dc -A10000> coastal.txt
gmt select trc.dat -Lcoastal.txt+d50k -fg > subset.d2

REM satellite tracks
gmt psxy -R -J ja.dat -W1p,black -O -K -t80>>%ps%
gmt psxy -R -J hy2.dat -W1p,red -O -K -t80>>%ps%

REM All CMONOC sites
gmt psxy -R -J trc.dat -Sc0.1 -Gblack -O -K -t70 >>%ps%

REM tide sites
gawk "{print $3, $4, 4, 0, 1, 1, $1}" yh_tide.txt2  >temp.d
REM Coastal GNSS sites
gawk "{print $1, $2, 4, 0, 1, 1, $4}" gnss.txt >temp.d01
gmt psxy temp.d01  -R  -J  -K -O  -Gblack -St0.05i -t70 >>%ps%

REM Plot the selected sites 
gawk "{print $1, $2, 4, 0, 1, 1, $4}" gnss2.txt | gmt psxy -R -J -O -Gred -Sc0.08i -K>>%ps%
gmt psxy -R -J subset.d2 -Ss0.15 -GDEEPSKYBLUE1 -O -K>>%ps%
gmt psxy -R -J gnss3.txt -Ss0.2 -Gred -O -K>>%ps%
gmt psxy -R -J gnss3.txt -Ss0.12 -Gblue -O -K >>%ps%

echo 1.2116450e+02,2.4953600e+01 | gmt psxy -R -J -Sa0.2 -Gblue -O -K>>%ps%
echo 118.38,24.46  | gmt psxy -R -J -Sa0.2 -Gblue -O -K>>%ps%

gawk "!/SDYT/ && !/SDQD/ {print $1, $2, $3}" gnss3.txt | gmt pstext  -R -J -F+f7p,black+jTL -O -K -Gwhite -D0.2/0.1>> %ps%
gawk "/SDYT/ || /SDQD/ {print $1, $2, $3}" gnss3.txt | gmt pstext  -R -J -F+f7p,black+jTL -O -K -Gwhite -D-0.5/0.4>> %ps%
gawk "{print $1, $2, $4}" gnss2.txt | gmt pstext -R -J -F+f7p,black+jTL -O -K -Gwhite -D0.2/0.1>> %ps%
echo 1.2116450e+02,2.4953600e+01,TWTF | gmt pstext -R -J -F+f7p,black+jTL -O -K -Gwhite -D0.2/0.1>> %ps%

REM add legend
echo S 0.1i c 0.05c black 1p,black 0.3i Total CMONOC >le.d
echo S 0.1i s 0.15c DEEPSKYBLUE1 .7p,DEEPSKYBLUE1 0.3i Coastal CMONOC >>le.d
echo S 0.1i s 0.2c blue 0.7p,red 0.3i Selected CMONOC >>le.d
echo S 0.1i t 0.05i gray - 0.3i Costal CGN>>le.d
echo S 0.1i c 0.08i  red - 0.3i Selected CGN >>le.d
echo S 0.1i a 0.15c blue 0.7p,blue 0.3i IGS >>le.d
gmt pslegend -D+w1.6i+jBL+o2.8i/0.2i -R -J -O -F+p1p+gwhite --FONT_ANNOT_PRIMARY=7p le.d >> %ps%

REM Creat kml file to be loaded in Google earth.
gmt 2kml subset.d2 -Gred+f -Fs > mypoints.kml

gmt psconvert %ps%  -A -P -Tf

200公里线绘制

使用GMT岸线绘制领海线、专属经济区等政治界线比较危险，因为上面的数据分析使用了低分辨率的岸线，并忽略了很多海岛。这里选择绘制距离大陆200km的等值线。

具体做法是对grdcontour的-C做修改，设置为-C200,。其他距离可以通过设定-C后面的数值实现。

Code

gmt grdcontour file2.nc -JM4.5i -R100/140/15/45 -C200, -D > out.d50
gmt makecpt -Ccopper -T-0/200/10> t.cpt

gmt grdimage file2.nc -JM4.5i -R100/130/15/45 -P -K -Ct.cpt > %ps%
gmt psxy -R -J out.d50 -W0.4p,green -O -K >>%ps%
gmt pscoast -R -JM4.5i -Dl -A10000/0/1 -B5g5 --FONT_TITLE=10p -Gwhite -W0.1p --MAP_ANNOT_OBLIQUE=45 -O -N1>> %ps%

结果

距离岸线50km的GNSS站点分布

距离岸线200km的等值线

坐标文件

陆态网GNSS站点坐标trc.dat：https://www.jianguoyun.com/p/DZNDa0kQ2PCQBxih5ogE
HY-2高度计卫星轨迹hy2.dat：https://www.jianguoyun.com/p/DZvezRQQ2PCQBxim5ogE
Jason-2高度计轨迹ja.dat：https://www.jianguoyun.com/p/Dfulo24Q2PCQBxiq5ogE
沿海验潮站概略坐标yh_tide.txt2：https://www.jianguoyun.com/p/DaYSJ8kQ2PCQBxis5ogE
沿海GNSS业务化系统概略坐标gnss.txt：https://www.jianguoyun.com/p/DWhRCEUQ2PCQBxiv5ogE
绘图用的个别GNSS站点：https://www.jianguoyun.com/p/DdUqCG4Q2PCQBxi05ogE； https://www.jianguoyun.com/p/DcBuNBsQ2PCQBxi15ogE

注意事项

gmt grdmath -R100/140/15/45 -A1000000/0/4 -Dc -I10m LDISTG = dist_to_gshhg.nc

这一行代码中的-Dc表示低分辨率海岸线，如果计算的区域比较大，使用低分辨率岸线可以加快时间。不然，计算速度可能非常慢。