搜索地球表面上的点（工程）

发表于 | 更新于
为了找到自己的位置，我们需要将全世界编码

这里尝试解决一个需求：寻找附近的 POI（Points of Interest，相关地点，下同）。一些数据库，比如 MangoDB，ElasticSearch 已经有现成的方案。但如果你用的数据库，比如 AWS DyanmoDB，Google Cloud Firestore，没有已经封装好的 API，可以按照本文的方式，通过地理索引，实现附近的 POI 查询。如果不想看关于地理索引的分析，可以直接跳到「代码实现」一节。

第一种方法，也是最直接、最简单的办法是：

  1. 确定中心点
  2. 计算数据库中各个点和中心点的距离
  3. 保留合适的点

当数据很少的时候，这样可能是最好的方法，但是如果 POI 数量增加，m 次查询的计算量会按照 O(m * n) 的程度增长。对于大数据量是不可行的。

上面的方法可以改良。如果这些点按照经纬度排序，那么可以通过坐标先排除一部分明显很远的点。于是，一个可能的方法是：

  1. 确定中心点
  2. 只查询附近的经纬度，取回满足初步条件的点（可称为候选点）
  3. 计算候选点到中心点的距离
  4. 保留合适的点

但这个方法也面临问题。

首先是基于经纬度建立地理索引的技术问题。我们考虑下面几个方案和他们的问题：

a. 建立经、纬度联合索引。但是，地理查询需要用到经纬度的上界和下界，在 SQL 中使用这种范围查询会导致联合索引失效。对于分布式 NoSQL，联合索引或许可以通过分区键（partition key）和排序键（sort key）来实现，但是如果以小数点后好几位的经纬度划分分区键，会导致分区过多；分区过多，会导致查询次数过多，增加数据库负担。

b. 把经、纬度 2 个值归结为一个索引值。生成这种索引的方法有很多，简单的字符串拼接、Geohash1等，都可以把经纬度两个值合并成一个值，而且连续的索引值代表的区域在地理上也大致接近。但由于地球是个球体，相同的经度差（∆lng）越靠近两极，就越只有很短距离；最极端的例子是在北纬或南纬 90 度的极点，从 0 度经线到 180 度经线的最短距离是零。也就是说，如果用同样的经纬度差值来划分地理区域，在海南能划出一大片区域，而在黑龙江只能划出很小的区域。这就导致，要想查询一个覆盖范围，比如周围 10km，那么，在海南地区只需要查询一个索引下的区域，而在黑龙江，则会击中很多个索引。这对于 SQL 数据库的影响是需要设计球面几何算法，来计算需要查询的地理索引的范围；对于分布式 NoSQL 的影响是，查询普遍击中一个分区会造成热键，或者需要查询过多分区时则会增加数据库负担。因此都需要平衡。

c. 找到一个让地理分区的面积大小基本接近的索引算法，解决 b 中关于球面几何计算的问题，也能平衡数据分区键下的数据量大小。s2 geometry 就是让同一级别的地理编码所代表的面积相近2，即使在南北极，一块分区的面积大小，也不会跟赤道附近的差太多。同样的算法还有 uber 开发的 h3。

本文就尝试使用 s2 geometry 来实现地理索引，解决查询附近的 POI 这个需求。

另外，一个小话题，如果你是看了 AWS 的博客或者 Google 搜索排名很靠前的一篇文章，而想在自己的 dynamodb 中用 Geohash 做数据分区，那么他们可能误导了你。这些博客中提到的代码库其实用的是 s2 geometry 做的地理编码。可以通过其 github 上的代码确认。

总之，到这里，我们的解决方案是：

  1. 利用合适的地理编码系统（s2 geometry cell）为 POI 的地理位置生成索引；同时，在实现过程中，让一个索引对应的地理区域大小要合适，避免热键和分区过多的问题
  2. 查询索引，找到初步候选点
  3. 计算候选点到中心点的距离
  4. 保留合适的点

代码实现

下面进入工程实现的部分（语言为 golang）。

  1. 生成中心点（Point 类型）
  2. 生成目标搜索范围：从中心点向外画一个圆盘（Cap 类型）
  3. 计算覆盖区域：目标范围覆盖了哪些地理分区（Cells）
  4. 计算数据分区：地理分区的编码就是数据分区的分区键
  5. 请求数据库，查询候选点
  6. 计算取回的候选点是否在目标区域内

第一步，生成中心点 center。这个中心点的左右可以由前端获取后传递到后端。

1
2
3
4
5
6
7
import "github.com/golang/geo/s2"
import "github.com/golang/geo/s1"
import "fmt"

func main() {
	latLng := s2.LatLngFromDegrees(39.916667, 116.383333) // Beijing
	center := s2.PointFromLatLng(latLng)

第二步，生成搜索范围，用 Cap 类型表示。Cap 的意思是帽子，如此命名可能是因为一个圆盘范围扣在地球上，就像一个帽子戴在头上。

By Wikipedia

1
2
3
const EarthRadius = 6378.137
cap := center.CapBound()
region := cap.Expanded(s1.Angle(distance / EarthRadius)) // distance is set by you

因为 s2 用角度大小表示距离3，所以需要先把我们所用的距离，比如 10 公里，转化为地球弧度角。在地球表面，10 公里（在过球心的圆面上）所形成的弧度（radian）大概是 10 / 6378.137，其中 6378.137 是地球在赤道上的半径距离（km）4。

第三步，计算搜索范围 region 覆盖了哪些地理分区（cells），因为数据分区是根据地理分区制定的，只有得到了地理分区，我们才能知道候选点存在哪个数据分区。

图片使用 s2.sidewalklabs.com 制作

1
2
3
rc := &s2.RegionCoverer{MinLevel: 5, MaxLevel: 20, MaxCells: 8}
covering := rc.Covering(region)
covering.Normalize()

其中，

  1. MaxLevel 代表最小的 cell，即覆盖区域精细到什么程度，如果数值太大，可能会返回一些细小的 cell；
  2. MinLevel 代表最大的 Cell，即不能使用更大的 cell 来代表覆盖区域，MinLevel 也不能小于数据分区使用的编码级别，因为级别太小的话（覆盖范围太大）就分不清具体指向哪个数据分区。
  3. Covering 方法返回的结果是一个 CellUnion 类型的值（包含一个或多个 cell id），而这个结果可能没有归一化（normalize），所以有时 4 个 cell id 已经构成一个完整的 cell（更大的 cell），但却没有被合并5，增加了接下来数据查询的复杂度，因此需要进一步归一化，即调用 Normalize 方法。

第四步，计算数据分区。上面得到的 cell ids 并不直接就是数据分区的分区键（partition key），因为有些 cell id 的范围比数据分区所代表的地理范围小。但通过找到 cell id 的父级，就能得到其所对应的数据分区。

1
pk := v.Parent(parentLevel) // parentLevel is set by you

如果有 cell ids 所对应的数据分区相同，则可以合并到同一个分区中搜索数据。并且，利用 sort key 可以只取回此数据分区中相关的候选点，类似于 cellID.RangeMin() <= sort_key <= cellID.RangeMax()，因为一个 cell 中的 cell ids 是连续的（编码连续，但数值间隔 2），所以在 min 和 max 之间的点肯定属于这个 cell。

第五步，数据查询。利用上面得到的数据分区键和排序键（sort key），取回一定数量的候选点。

最后，再次计算候选点是否处在我们预想的范围 region 内：

1
2
3
4
5
for _, v := range pointsYouRetrievedFromDatabase {
  if region.ContainsPoint(v) {
    // do something ...
  }
}

可能的问题

  1. 人口扎堆在大城市，对应的数据分区更容易被击中，成为热键（如果有很多人用的话）
  2. 如果开发者想在前端对经纬度进行脱敏，可以使用 js 的 s2 geometry 开发包，生成一个模糊的地理区域再传入后端。
  3. 需要平衡计算覆盖区域、数据查询、再次确认距离（即，第四、五、六步）之间的时间开销，在某一步降低消耗可能会导致剩下的步骤产生更多消耗。比如，计算覆盖区域时不使用精确计算，可能导致数据查询时返回更多候选点，进而需要分页处理以及更多次确认距离。监控、平衡各步的消耗，降低总体消耗，可能是个合理的选择。
1. Wikipedia Contributors. (2020, March 19). Geohash. Retrieved April 10, 2020, from Wikipedia website: https://en.wikipedia.org/wiki/Geohash#Non-linearity‌ ↩
2. S2 Cells. (2020). Retrieved April 10, 2020, from S2Geometry website: https://s2geometry.io/devguide/s2cell_hierarchy#coordinate-systems‌ ↩
3. Quick Start. (2020). Retrieved April 11, 2020, from S2Geometry website: https://s2geometry.io/devguide/cpp/quickstart#s2closestpointquery‌ ↩
4. Earth Fact Sheet. (2020). Retrieved April 11, 2020, from Nasa.gov website: https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html‌ ↩
5. S2 Cells. (2020). Retrieved April 12, 2020, from S2Geometry website: https://s2geometry.io/devguide/s2cell_hierarchy#approximating-regions‌ ↩