空间索引之 Google S2

在空间索引类问题中，一个最普遍而又最重要的问题是：给定你某个点的坐标，你如何能够在海量的数据点中找到他所在的区域以及最靠近他的点？”，比方说客户在路上突然想吃饭了，那么就要根据他的位置查询最近的餐馆并做出推荐。

通常情况下，一提到查找类问题，我们就会想到二分查找或者是B树查找。但是问题在于我们不仅要找到这个点，而且要找到这个点附近的点。因此对于以经纬度来确定的坐标又不好直接进行二分查找。通常情况下我们会用R树、KD树或者是四叉树之类的数据结构来存储这些点从而高效的做到临近点的查找。但是这些数据结构通常都会存在数据冗余，以及不稳定的查改效率；况且抛开他们的时间效率、空间效率以及算法复杂度不谈，用了这些数据结构也就意味这我们放弃了使用现成强大的数据库而自己编写数据查改系统，这显然是繁琐而又没有必要的过程。而且当系统需要扩展到分布式计算的时候就更不如使用那些分布式的数据库了。

这时我们就会想到，如果能够把一个二维的信息转化为一维的数据加以存储，那么我们不就可以直接存储到数据库中做到快速的查找了么？

GeoHash 算法

GeoHash是由Gustavo Niemeyer提出的，目的原本是为地球上的每一个点（根据经纬度）确定一条短的URL作为唯一标识。只是后来被广泛的应用到空间检索方面。GeoHash所做的事就是把一个坐标点映射到一个字符串上，每一个字符串代表的就是一个以经纬度划分的矩形区域。Geohash是一种分级的数据结构，把空间划分为网格。Geohash 属于空间填充曲线中的 Z 阶曲线（Z-order curve）的实际应用。何为 Z 阶曲线？

上图就是 Z 阶曲线。这个曲线比较简单，生成它也比较容易，只需要把每个 Z 首尾相连即可。Geohash算法的理论基础就是基于 Z 曲线的生成原理。Z 阶曲线同样可以扩展到三维空间。

Geohash 能够提供任意精度的分段级别。一般分级从 1-12 级。

我们可以利用 Geohash 的字符串长短来决定要划分区域的大小。这个对应关系可以参考上面表格里面 cell 的宽和高。一旦选定 cell 的宽和高，那么 Geohash 字符串的长度就确定下来了。这样我们就把地图分成了一个个的矩形区域了。关于上表中的Bounding Boxes随着迭代过程的变化，可直观的体现在下图中：

地图上虽然把区域划分好了，但是还有一个问题没有解决，那就是如何快速的查找一个点附近邻近的点和区域呢？Geohash 有一个和 Z 阶曲线相关的性质，那就是一个点附近的地方(但不绝对) hash 字符串总是有公共前缀，并且公共前缀的长度越长，这两个点距离越近。

GeoHash的编码流程

对一个地理坐标编码时，按照初始区间范围纬度[-90,90]和经度[-180,180]，计算目标经度和纬度分别落在左区间还是右区间。落在左区间则取0，右区间则取1。然后，对上一步得到的区间继续按照此方法对半查找，得到下一位二进制编码。当编码长度达到业务的进度需求后，根据“偶数位放经度，奇数位放纬度”的规则，将得到的二进制编码穿插组合，得到一个新的二进制串。最后，根据base32的对照表，将二进制串翻译成字符串，即得到地理坐标对应的目标GeoHash字符串。

以坐标“30.280245,120.027162”为例，计算其GeoHash字符串。首先对纬度做二进制编码：

将[-90,90]平分为2部分，“30.280245”落在右区间(0,90]，则第一位取1。
将(0,90]平分为2分，“30.280245”落在左区间(0,45]，则第二位取0。

不断重复以上步骤，得到的目标区间会越来越小，区间的两个端点也越来越逼近“30.280245”。

下图的流程详细地描述了前几次迭代的过程：

按照上面的流程，继续往下迭代，直到编码位数达到我们业务对精度的需求为止。完整的15位二进制编码迭代表格如下：

得到的纬度二进制编码为10101 01100 01000。

按照同样的流程，对经度做二进制编码，具体迭代详情如下：

得到的经度二进制编码为11010 10101 01101。

按照“偶数位放经度，奇数位放纬度”的规则，将经纬度的二进制编码穿插，得到完成的二进制编码为：11100 11001 10011 10010 00111 00010。由于后续要使用的是base32编码，每5个二进制数对应一个32进制数，所以这里将每5个二进制位转换成十进制位，得到28,25,19,18,7,2。对照base32编码表，得到对应的编码为：wtmk72。

GeoHash字符串的长度与精度的对应关系如下：

Geohash 和 Z 阶曲线的关系

回顾最后一步合并经纬度字符串的规则，“偶数位放经度，奇数位放纬度”这个规则就是 Z 阶曲线。

下图显示整数坐标（x轴就是纬度，y轴就是经度）$0 \leq x \leq 7,0 \leq y \leq 7$的二维情况下的 z 值(以十进制和二进制表示)。交错的二进制坐标值生成二进制 z 值，如下所示。将 z 值按其数字顺序连接起来就得到了递归的 z 形曲线。二维 z 值也称为四键值。

Geohash 的优缺点

Geohash的优点

纬度和经度是二维的数据，建索引的时候需要对纬度和经度同时建索引，Geohash 能将二维的数据转化为一维的数据，可以用B树等索引
Geohash 生成的字符串代表的不是地图上的一个点，而是地图上一个矩形的区域，在一定程度上能保证隐私
字符串越长，表示的范围能够更加的精确。 8位的geohash 编码的经度能达到19米左右，9位的geohash经度能达到米级，一般情况下能满足我们大部分的需求
GeoHash值的前缀相同的位数越多，代表的位置越接近，可以方便索引。（反之不成立，位置接近的GeoHash值不一定相似）

Geohash的缺点

Geohash 的缺点之一也来自 Z 阶曲线。Z 阶曲线有一个比较严重的问题，虽然有局部保序性，但是它也有突变性。在每个 Z 字母的拐角，都有可能出现顺序的突变。

空间填充曲线

在数学分析中，有这样一个难题：能否用一条无限长的线，穿过任意维度空间里面的所有点？1890年，意大利数学家Giuseppe Peano发明能填满一个正方形的曲线，叫做皮亚诺曲线，其构造方法如下：取一个正方形并且把它分出九个相等的小正方形，然后从左下角的正方形开始至右上角的正方形结束，依次把小正方形的中心用线段连接起来；下一步把每个小正方形分成九个相等的正方形，然后上述方式把其中中心连接起来，将这种操作手续无限进行下去，最终得到的极限情况的曲线就被称作皮亚诺曲线。

一年后，即1891年，希尔伯特就作出了这条曲线，叫希尔伯特曲线（Hilbert curve）。

希尔伯特曲线的特点：

降维。作为空间填充曲线，希尔伯特曲线可以对多维空间有效的降维。
稳定。当n阶希尔伯特曲线，n趋于无穷大的时候，曲线上的点的位置基本上趋于稳定。
连续。希尔伯特曲线是连续的，所以能保证一定可以填满空间。

之后还有很多变种的空间填充曲线，龙曲线(Dragon curve)、高斯帕曲线(Gosper curve)、Koch曲线(Koch curve)、摩尔定律曲线(Moore curve)、谢尔宾斯基曲线(Sierpiński curve)、奥斯古德曲线(Osgood curve)。

参考链接：https://en.wikipedia.org/wiki/Space-filling_curve

Google S2算法

S2其实是来自几何数学中的一个数学符号 S²，它表示的是单位球。S2 这个库其实是被设计用来解决球面上各种几何问题的。接下来就看看怎么用 S2 来解决多维空间点索引的问题的。

球面坐标转换

球面上的一个点，在直角坐标系中，可以这样表示：

x = r * sin θ * cos φ
y = r * sin θ * sin φ 
z = r * cos θ

通常地球上的点我们会用经纬度来表示：

基于地球上任意的一个经纬度的点S(lat, lng)，都可以转换成 f(x,y,z)。

球面变平面

S2 是如何把球面碾成平面的？首先在地球外面套了一个外切的正方体，如下图。

从球心向外切正方体6个面分别投影。S2 是把球面上所有的点都投影到外切正方体的6个面上。

上图左边的是投影到正方体一个面的示意图，实际上影响到的球面是右边那张图。从侧面看，其中一个球面投影到正方体其中一个面上。我们可以在球的6个方向上，把45°的辅助圆画出来，见下图左边。

上图左边的图画了6个辅助线，蓝线是前后一对，红线是左右一对，绿线是上下一对。这样我们就可以把投影到外切正方体6个面上的球面画出来，见上图右边。投影到正方体以后，我们就可以把这个正方体展开。一个正方体的展开方式有很多种。不管怎么展开，最小单元都是一个正方形。在 Google S2 中，它是把地球展开成如下的样子：

S2是用的正方形。这样第一步的球面坐标进一步的被转换成 f(x,y,z) -> g(face,u,v)，face是正方形的六个面，u，v对应的是六个面中的一个面上的x，y坐标。

球面矩形投影修正

上一步我们把球面上的球面矩形投影到正方形的某个面上，形成的形状类似于矩形，但是由于球面上角度的不同，最终会导致即使是投影到同一个面上，每个矩形的面积也不大相同。

上图就表示出了球面上个一个球面矩形投影到正方形一个面上的情况。

经过实际计算发现，最大的面积和最小的面积相差5.2倍。见上图左边。相同的弧度区间，在不同的纬度上投影到正方形上的面积不同。现在就需要修正各个投影出来形状的面积。如何选取合适的映射修正函数就成了关键。目标是能达到上图右边的样子，让各个矩形的面积尽量相同。

线性变换是最快的变换，但是变换比最小。tan() 变换可以使每个投影以后的矩形的面积更加一致，最大和最小的矩形比例仅仅只差0.414。可以说非常接近了。但是 tan() 函数的调用时间非常长。如果把所有点都按照这种方式计算的话，性能将会降低3倍。最后谷歌选择的是二次变换，这是一个近似切线的投影曲线。它的计算速度远远快于tan()，大概是tan()计算的3倍速度。生成的投影以后的矩形大小也类似。不过最大的矩形和最小的矩形相比依旧有2.082的比率。

至此，球面上的点S(lat,lng) -> f(x,y,z) -> g(face,u,v) -> h(face,s,t)。目前总共转换了4步，球面经纬度坐标转换成球面xyz坐标，再转换成外切正方体投影面上的坐标，最后变换成修正后的坐标。

点与坐标轴点相互转换

在 S2 算法中，默认划分 Cell 的等级是30，也就是说把一个正方形划分为 2^30 * 2^30个小的正方形。那么上一步的s,t映射到这个正方形上面来。到这一步，是h(face,s,t) -> H(face,i,j)。

坐标轴点与希尔伯特曲线 Cell ID 相互转换

最后一步，如何把 i，j 和希尔伯特曲线上的点关联起来呢？

在变换之前，先来解释一下定义的一些变量。posToIJ 代表的是一个矩阵，里面记录了一些单元希尔伯特曲线的位置信息。把 posToIJ 数组里面的信息用图表示出来，如下图：

同理，把 ijToPos 数组里面的信息用图表示出来，如下图：

posToOrientation 数组里面装了4个数字，分别是1,0,0,3。lookupIJ 和 lookupPos 分别是两个容量为1024的数组。这里面分别对应的就是希尔伯特曲线 ID 转换成坐标轴 IJ 的转换表，和坐标轴 IJ 转换成希尔伯特曲线 ID 的转换表。

上图是一个4阶希尔伯特曲线。初始化的实际过程就是初始化4阶希尔伯特上的1024个点的坐标与坐标轴上的x，y轴的对应关系表。

举个例子，下表是i，j在递归过程中产生的中间过程。下表是lookupPos 表计算过程：

至此，整个球面坐标的坐标映射就已经完成了。球面上的点S(lat,lng) -> f(x,y,z) -> g(face,u,v) -> h(face,s,t) -> H(face,i,j) -> CellID。目前总共转换了6步，球面经纬度坐标转换成球面xyz坐标，再转换成外切正方体投影面上的坐标，最后变换成修正后的坐标，再坐标系变换，映射到 [0,2^30^-1]区间，最后一步就是把坐标系上的点都映射到希尔伯特曲线上。

S2 Cell ID 数据结构

S2 Cell ID 数据结构直接关系到不同 Level 对应精度的问题。

上图左图中对应的是 Level 30 的情况，右图对应的是 Level 24 的情况。

在 S2 中，每个 CellID 是由64位的组成的。可以用一个 uint64 存储。开头的3位表示正方体6个面中的一个，取值范围[0,5]。3位可以表示0-7，但是6，7是无效值。64位的最后一位是1，这一位是特意留出来的。用来快速查找中间有多少位。从末尾最后一位向前查找，找到第一个不为0的位置，即找到第一个1。这一位的前一位到开头的第4位（因为前3位被占用）都是可用数字。绿色格子有多少个就能表示划分多少格。不同 level 可以代表的网格的面积：

S2 与 Geohash 对比

Geohash 有12级，从5000km 到7cm。中间每一级的变化比较大。有时候可能选择上一级会大很多，选择下一级又会小一些。这种情况选择多长的 Geohash 字符串就比较难选。选择不好，每次判断可能就还需要取出周围的8个格子再次进行判断。S2 有30级，从 0.7cm² 到 85,000,000km² 。中间每一级的变化都比较平缓，接近于4次方的曲线。所以选择精度不会出现 Geohash 选择困难的问题。
Geohash 需要 12 bytes 存储。S2 的存储只需要一个 uint64 即可存下。
S2 库里面不仅仅有地理编码，还有其他很多几何计算相关的库。地理编码只是其中的一小部分。本文没有介绍到的 S2 的实现还有很多很多，各种向量计算，面积计算，多边形覆盖，距离问题，球面球体上的问题，它都有实现。

Python环境下使用Google S2

import s2sphere


def get_cellid_from_latlng(lat, lng, level=20):
    ll = s2sphere.LatLng.from_degrees(lat, lng)
    cell = s2sphere.CellId().from_lat_lng(ll)
    return cell.parent(level).to_token()


def lat_lng_to_cell_id(lat, lng, level=10):
    region_cover = s2sphere.RegionCoverer()
    region_cover.min_level = level
    region_cover.max_level = level
    region_cover.max_cells = 1
    p1 = s2sphere.LatLng.from_degrees(lat, lng)
    p2 = s2sphere.LatLng.from_degrees(lat, lng)
    covering = region_cover.get_covering(
        s2sphere.LatLngRect.from_point_pair(p1, p2))
    # we will only get our desired cell ;)
    return covering[0].id()


def middle_of_cell(cell_id):
    cell = s2sphere.CellId(cell_id)
    lat_lng = cell.to_lat_lng()
    return lat_lng.lat().degrees, lat_lng.lng().degrees


def coords_of_cell(cell_id):
    cell = s2sphere.Cell(s2sphere.CellId(int(cell_id)))
    coords = []
    for v in range(0, 4):
        vertex = s2sphere.LatLng.from_point(cell.get_vertex(v))
        coords.append([vertex.lat().degrees, vertex.lng().degrees])
    return coords


def get_position_from_cell(cell_id):
    cell = s2sphere.CellId(id_=int(cell_id)).to_lat_lng()
    return s2sphere.math.degrees(cell.lat().radians), s2sphere.math.degrees(cell.lng().radians), 0


if __name__ == "__main__":
    print(get_cellid_from_latlng(39.993379, 116.513977, 10))
    print(get_cellid_from_latlng(39.993379, 116.513977, 11))
    print(lat_lng_to_cell_id(39.993379, 116.513977, 20))
    print(middle_of_cell(lat_lng_to_cell_id(39.993379, 116.513977, 20)))
    print(coords_of_cell(lat_lng_to_cell_id(39.993379, 116.513977, 20)))
    print(get_position_from_cell(lat_lng_to_cell_id(39.993379, 116.513977, 20)))

参考链接：

空间索引之 Google S2