搜索地球表面上的点(工程)

说明:

  1. 这里讲不了太多原理,因为我也还不懂
  2. 地球覆盖范围计算基于 Google S2 Geometry 软件库
  3. 数据存储、查询基于 DynamoDB,迁移到其他 NoSQL 应该差不多

这里尝试解决一个需求:寻找附近的 POI(Points of Interest,相关地点,下同)。

第一种方法,也是最直接、最简单的办法是:

  1. 确定中心点
  2. 计算数据库中各个点和中心点的距离
  3. 保留距离在要求范围内(假设为 [0, D])的点

当数据很少的时候,这样可能是最好的方法,但是如果数据增加,计算量会按照 O(n) 函数线性增加。对于大数据量是不可行的。

上面的方法可以改良。如果这些点按照经纬度排序,那么可以通过坐标先排除一部分明显很远的点。如果纬度或经度的存储使用 B+ tree 等 index,定位某个范围的时间可能更短。于是,一个可能的方法是:

  1. 搜索临近经纬度,只取回满足初步条件的点(可称为候选点)
  2. 计算候选点到中心点的距离
  3. 保留距离在 [0, D] 内的点

但这个方法也面临问题。

首先,DynamoDB 等分布式数据库是按分区存储数据的,如果通过划分经纬度来制定数据分区键(partition key),那么有些目标区域可能跨越 N 个数据分区,需要做 N 次查询请求;分区越细,目标区域跨越的分区就越多,N 就越大。如果把所有数据都存在一个分区内,可以利用纬度和经度的组合字段做 index(类似于关系型数据库中的多字段 index),但又容易造成此分区过于繁忙(热键问题)。当然,某些情况下这些处理方法是可以接受的。

其次,虽然相同的纬度差(∆lat)代表的距离总是大致相等,但是相同的经度差(∆lng)却可能代表非常不同的距离。相同的 ∆lng 在赤道附近有更长的距离,越靠近南北极越短,因此需要调整经度区间的大小,才能保持搜索范围的大小和靠近赤道的区域一致。

对于上述之第一个问题,可以把经、纬度 2 个值编码为一个值,最简单的方法就是拼接,比如,把 N30° - N31° 之间且 E110° - E111° 之间的区域编码为 N30.5-E110.5。这样,此范围内的点都可以对应相同的地理编码,进而存在同一个数据分区内(以地理编码做分区键)。如果某个中心点之外,半径 D 之内没有超过 N30.5-E110.5 编码所表示的区域,那么仅在一个数据分区内就能取回所有的候选点,只做一次请求。而且,数据分区内还可以利用 sort key 等方法进一步缩小候选点的范围。

一个地理编码表示的范围越大,你就越有可能用 1 个请求取回所有的候选点。但如果编码覆盖的区域太大,很多查询可能同时击中一个分区,同样热键问题。这里需要看怎么平衡两种消耗(请求次数,热键的可能性)。

上面说的地理编码方式很简单。另外也有很多编码规范供你选择,Geohash 就是一种,应该算是比较知名的。

然而 Geohash 并没有解决上述之第二个问题:单位经度差所代表的实际距离随着纬度升高而大幅减小。所以,同一级别的地理编码,在赤道代表的面积比在南北极代表的面积大很多1。虽然在平面地图上看起来南北极附近的编码区域好像很大,但那是因为平面地图把极点附近聚拢的区域展开成了跟赤道一样宽。(类似地,俄罗斯、南极洲都没有平面地图上显示的比例那么大;根据 thetruesize.com 计算,如果把中国挪到俄罗斯的纬度位置,在平面地图上看起来会比在中纬度大得多)

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

s2 geometry 让同一级别的地理编码所代表的面积相近2,即使在南北极也不会和赤道的差距过于悬殊。所以,到这里,我们基本上可以解决上述提到的问题:

  1. 利用地理区域筛选候选点,解决逐一计算距离则计算量过大的问题
  2. 利用合适的地理编码系统(s2 geometry cell)解决地理区域划分的部分难点
  3. 调整地理编码所代表的区域大小,使请求次数和产生热键的概率都保持在可接受的范围内

至于为什么一定要用 s2 geometry cell,其实不一定,但我只了解过 s2。其实还有其他的地理编码系统,比如 uber 开发的 h3,但我还不是很了解。

代码实现

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

  1. 生成中心点(Cell 类型)
  2. 扩大中心点成一个圆盘范围(Cap 类型)
  3. 计算覆盖区域:圆盘范围覆盖了哪些地理分区(Cells)
  4. 计算数据分区:地理分区的编码就是数据分区的分区键
  5. 请求数据库,查询候选点
  6. 计算取回的候选点是否在目标区域内

作为事先准备,我们先生成一个 cell id。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
import "github.com/golang/geo/s2"
import "github.com/golang/geo/s1"
import "fmt"

func main() {
latLng := s2.LatLngFromDegrees(39.916667, 116.383333) // Beijing
cellID := s2.CellIDFromLatLng(latLng)
cellIDUint := uint64(cellID)
base2 := strconv.FormatUint(cellIDUint, 2)
base10 := strconv.FormatUint(cellIDUint, 10)
base16 := strconv.FormatUint(cellIDUint, 16)
fmt.Printf("cell id:\nbase2(%s)\nbase10(%s)\nbase16(%s)\ntoken(%s)\n",
base2, base10, base16, cellID.ToToken(),
)
// cell id:
// base2(11010111110000010100101111001001110010010000011111100101111001)
// base10(3886697679673227641)
// base16(35f052f27241f979)
// token(35f052f27241f979)
}

因此,我们用作实验的 cell id 是:35f052f27241f979

搜索工程就从接收这个 cell id 开始。至于为什么要传入 cell id,而不是经纬度。因为:

  1. cell id 也可以用来表示地图上的点,和传入经纬度坐标是一样的。
  2. cell id 有现成的方法可以控制精度,不需要精确定位时,可以使用低级别的 cell id。

当然,如果你想使用精确的经纬度来生成点(S2Point)也是可以的。

第一步,生成中心点,我用 Cell 类型的值表示。

1
2
3
// generate a cell given a token of cell id
cellID := "35f052f27241f979"
cell := s2.CellFromCellID(s2.CellIDFromToken(cellID))

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

By Wikipedia

1
2
3
const EarthRadius = 6378.137
cap := cell.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 的半径,比 distance 大一点,因为 cell 本身也有一定宽度,如果不能接受这样的误差,可以主动调整 distance,抵消 cell 的宽度。

第三步,计算搜索范围 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. cell id 的发送方(比如网页前端)和接收方(比如服务器后端)对于搜索范围的确定可能有偏差,因为前端标示距离的算法和后端使用的 s2 geometry 计算距离的算法可能不一样;而且使用较模糊的 cell id(level 较低的 cell),本来就难以保证前后端中心点一致

  3. 需要平衡计算覆盖区域、数据查询、再次确认距离(即,第四、五、六步)之间的时间开销,在某一步降低消耗可能会导致剩下的步骤产生更多消耗。比如,计算覆盖区域时不使用精确计算,可能导致数据查询时返回更多候选点,进而需要分页处理以及更多次确认距离。监控、平衡各步的消耗,降低总体消耗,可能是个合理的选择。

最后,s2 geometry 连接各地理分区所用的空间填充曲线很有意思,让人怀疑这个世界到底是什么样子。

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‌