圣盖博山脉的难以到达极点

# 圣盖博山脉中的难达极点 来源：https://notes.secretsauce.net/notes/2015/05/06_poles-of-inaccessibility-in-the-san-gabriel-mountains.html 话说我和一位朋友外出徒步时，聊起了一个问题：圣盖博山脉（San Gabriel Mountains）中最难以到达的地方在哪里？这里“可到达”的定义是附近有道路或小径。我决定回答这个问题。相关代码放在一个新仓库中：https://github.com/dkogan/inaccessibility。 这被称为“难达极点”（Pole of Inaccessibility，http://en.wikipedia.org/wiki/Pole_of_inaccessibility）：即距离给定对象集合尽可能远的点。这类极点已知的有地球上最内陆的地点，或离陆地最远的地点。这里我们把范围限制在圣盖博山脉，并尽量远离道路和小径。 ## 方法 ### 输入数据处理 OpenStreetMap（http://www.openstreetmap.org/）拥有开放数据，可以用来绘制所有道路和小径。这就是输入数据集。 对于二维几何，计算难达极点最好的方法似乎是构建我们试图远离的几何对象的 Voronoi 图（http://en.wikipedia.org/wiki/Voronoi_diagram），然后找到对应最远点的 Voronoi 顶点。 我们的世界并不是二维的；事实上，地球是一个椭球体，上面还有起伏的地形。这里的主要目的是算出一个身强力壮的人能去的地方，并告诉所有人他们做到了，因此不需要极高的精度。所以我认为，假设世界是局部平坦的，并用基于 Voronoi 图的方法就足够了。于是，我构造一个最能描述查询区域的平面，把所有输入点投影到这个平面上。这个平面在查询区域中心与地球表面相切。如果要找像海洋那样大范围的难达极点，这显然行不通，但在这里没问题。 为了计算切平面，我假设地球是球形的。随着我沿着切平面远离切点，高程误差会增大： E = sqrt(R_earth² + d²) - R_earth 圣盖博山脉跨度约80公里，而切平面位于中间，所以最坏情况下 d = 40公里，误差约为125米。这已经足够好了。图示（源文件：https://notes.secretsauce.net/notes/2015/05/06_poles-of-inaccessibility-in-the-san-gabriel-mountains/plot_flat_earth_error.gp）： plot_flat_earth_error.svg 我完全忽略地球的椭球形状，也忽略地形；因为如果把地形纳入距离度量，就需要比生成 Voronoi 图更复杂的算法，而且会把“难以到达”这个概念弄得更模糊。 ### 难达极点的计算 我想使用最基础的 Voronoi 算法，所以输入只表示为点的集合，而不是线段。为了保证合理的精度，我确保每条道路至少每100米采样一次。 现在，在获得了足够密集的二维点集之后，我构造 Voronoi 图。如果不加约束，最远的点会跑到某一侧无穷远的地方，所以人们通常会把解限制在输入点的凸包内。这意味着难达极点要么位于某个 Voronoi 顶点，要么位于 Voronoi 边与输入凸包的交点处。在我的情形中，查询区域边缘的道路通常比内部多（山里的东西比平地少），所以我直接假设难达极点不在凸包上。这简化了实现，因为我只需要忽略查询区域之外的所有 Voronoi 顶点。 因此，我需要检查每个 Voronoi 顶点，计算它与相邻输入点之间的距离，并返回距离最大的那个顶点。 ## 实现 流程中的每一步都是一个独立的程序。这简化了实现，也方便分别处理每个部分。 ### 数据导入 首先查询 OSM。这是通过 `query.sh` 脚本完成的。它接收查询区域的角点