要计算某个点到海岸线的距离，思路很简单，先获取离散的海岸线数据，然后计算该点到所有海岸线数据点的距离，取其中的最小值即可认为是海岸线的距离。

获取海岸线数据

海岸线数据可以从 GMT 的 GSHHG 数据中提取，命令如下:

gmt pscoast -Di -Rd -W -M > coast.txt

简单解释一下:

  • GSHHG 包含了 5 种不同精度的数据，分别为 f 、 h 、 i 、 l 、 c，过高精度的数据包含了太多的细节、过低精度的数据丢失了太多细节。选择哪一种精度的数据需要仔细权衡。
  • -R 限定了感兴趣的海岸线数据的范围
  • -W 表示提取海岸线数据，还可以使用 -I 提取河流数据，使用 -N 提取国界数据
  • -M 表示将数据以 ASCII 格式输出（GMT4 里是 -m ）

该命令提取出的海岸线数据包含了全部 4 个 level 的数据，即海岸线数据、湖岸线数据、湖内岛数据以及湖内岛的内湖数据，如只需要 Level 1 即海岸线的数据，则需要使用 gmtconvert 对数据进行进一步处理:

gmt convert coast.txt -S"Level 1" > coast_1.txt

该命令的 -S"Level 1" 表示仅提头段记录中包含字符 “Level 1” 的线段。

在 GMT4 中，该命令应该这样写:

pscoast -Di -Rd -W -m > coast.txt
gmtconvert coast.txt -m -S"Level 1" > coast_1.txt

计算一点到海岸线的距离

计算一点到海岸线的距离，可以使用 mapproject 实现。例如要计算 (240, 30) 这个点到海岸线的距离，可以使用如下命令:

$ echo 240 30 | gmt mapproject -Lcoast_1.txt
240      30    183529.043998     -118.367666133   29.159609369

输出的 5 列数据为：要计算的点的经纬度、离海岸线的最小距离以及最小距离所对应的海岸线上的点。其中，最小距离的默认单位是 m 。也可以在 -Lcoast_1.txt 后加上不同的单位以指定最小距离的单位。比如:

$ echo 240 30 | gmt mapproject -Lcoast_1.txt/k
240      30    183.529043998     -118.367666133   29.159609369

的输出结果中第三列的单位就是 km。

完成了一个点的计算，多个点也没有什么难度。

计算网格点到海岸线的距离

要计算一个网格内所有点到海岸线的距离，也可以用上面的方法计算，当然有更简单的方式。

grdmath

GMT 的 grdmath 命令可以对网格数据进行多种数学操作，其中之一就是计算点到海岸线的距离。

对于网格的每个点，计算该点到海岸线数据的所有点的距离，找出离海岸线最小距离，将该最小距离作为该网格点的 Z 值。

grdmath 示例

$ gmt grdmath -R240/242/30/32 -I1 coast_1.txt LDIST -fg = dist.grd
$ gmt grd2xyz dist.grd
240 32  143.028015137
241 32  103.875091553
242 32  82.0870361328
240 31  249.026535034
241 31  206.842697144
242 31  138.674285889
240 30  183.548187256
241 30  111.494659424
242 30  92.4527893066

说明：

  • 此处通过 -R 和 -I 选项临时生成了一个网格文件，也可以直接指定一个网格文件
  • 参数格式为 操作数 操作符 = 输出网格文件名 ，操作符 LDIST 需要一个操作数 coast_1.txt 。注意这里 coast_1.txt 的位置必须在 LDIST 之前。
  • 默认情况下会计算网格点与线段之间的笛卡尔距离， -fg 选项表明输入数据为地理坐标，此时的输出为千米。
  • 通过 grd2xyz 可以看到计算的结果，这里计算出了 9 个点与海岸线的距离，单位千米。

附录

关于计算一点到海岸线的距离，前面使用了 mapproject 命令来实现。当初写博文的时候没有找到此方法，因而介绍了一个更复杂的方法。虽然本方法很复杂，没有必要再介绍，但考虑到其他有一些有意思的细节，故而保留，仅供参考。

上面的例子中， -R240/242/30/32 -I1 指定的网格中包含了 9 个点，然后计算出了 9 个点到海岸线的距离。当只需要计算某个点到海岸线的距离时，理论上只要网格中只包含一个点就可以了。

想要构建一个只含一个点的网格，网格必须使用像素网格配准方式。

使用默认的网格线配准方式，网格中至少会有四个网格点:

$ grdmath -R240/241/30/31 -I1 coast.txt LDIST -fg = dist.grd
$ grd2xyz dist.grd
240 31  249.291519165
241 31  207.351837158
240 30  183.434417725
241 30  111.675422668

通过 -F 使用像素配准方式，则生成的网格中只有一个点:

$ grdmath -R240/241/30/31 -I1 -F -fg coast_1.txt LDIST = dist.grd
$ grd2xyz dist.grd
240.5   30.5    184.780654907

想要计算点 (240.111, 30.222) 到海岸线的距离:

$ grdmath -R240.011/240.211/30.122/30.322 -I0.2 -F -fg coast_1.txt LDIST = dist.grd
$ grd2xyz dist.grd
240.111 30.222  188.54675293

注意理解 -R 和 -I 是如何选取的！

在 GMT5 下命令如下:

$ gmt grdmath -R240.011/240.211/30.122/30.322 -I0.2 -r -fg coast_1.txt LDIST = dist.grd
$ gmt grd2xyz dist.grd
240.111 30.222  188.537643433

参考

修订历史

  • 2014-02-01：初稿；
  • 2015-02-27：使用 gmtconvert 提取指定 Level 的数据； grdmath 中应使用 -fg 选项才会计算得到笛卡尔距离；Thanks to Jiangbo；
  • 2015-10-29：使用 mapproject 计算任意点到海岸线的距离；Thanks to Chuanxu Chen；