0%

GMT 计算任意点到海岸线的最短距离

发表于 更新于 阅读次数： Valine：

问：我想计算海洋上的任意点到海岸线的（最短）距离，并制作一个网格文件存储这些距离，请问如何？

答：使用grdmath可获得全球任意点到海岸线距离的网格。这里的距离就是垂直最短距离。

问：我想绘制一条线，表示到海岸线的距离。

答：通过设置grdcontour可以导出指定的等值线坐标到文本，比如我们要绘制距离海岸线50km的线，或者专属经济区的线等。

问：我想根据到海岸线的距离提取点位，可以吗？

答：可以。下面的例子中有这个功能的使用。如从我国GNSS陆态网站点中提取出距离海岸线50km的点，并绘制。

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
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 

echo GMT EXAMPLE xxxx
set ps=example_xx.ps

gmt grdmath -R100/170/00/55 -A0/0/4 -Dc -I30m LDISTG = dist_to_gshhg.nc
gmt grdlandmask -R -Dc -I30m -N1/-1 -Gland_mask.nc -V
gmt grdmath dist_to_gshhg.nc land_mask.nc MUL = file.nc
gmt makecpt -Crainbow -T-100/2500/20> t.cpt
REM gmt grdimage dist_to_gshhg.nc -JM4.5i -R110/170/20/55 -P -K -Ct.cpt -X0.75i -Y2i > %ps%
REM gmt grdcontour dist_to_gshhg.nc  -JM4.5i -R -C50 -A10+f10p,Helvetica,black -L00/20000 -GL150/45/165/35 -Wa0.75p,black -Wc0.25p,black -K > %ps%
gmt grdcontour file.nc  -JM4.5i -R -C50, -A50,+f10p,Helvetica,black -Wa0.75p,black -Wc0.25p,black -K > %ps%
gmt pscoast -R -J -O -Df -EJP,CN,TW+gred@60 -ECN.37+gyellow@60 -A10000/0/1 -B5g5 -B+t"Distances from GSHHG crude coastlines" --FONT_TITLE=10p -K -Glightyellow -Slightblue --MAP_ANNOT_OBLIQUE=45>> %ps%
gmt grdcontour file.nc -JM4.5i -R -C50, -D > out.d50
gmt select trc.dat -Lout.d50+d100k  -fg > subset.d

REM 使用GSHHG 提取海岸线，需要使用GSHHG的bin文件，此方法可以提取出闭合的全球地理数据 
rem 进入 gshhg bin文件夹，执行：
REM gmt gshhg gshhs_i.b --IO_SEGMENT_MARKER=N > gshhs_i.txt # http://gmt.soest.hawaii.edu/boards/1/topics/6324 , http://www.soest.hawaii.edu/pwessel/gshhg/
REM 使用pscoast 提取海岸线
gmt pscoast -R100/170/00/55 -M -W -Dc -A10000> coastal.txt
REM another way to select data
gmt select trc.dat -Lcoastal.txt+d50k -fg > subset.d2

gmt psxy -R -J trc.dat -Sc0.05 -Gblack -O -K >>%ps%
gmt psxy -R -J subset.d -Sc0.2 -Gblue -O -K>>%ps%
gmt psxy -R -J subset.d2 -Sc0.1 -Gred -O -K>>%ps%
gmt psxy -R -J out.d50 -W0.4p,green -O >>%ps%
REM gmt psxy -R -J coastal.txt -W0.4p,green -O >>%ps%
gmt psconvert example_xx.ps -A -P -Tg

结果

获取数据

sites data trc.dat

注意点

  • grdmath的网格设置如果太密，范围太大，则计算速度较慢。
  • gshhg 提取岸线，需进入数据文件夹。