高效的空间索引算法——Geohash 和 Google S2_python geohash

1 前言

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

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

 这时我们就会想到，如果能够把一个二维的信息转化为一维的数据加以存储，那么我们不就可以直接存储到数据库中做到快速的查找了么？本篇文章就来介绍2种比较通用的空间点索引算法：Geohash和Google S2。

2 Geohash算法

2.1 Geohash理论介绍

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

 上图就是 Z 阶曲线。这个曲线比较简单，生成它也比较容易，只需要把每个 Z 首尾相连即可。Z 阶曲线同样可以扩展到三维空间。只要 Z 形状足够小并且足够密，也能填满整个三维空间。

Geohash算法的理论基础就是基于 Z 曲线的生成原理。Geohash 能够提供任意精度的分段级别。一般分级从 1-12 级。

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

 地图上虽然把区域划分好了，但是还有一个问题没有解决，那就是如何快速的查找一个点附近邻近的点和区域呢？Geohash 有一个和 Z 阶曲线相关的性质，那就是一个点附近的地方(但不绝对) hash 字符串总是有公共前缀，并且公共前缀的长度越长，这两个点距离越近。由于这个特性，Geohash 就常常被用来作为唯一标识符。用在数据库里面可用 Geohash 来表示一个点。Geohash 这个公共前缀的特性就可以用来快速的进行邻近点的搜索。越接近的点通常和目标点的 Geohash 字符串公共前缀越长（但是这不一定，也有特殊情况，下面举例会说明）。

 Geohash 有几种编码形式，常见的有2种：base 32 和 base 36。

• base 32 的版本用32个字符，“0123456789bcdefghjkmnpqrstuvwxyz”，编码对照表如下。

• base 36 的版本对大小写敏感，用了36个字符，“23456789bBCdDFgGhHjJKlLMnNPqQrRtTVWX”，编码对照表如下。

2.2 Geohash的编码流程

举例说明一

 以 base-32为例，对一个地理坐标编码时，按照初始区间范围纬度[-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。通过查表，我们选取字符串长度为6的矩形来网格化这张地图。经过查询，美罗城的经纬度是[31.1932993, 121.43960190000007]。

 先处理纬度。地球的纬度区间是[-90,90]。把这个区间分为2部分，即[-90,0)，[0,90]。31.1932993位于(0,90]区间，即右区间，标记为1。然后继续把(0,90]区间二分，分为[0,45)，[45,90]，31.1932993位于[0,45)区间，即左区间，标记为0。一直划分下去。

 再处理经度，一样的处理方式。地球经度区间是[-180,180]

 纬度产生的二进制是101011000101110，经度产生的二进制是110101100101101，按照“偶数位放经度，奇数位放纬度”的规则，重新组合经度和纬度的二进制串，生成新的：111001100111100000110011110110，最后一步就是把这个最终的字符串转换成字符，对应需要查找 base-32 的表。11100 11001 11100 00011 00111 10110 转换成十进制是 28 25 28 3 7 22，查表编码得到最终结果：wtw37q。

 我们还可以把这个网格周围8个各自都计算出来，然后可视化到地图上。可以看出，这邻近的9个格子，前缀都完全一致。都是wtw37。

 如果我们把字符串再增加一位，会有什么样的结果呢？Geohash 增加到7位。从下图可见，当Geohash 增加到7位的时候，网格更小了，美罗城的 Geohash 变成了 wtw37qt。

 把6位和7位都组合到一张图上面来看。可以看到中间大格子的 Geohash 的值是 wtw37q，那么它里面的所有小格子前缀都是 wtw37q。可以想象，当 Geohash 字符串长度为5的时候，Geohash 肯定就为 wtw37 了。

精度的问题

 GeoHash字符串的长度与精度的对应关系如下（数据一部分来自 Wikipedia）：

2.3 Geohash 和 Z 阶曲线的关系

 回顾上述编码流程最后一步合并经纬度二进制字符串的规则，“偶数位放经度，奇数位放纬度” 这个规则其实就是 Z 阶曲线。

 下图显示整数坐标（x轴就是纬度，y轴就是经度）$0\le x\le 7,0\le y\le 7$ 的二维情况下的 z 值(以十进制和二进制表示)。 交错的二进制坐标值生成二进制 z 值，如下所示。

 将 z 值按其数字顺序连接起来就得到了递归的 z 形曲线。 二维 z 值也称为四键值（下图只是演示，实际上数字是从0开始的）。

2.4 Geohash 的优缺点

Geohash的优点

 Geohash 的优点很明显，它利用 Z 阶曲线进行编码。而 Z 阶曲线可以将二维或者多维空间里的所有点都转换成一维曲线。在数学上成为分形维。并且 Z 阶曲线还具有局部保序性。Z 阶曲线通过交织点的坐标值的二进制表示来简单地计算多维度中的点的z值。一旦将数据被加到该排序中，任何一维数据结构，例如二叉搜索树，B树，跳跃表或（具有低有效位被截断）哈希表 都可以用来处理数据。通过 Z 阶曲线所得到的顺序可以等同地被描述为从四叉树的深度优先遍历得到的顺序。这也是 Geohash 的另外一个优点，搜索查找邻近点比较快。

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

Geohash的缺点

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

 看上图中标注出来的蓝色的点点。每两个点虽然是相邻的，但是距离相隔很远。看右下角的图，两个数值邻近红色的点两者距离几乎达到了整个正方形的边长。两个数值邻近绿色的点也达到了正方形的一半的长度。

 Geohash 的另外一个缺点是，如果选择不好合适的网格大小，判断邻近点可能会比较麻烦。
 看上图，如果选择 Geohash 字符串为6的话，就是蓝色的大格子。红星是美罗城，紫色的圆点是搜索出来的目标点。如果用 Geohash 算法查询的话，距离比较近的可能是wtw37p、wtw37r、wtw37w、wtw37m。但是其实距离最近的点就在 wtw37q。如果选择这么大的网格，就需要再查找周围的8个格子。

 如果选择 Geohash 字符串为7的话，那变成黄色的小格子。这样距离红星星最近的点就只有一个了。就是 wtw37qw。

 所以说，如果网格大小、精度选择的不好，那么查询最近点还需要再次查询周围8个点。

2.5 Python环境下使用Geohash

 geohash的最大用途就是附近地址搜索了。不过，从geohash的编码算法中可以看出它的一个缺点：位于格子边界两侧的两点，虽然十分接近，但编码会完全不同。实际应用中，可以同时搜索当前格子周围的8个格子，即可解决这个问题。

 使用python来实现Geohash的库已经有多个，这里整理记录几个常用的库：

• python-geohash

  项目地址：

  https://github.com/hkwi/python-geohash
  python-geohash · PyPI

  官方没有说明文档，以下为使用示例：

  import geohash
  
  print(geohash.encode(30.723514, 104.123245, precision=5)) # wm6nc
  
  print(geohash.decode("wm6nc")) # (30.73974609375, 104.12841796875) #中心经纬度
  
  print(geohash.decode("wm6nc", delta=True)) # (30.73974609375, 104.12841796875, 0.02197265625, 0.02197265625) #中心及偏差
  
  print(geohash.decode_exactly("wm6nc")) # (30.73974609375, 104.12841796875, 0.02197265625, 0.02197265625)
  
  print(geohash.bbox("wm6nc")) # {'s': 30.7177734375, 'w': 104.1064453125, 'n': 30.76171875, 'e': 104.150390625}
  
  print(geohash.neighbors("wm6nc")) # ['wm6nb', 'wm6nf', 'wm6n9', 'wm6n8', 'wm6nd', 'wm6p1', 'wm6p0', 'wm6p4']
  
  print(geohash.expand("wm6nc")) # ['wm6nb', 'wm6nf', 'wm6n9', 'wm6n8', 'wm6nd', 'wm6p1', 'wm6p0', 'wm6p4', 'wm6nc']
  1
  2
  3
  4
  5
  6
  7
  8
  9
  10
  11
  12
  13
  14
  15

• Geohash

  项目地址：

  https://github.com/vinsci/geohash
  Geohash · PyPI

  使用示例：

  import Geohash
  
  print(Geohash.encode(30.723514, 104.123245, precision=5)) # wm6nc
  
  print(Geohash.decode("wm6nc")) # (30.73974609375, 104.12841796875)
  
  print(Geohash.decode("wm6nc", delta=True)) # (30.73974609375, 104.12841796875, 0.02197265625, 0.02197265625)
  
  print(Geohash.decode_exactly("wm6nc")) # (30.73974609375, 104.12841796875, 0.02197265625, 0.02197265625)
  1
  2
  3
  4
  5
  6
  7
  8
  9

  注意：安装后导出可能会报错，解决办法可以参考 这篇文章

• geohash2

  项目地址：

  https://github.com/dbarthe/geohash/
  geohash2 · PyPI

  使用示例：

  import geohash2
  
  print(geohash2.encode(30.723514, 104.123245, precision=5)) # wm6nc
  
  print(geohash2.decode("wm6nc")) # ('30.7', '104.1')
  
  print(geohash2.decode_exactly("wm6nc")) # (30.73974609375, 104.12841796875, 0.02197265625, 0.02197265625)
  
  print(geohash2.encode(30.7, 104.1, precision=5)) # wm6n8
  1
  2
  3
  4
  5
  6
  7
  8
  9

  注意：从上面的数据可以看到decode后获取的不是中心点的坐标，而是对坐标进行了精度的调整，导致上述问题，所以不推荐使用。

• geolib

  项目地址：

  https://github.com/joyanujoy/geolib
  geolib · PyPI

  使用示例：

  from geolib import geohash
  
  print(geohash.encode(30.723514, 104.123245, precision=5)) # wm6nc
  
  print(geohash.decode("wm6nc")) # Point(lat=Decimal('30.73974609375'), lon=Decimal('104.12841796875'))
  
  print(geohash.bounds("wm6nc")) # Bounds(sw=SouthWest(lat=30.7177734375, lon=104.1064453125), ne=NorthEast(lat=30.76171875, lon=104.150390625))
  
  print(geohash.neighbours("wm6nc")) # Neighbours(n='wm6p1', ne='wm6p4', e='wm6nf', se='wm6nd', s='wm6n9', sw='wm6n8', w='wm6nb', nw='wm6p0')
  
  print(geohash.adjacent("wm6nc", 'n')) # wm6p1
  1
  2
  3
  4
  5
  6
  7
  8
  9
  10
  11

• ProximityHash

  项目地址：

  https://github.com/ashwin711/proximityhash
  proximityhash · PyPI

  简介：给定中心坐标和半径，ProximityHash生成一组覆盖圆形区域的地理HASH。它还可以使用GeoRaptor创建不同级别的geohash的最佳组合来表示圆，从最高级别开始迭代，直到绘制出最佳混合。结果精度与初始geohash级别的精度相同，但数据大小大大减小，从而提高了速度和性能。

• georaptor

  项目地址：

  https://github.com/ashwin711/georaptor
  georaptor · PyPI

  简介：GeoRaptor通过从最高级别开始迭代直到绘制出最佳混合，创建跨不同级别的geohash的最佳组合来表示多边形。结果精度与初始geohash级别相同，但对于大型多边形，数据大小会显著减小，从而提高速度和性能。

• TransBigData

  项目地址：

  https://github.com/ni1o1/transbigdata
  transbigdata· PyPI
  TransBigData – for Transportation Spatio-Temporal Big Data

  使用示例：

  # 输入经纬度与精度，输出geohash编码
  transbigdata.geohash_encode(lon, lat, precision=12)
  # 输入geohash编码，输出经纬度
  transbigdata.geohash_decode(geohash)
  # 输入geohash编码，输出geohash网格的地理信息图形Series列（geohash编码生成栅格矢量图形）
  transbigdata.geohash_togrid(geohash)
  1
  2
  3
  4
  5
  6

  简介：TransBigData是为交通时空大数据处理和分析而开发的Python包。TransBigData为出租车GPS数据、共享单车GPS数据、公交GPS数据等常见交通时空大数据提供快速、简洁的处理方法。它包括栅格化、数据质量分析、数据预处理、数据集计数、轨迹分析、GIS处理、地图底图加载、坐标和距离计算、数据可视化等一般方法。

• 其他项目

  https://github.com/mathieuripert/geoh
  https://github.com/Bonsanto/polygon-geohasher
  https://github.com/wmorin/geohashit
  https://github.com/derrickpelletier/geohash-poly
  https://github.com/tammoippen/geohash-hilbert
  https://github.com/kylebebak/py-geohash-any

• GEOHASH的可视化项目

  https://github.com/missinglink/leaflet-spatial-prefix-tree
  https://github.com/rzanato/geohashgrid
  Geohash encoding/decoding
  https://github.com/soorajb/geohash-grid

3 空间填充曲线和分形

 在介绍第二种多维空间点索引算法之前，要先谈谈 空间填充曲线(Space-filling curve) 和 分形。解决多维空间点索引需要解决2个问题：第一，如何把多维降为低维或者一维？ 第二，一维的曲线如何分形？

3.1 空间填充曲线

 在数学分析中，有这样一个难题：能否用一条无限长的线，穿过任意维度空间里面的所有点？

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

 一年后，即1891年，希尔伯特就作出了这条曲线，叫 希尔伯特曲线（Hilbert curve）。下图就是1-6阶的希尔伯特曲线。具体构造方式在下一章再说。

 下图是希尔伯特曲线填充满3维空间。

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

 在数学分析中，空间填充曲线是一个参数化的注入函数，它将单位区间映射到单位正方形，立方体，更广义的，n维超立方体等中的连续曲线，随着参数的增加，它可以任意接近单位立方体中的给定点。除了数学重要性之外，空间填充曲线也可用于降维，数学规划，稀疏多维数据库索引，电子学和生物学。空间填充曲线的现在被用在互联网地图中。

3.2 分形

皮亚诺曲线的出现，说明了人们对维数的认识是有缺陷的，有必要重新考察维数的定义。这就是分形几何考虑的问题。在分形几何中，维数可以是分数叫做 分维。

多维空间降维以后，如何 分形，也是一个问题。分形的方式有很多种，这里有一个 列表 ，可以查看如何分形，以及每个分形的分形维数，即 豪斯多夫分形维(Hausdorff fractals dimension)和拓扑维数。这里就不细说分形的问题了，感兴趣的可以仔细阅读链接里面的内容。

3.3 希尔伯特曲线 （Hilbert Curve）

 后面要介绍多维空间点索引算法的理论基础来自希尔伯特曲线，先来仔细说说希尔伯特曲线。

3.3.1 希尔伯特曲线的定义

 希尔伯特曲线一种能填充满一个平面正方形的分形曲线（空间填充曲线），由大卫·希尔伯特在1891年提出。由于它能填满平面，它的豪斯多夫维是2。取它填充的正方形的边长为1，第n步的希尔伯特曲线的长度是 2ⁿ - 2^-n。

3.3.2 希尔伯特曲线的构造方法

 一阶的希尔伯特曲线，生成方法就是把正方形四等分，从其中一个子正方形的中心开始，依次穿线，穿过其余3个正方形的中心。

 二阶的希尔伯特曲线，生成方法就是把之前每个子正方形继续四等分，每4个小的正方形先生成一阶希尔伯特曲线。然后把4个一阶的希尔伯特曲线首尾相连。

 三阶的希尔伯特曲线，生成方法就是与二阶类似，先生成二阶希尔伯特曲线。然后把4个二阶的希尔伯特曲线首尾相连。

 n阶的希尔伯特曲线的生成方法也是递归的，先生成n-1阶的希尔伯特曲线，然后把4个n-1阶的希尔伯特曲线首尾相连。

3.3.3 希尔伯特曲线的特性

 看到这里可能就有读者有疑问了，这么多空间填充曲线，为何要选希尔伯特曲线？因为希尔伯特曲线有非常好的特性。

• 降维。作为空间填充曲线，希尔伯特曲线可以对多维空间有效的降维。

  下图就是希尔伯特曲线在填满一个平面以后，把平面上的点都展开成一维的线了。

  上图里面的希尔伯特曲线只穿了16个点，怎么能代表一个平面呢？当n趋近于无穷大的时候，n阶希尔伯特曲线就可以近似填满整个平面了。

• 稳定。当n阶希尔伯特曲线，n趋于无穷大的时候，曲线上的点的位置基本上趋于稳定。

  举个例子：下图左边是希尔伯特曲线，右边是蛇形的曲线。当n趋于无穷大的时候，两者理论上都可以填满平面。

  在蛇形曲线上给定一个点，当n趋于无穷大的过程中，这个点在蛇形曲线上的位置是时刻变化的。

  但是为何希尔伯特曲线更加优秀呢？再看看希尔伯特曲线，同样是一个点，在n趋于无穷大的情况下：点的位置几乎没有怎么变化。所以希尔伯特曲线更加优秀。

• 连续。希尔伯特曲线是连续的，所以能保证一定可以填满空间。

   连续性是需要数学证明的。具体证明方法这里就不细说了，感兴趣的同学可以看 希尔伯特曲线论文。

4 Google S2 算法

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

4.1 球面坐标转换

 按照之前我们处理多维空间的思路，先考虑如何降维，再考虑如何分形。地球是近似一个球体。球体是一个三维的，如何把三维降成一维呢？

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

						x = r * sin θ * cos φ
						y = r * sin θ * sin φ 
						z = r * cos θ
1
2
3

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

 再进一步，我们可以和球面上的经纬度联系起来。不过这里需要注意的是，纬度的角度 α 和直角坐标系下的球面坐标 θ 加起来等于 90°。所以三角函数要注意转换。于是地球上任意的一个经纬度的点，就可以转换成 f(x,y,z)。在 S2 中，地球半径被当成单位 1 了。所以半径不用考虑。x，y，z 的值域都被限定在了 [-1,1] x [-1,1] x [-1,1] 这个区间之内了。

4.2 球面变平面

 接下来一步 S2 把球面碾成平面。怎么做的呢？

 S2 的投影方案

 首先在地球外面套了一个外切的正方体，如下图：

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

 这里简单的画了一个投影图，上图左边的是投影到正方体一个面的示意图，实际上影响到的球面是右边那张图。

  从侧面看，其中一个球面投影到正方体其中一个面上，边缘与圆心的连线相互之间的夹角为90°，但是和x，y，z轴的角度是45°。我们可以在球的6个方向上，把45°的辅助圆画出来，见下图左边。

 上图左边的图画了6个辅助线，蓝线是前后一对，红线是左右一对，绿线是上下一对。分别都是45°的地方和圆心连线与球面相交的点的轨迹。这样我们就可以把投影到外切正方体6个面上的球面画出来，见上图右边。投影到正方体以后，我们就可以把这个正方体展开了。

 一个正方体的展开方式有很多种。不管怎么展开，最小单元都是一个正方形。

 以上就是 S2 的投影方案。接下来讲讲其他的投影方案。

 首先有下面一种方式，三角形和正方形组合

 这种方式展开图如下图。

 这种方式其实很复杂，构成子图形由两种图形构成。坐标转换稍微复杂一点。

 再还有一种方式是全部用三角形组成

 这种方式三角形个数越多，就能越切近于球体。

 上图最左边的图，由20个三角形构成，可以看的出来，菱角非常多，与球体相差比较大，当三角形个数越来越多，就越来越贴近球体。20个三角形展开以后就可能变成这样。

 最后一种方式可能是目前最好的方式，不过也可能是最复杂的方式。按照六边形来投影。

 六边形的菱角比较少，六个边也能相互衔接其他的六边形。看上图最后边的图可以看出来，六边形足够多以后，非常近似球体。

Uber 在一个公开分享上提到了他们用的是六边形的网格，把城市划分为很多六边形。这块应该是他们自己开发的。
也许滴滴也是划分六边形，也许滴滴有更好的划分方案也说不定。

 在 Google S2 中，它是把地球展开成如下的样子：

 如果上面展开的6个面，假设都用5阶的希尔伯特曲线表示出来的话，6个面会是如下的样子：

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

4.3 球面矩形投影修正

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

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

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

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

 上表中，ToPoint 和 FromPoint 分别是把单位向量转换到 Cell ID 所需要的毫秒数、把 Cell ID 转换回单位向量所需的毫秒数（Cell ID 就是投影到正方体六个面，某个面上矩形的 ID，矩形称为 Cell，它对应的 ID 称为 Cell ID）。ToPointRaw 是某种目的下，把 Cell ID 转换为非单位向量所需的毫秒数。

 在 S2 中默认的转换是二次转换。

#define S2_PROJECTION S2_QUADRATIC_PROJECTION

1
2
3

 详细看看这三种转换到底是怎么转换的。

#if S2_PROJECTION == S2_LINEAR_PROJECTION

inline double S2::STtoUV(double s) {
  return 2 * s - 1;
}

inline double S2::UVtoST(double u) {
  return 0.5 * (u + 1);
}

#elif S2_PROJECTION == S2_TAN_PROJECTION

inline double S2::STtoUV(double s) {
  // Unfortunately, tan(M_PI_4) is slightly less than 1.0.  This isn't due to
  // a flaw in the implementation of tan(), it's because the derivative of
  // tan(x) at x=pi/4 is 2, and it happens that the two adjacent floating
  // point numbers on either side of the infinite-precision value of pi/4 have
  // tangents that are slightly below and slightly above 1.0 when rounded to
  // the nearest double-precision result.

  s = tan(M_PI_2 * s - M_PI_4);
  return s + (1.0 / (GG_LONGLONG(1) << 53)) * s;
}

inline double S2::UVtoST(double u) {
  volatile double a = atan(u);
  return (2 * M_1_PI) * (a + M_PI_4);
}

#elif S2_PROJECTION == S2_QUADRATIC_PROJECTION

inline double S2::STtoUV(double s) {
  if (s >= 0.5) return (1/3.) * (4*s*s - 1);
  else          return (1/3.) * (1 - 4*(1-s)*(1-s));
}

inline double S2::UVtoST(double u) {
  if (u >= 0) return 0.5 * sqrt(1 + 3*u);
  else        return 1 - 0.5 * sqrt(1 - 3*u);
}

#else

#error Unknown value for S2_PROJECTION

#endif

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

 上面有一处对 tan(M_PI_4) 的处理，是因为精度的原因，导致略小于1.0 。所以投影之后的修正函数三种变换应该如下：

// 线性转换
u = 0.5 * ( u + 1)

// tan() 变换
u = 2 / pi * (atan(u) + pi / 4) = 2 * atan(u) / pi + 0.5

// 二次变换
u >= 0，u = 0.5 * sqrt(1 + 3*u)
u < 0,    u = 1 - 0.5 * sqrt(1 - 3*u)

1
2
3
4
5
6
7
8
9
10
11

 注意上面虽然变换公式只写了u，不代表只变换u。实际使用过程中，u，v都分别当做入参，都会进行变换。这块修正函数在 Go 的版本里面就直接只实现了二次变换，其他两种变换方式找遍整个库，根本没有提及。

// stToUV converts an s or t value to the corresponding u or v value.
// This is a non-linear transformation from [-1,1] to [-1,1] that
// attempts to make the cell sizes more uniform.
// This uses what the C++ version calls 'the quadratic transform'.
func stToUV(s float64) float64 {
    if s >= 0.5 {
        return (1 / 3.) * (4*s*s - 1)
    }
    return (1 / 3.) * (1 - 4*(1-s)*(1-s))
}

// uvToST is the inverse of the stToUV transformation. Note that it
// is not always true that uvToST(stToUV(x)) == x due to numerical
// errors.
func uvToST(u float64) float64 {
    if u >= 0 {
        return 0.5 * math.Sqrt(1+3*u)
    }
    return 1 - 0.5*math.Sqrt(1-3*u)
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

 经过修正变换以后，u，v都变换成了s，t。值域也发生了变化。u，v的值域是[-1,1]，变换以后，是s，t的值域是[0,1]。

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

 从所述可以得出，S2 可以优化的点有两处，一是投影的形状能否换成六边形？二是修正的变换函数能否找到一个效果和 tan() 类似的函数，但是计算速度远远高于 tan()，以致于不会影响计算性能？

4.4 点与坐标轴点相互转换

 在 S2 算法中，默认划分 Cell 的等级是30，也就是说把一个正方形划分为 2³⁰ * 2³⁰个小的正方形。那么上一步的s，t映射到这个正方形上面来，对应该如何转换呢？

 原本s，t的值域是[0,1]，现在值域要扩大到[0,2³⁰-1]。

// stToIJ converts value in ST coordinates to a value in IJ coordinates.
func stToIJ(s float64) int {
    return clamp(int(math.Floor(maxSize*s)), 0, maxSize-1)
}
1
2
3
4

 C ++ 的实现版本也一样

inline int S2CellId::STtoIJ(double s) {
  // Converting from floating-point to integers via static_cast is very slow
  // on Intel processors because it requires changing the rounding mode.
  // Rounding to the nearest integer using FastIntRound() is much faster.
  // 这里减去0.5是为了四舍五入
  return max(0, min(kMaxSize - 1, MathUtil::FastIntRound(kMaxSize * s - 0.5)));
}
1
2
3
4
5
6
7

 到这一步，是 h(face,s,t) ==> H(face,i,j) 。

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

 最后一步，如何把 i，j 和希尔伯特曲线上的点关联起来呢？在变换之前，先来解释一下定义的一些变量。

const (
    lookupBits = 4
    swapMask   = 0x01
    invertMask = 0x02
)

var (
    ijToPos = [4][4]int{
        {0, 1, 3, 2}, // canonical order
        {0, 3, 1, 2}, // axes swapped
        {2, 3, 1, 0}, // bits inverted
        {2, 1, 3, 0}, // swapped & inverted
    }
    posToIJ = [4][4]int{
        {0, 1, 3, 2}, // canonical order:    (0,0), (0,1), (1,1), (1,0)
        {0, 2, 3, 1}, // axes swapped:       (0,0), (1,0), (1,1), (0,1)
        {3, 2, 0, 1}, // bits inverted:      (1,1), (1,0), (0,0), (0,1)
        {3, 1, 0, 2}, // swapped & inverted: (1,1), (0,1), (0,0), (1,0)
    }
    posToOrientation = [4]int{swapMask, 0, 0, invertMask | swapMask}
    lookupIJ         [1 << (2*lookupBits + 2)]int
    lookupPos        [1 << (2*lookupBits + 2)]int
)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

 posToIJ 代表的是一个矩阵，里面记录了一些单元希尔伯特曲线的位置信息。把 posToIJ 数组里面的信息用图表示出来，如下图：

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

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

func init() {
    initLookupCell(0, 0, 0, 0, 0, 0)
    initLookupCell(0, 0, 0, swapMask, 0, swapMask)
    initLookupCell(0, 0, 0, invertMask, 0, invertMask)
    initLookupCell(0, 0, 0, swapMask|invertMask, 0, swapMask|invertMask)
}

1
2
3
4
5
6
7
8

 以上是初始化的递归函数，在希尔伯特曲线的标准顺序中可以看到是有4个格子，并且格子都有顺序的，所以初始化要遍历满所有顺序。入参的第4个参数，就是从0 - 3 。

// initLookupCell initializes the lookupIJ table at init time.
func initLookupCell(level, i, j, origOrientation, pos, orientation int) {

    if level == lookupBits {
        ij := (i << lookupBits) + j
        lookupPos[(ij<<2)+origOrientation] = (pos << 2) + orientation
        lookupIJ[(pos<<2)+origOrientation] = (ij << 2) + orientation
    
        return
    }

    level++
    i <<= 1
    j <<= 1
    pos <<= 2
    
    r := posToIJ[orientation]
    
    initLookupCell(level, i+(r[0]>>1), j+(r[0]&1), origOrientation, pos, orientation^posToOrientation[0])
    initLookupCell(level, i+(r[1]>>1), j+(r[1]&1), origOrientation, pos+1, orientation^posToOrientation[1])
    initLookupCell(level, i+(r[2]>>1), j+(r[2]&1), origOrientation, pos+2, orientation^posToOrientation[2])
    initLookupCell(level, i+(r[3]>>1), j+(r[3]&1), origOrientation, pos+3, orientation^posToOrientation[3])
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

 上面这个函数是生成希尔伯特曲线的。我们可以看到有一处对pos << 2的操作，这里是把位置变换到第一个4个小格子中，所以位置乘以了4。由于初始设置的 lookupBits = 4，所以i，j的变化范围是从[0,15]，总共有1616=256个，然后i，j坐标是表示的4个格子，再细分，lookupBits = 4这种情况下能表示的点的个数就是2564=1024个。这也正好是 lookupIJ 和 lookupPos 的总容量。

 画一个局部的图，i，j从0-7变化。

 上图是一个4阶希尔伯特曲线。初始化的实际过程就是初始化4阶希尔伯特上的1024个点的坐标与坐标轴上的x，y轴的对应关系表。举个例子，下表是i，j在递归过程中产生的中间过程。下表是lookupPos 表计算过程。

 取出一行详细分析一下计算过程：

假设当前(i,j)=(0,2)，ij 的计算过程是把 i 左移4位再加上 j，整体结果再左移2位。目的是为了留出2位的方向位置。ij的前4位是i，接着4位是j，最后2位是方向。这样计算出ij的值就是8 。

接着计算lookupPos[i j]的值。从上图中可以看到(0,2)代表的单元格的4个数字是16，17，18，19 。计算到这一步，pos的值为4（pos是专门记录生成格子到第几个了，总共pos的值会循环0-255）。pos代表的是当前是第几个格子(4个小格子组成)，当前是第4个，每个格子里面有4个小格子。所以44就可以偏移到当前格子的第一个数字，也就是16 。posToIJ 数组里面会记录下当前格子的形状。从这里我们从中取出 orientation 。

看上图，16，17，18，19对应的是 posToIJ 数组轴旋转的情况，所以17是位于轴旋转图的数字1代表的格子中。这时 orientation = 1 。

这样 lookupPos[i j] 表示的数字就计算出来了，44+1=17 。这里就完成了i，j与希尔伯特曲线上数字的对应。

 那如何由希尔伯特曲线上的数字对应到实际的坐标呢？

lookupIJ 数组里面记录了反向的信息。lookupIJ 数组 和 lookupPos 数组存储的信息正好是反向的。lookupIJ 数组 下表存的值是 lookupPos 数组 的下表。我们查 lookupIJ 数组 ，lookupIJ[17]的值就是8，对应算出来(i,j)=(0,2)。这个时候的i，j还是大格子。还是需要借助 posToIJ 数组 里面描述的形状信息。当前形状是轴旋转，之前也知道 orientation = 1，由于每个坐标里面有4个小格子，所以一个i，j代表的是2个小格子，所以需要乘以2，再加上形状信息里面的方向，可以计算出实际的坐标 (0 * 2 + 1 , 2 * 2 + 0) = ( 1，4) 。

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

4.6 S2 Cell ID 数据结构

 最后需要来谈谈 S2 Cell ID 数据结构，这个数据结构直接关系到不同 Level 对应精度的问题。

 上图左图中对应的是 Level 30 的情况，右图对应的是 Level 24 的情况。(2的多少次方，角标对应的也就是 Level 的值)

在 S2 中，每个 CellID 是由64位的组成的。可以用一个 uint64 存储。开头的3位表示正方体6个面中的一个，取值范围[0,5]。3位可以表示0-7，但是6，7是无效值。

64位的最后一位是1，这一位是特意留出来的。用来快速查找中间有多少位。从末尾最后一位向前查找，找到第一个不为0的位置，即找到第一个1。这一位的前一位到开头的第4位（因为前3位被占用）都是可用数字。

绿色格子有多少个就能表示划分多少格。上图左图，绿色的有60个格子，于是可以表示[0,230 -1] * [0,230 -1]这么多个格子。上图右图中，绿色格子只有48个，那么就只能表示[0,224 -1]*[0,224 -1]这么多个格子。

 那么不同 level 可以代表的网格的面积究竟是多大呢？

 由上一章我们知道，由于投影的原因，所以导致投影之后的面积依旧有大小差别。这里推算的公式比较复杂，就不证明了，具体的可以看文档。

MinAreaMetric = Metric{2, 8 * math.Sqrt2 / 9} 
AvgAreaMetric = Metric{2, 4 * math.Pi / 6} 
MaxAreaMetric = Metric{2, 2.635799256963161491}

1
2
3
4
5

 下图就是最大面积、最小面积、平均面积的倍数关系。(单位km2是平方公里)

 level 0 就是正方体的六个面之一。地球表面积约等于510,100,000 km²。level 0 的面积就是地球表面积的六分之一。level 30 能表示的最小的面积0.48cm²，最大也就0.93cm² 。

4.7 S2 与 Geohash 对比

 Geohash 有12级，从5000km 到 3.7cm。中间每一级的变化比较大。有时候可能选择上一级会大很多，选择下一级又会小一些。比如选择字符串长度为4，它对应的 cell 宽度是39.1km，需求可能是50km，那么选择字符串长度为5，对应的 cell 宽度就变成了156km，瞬间又大了3倍了。这种情况选择多长的 Geohash 字符串就比较难选。选择不好，每次判断可能就还需要取出周围的8个格子再次进行判断。Geohash 需要 12 bytes 存储。

 S2 有30级，从 0.7cm² 到 85,000,000km² 。中间每一级的变化都比较平缓，接近于4次方的曲线。所以选择精度不会出现 Geohash 选择困难的问题。S2 的存储只需要一个 uint64 即可存下。

 S2 库里面不仅仅有地理编码，还有其他很多几何计算相关的库。地理编码只是其中的一小部分。本文没有介绍到的 S2 的实现还有很多很多，各种向量计算，面积计算，多边形覆盖，距离问题，球面球体上的问题，它都有实现。

 S2 还能解决多边形覆盖的问题。比如给定一个城市，求一个多边形刚刚好覆盖住这个城市。

 如上图，生成的多边形刚刚好覆盖住下面蓝色的区域。这里生成的多边形可以有大有小。不管怎么样，最终的结果也是刚刚覆盖住目标物。

 用相同的 Cell 也可以达到相同的目的，上图就是用相同 Level 的 Cell 覆盖了整个圣保罗城市。
 这些都是 Geohash 做不到的。
 多边形覆盖利用的是近似的算法，虽然不是严格意义上的最优解，但是实践中效果特别好。
 额外值得说明的一点是，Google 文档上强调了，这种多边形覆盖的算法虽然对搜索和预处理操作非常有用，但是“不可依赖”的。理由也是因为是近似算法，并不是唯一最优算法，所以得到的解会依据库的不同版本而产生变化。

4.8 S2 Cell 举例

（1）先来看看经纬度和 CellID 的转换，以及矩形面积的计算

latlng := s2.LatLngFromDegrees(31.232135, 121.41321700000003)
cellID := s2.CellIDFromLatLng(latlng)
cell := s2.CellFromCellID(cellID) //9279882742634381312

// cell.Level()
fmt.Println("latlng = ", latlng)
fmt.Println("cell level = ", cellID.Level())
fmt.Printf("cell = %d\n", cellID)
// 这里 Parent 方法参数可以直接指定返回改点的对应 level 的 CellID。
smallCell := s2.CellFromCellID(cellID.Parent(10))
fmt.Printf("smallCell level = %d\n", smallCell.Level())
fmt.Printf("smallCell id = %b\n", smallCell.ID())
fmt.Printf("smallCell ApproxArea = %v\n", smallCell.ApproxArea())
fmt.Printf("smallCell AverageArea = %v\n", smallCell.AverageArea())
fmt.Printf("smallCell ExactArea = %v\n", smallCell.ExactArea())

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

 上面那些方法打印出来的结果如下：

latlng =  [31.2321350, 121.4132170]
cell level =  30
cell = 3869277663051577529

**** Parent **** 10000000000000000000000000000000000000000
smallCell level = 10
smallCell id = 11010110110010011011110000000000000000000000000000000000000000
smallCell ApproxArea = 1.9611002454714756e-06
smallCell AverageArea = 1.997370817559429e-06
smallCell ExactArea = 1.9611009480261058e-06

1
2
3
4
5
6
7
8
9
10
11
12

（2）再举一个覆盖多边形的例子

 我们先随便创建一个区域。

rect = s2.RectFromLatLng(s2.LatLngFromDegrees(48.99, 1.852))
rect = rect.AddPoint(s2.LatLngFromDegrees(48.68, 2.75))

rc := &s2.RegionCoverer{MaxLevel: 20, MaxCells: 10, MinLevel: 2}
r := s2.Region(rect.CapBound())
covering := rc.Covering(r)

1
2
3
4
5
6
7
8

 覆盖参数设置成 level 是 2 ~20，最多的 Cell 的个数是10个。

 接着我们把 Cell 至多改成20个。

 最后再改成30个。

 可以看到相同的 level 的范围，cell 个数越多越精确目标范围。以上是匹配矩形区域，匹配圆形区域也同理。

 代码就不贴了，与矩形类似。这种功能 Geohash 就做不到，需要自己手动实现了。

（3）最后举一个多边形匹配的例子。

func testLoop() {

    ll1 := s2.LatLngFromDegrees(31.803269, 113.421145)
    ll2 := s2.LatLngFromDegrees(31.461846, 113.695803)
    ll3 := s2.LatLngFromDegrees(31.250756, 113.756228)
    ll4 := s2.LatLngFromDegrees(30.902604, 113.997927)
    ll5 := s2.LatLngFromDegrees(30.817726, 114.464846)
    ll6 := s2.LatLngFromDegrees(30.850743, 114.76697)
    ll7 := s2.LatLngFromDegrees(30.713884, 114.997683)
    ll8 := s2.LatLngFromDegrees(30.430111, 115.42615)
    ll9 := s2.LatLngFromDegrees(30.088491, 115.640384)
    ll10 := s2.LatLngFromDegrees(29.907713, 115.656863)
    ll11 := s2.LatLngFromDegrees(29.783833, 115.135012)
    ll12 := s2.LatLngFromDegrees(29.712295, 114.728518)
    ll13 := s2.LatLngFromDegrees(29.55473, 114.24512)
    ll14 := s2.LatLngFromDegrees(29.530835, 113.717776)
    ll15 := s2.LatLngFromDegrees(29.55473, 113.3772)
    ll16 := s2.LatLngFromDegrees(29.678892, 112.998172)
    ll17 := s2.LatLngFromDegrees(29.941039, 112.349978)
    ll18 := s2.LatLngFromDegrees(30.040949, 112.025882)
    ll19 := s2.LatLngFromDegrees(31.803269, 113.421145)

    point1 := s2.PointFromLatLng(ll1)
    point2 := s2.PointFromLatLng(ll2)
    point3 := s2.PointFromLatLng(ll3)
    point4 := s2.PointFromLatLng(ll4)
    point5 := s2.PointFromLatLng(ll5)
    point6 := s2.PointFromLatLng(ll6)
    point7 := s2.PointFromLatLng(ll7)
    point8 := s2.PointFromLatLng(ll8)
    point9 := s2.PointFromLatLng(ll9)
    point10 := s2.PointFromLatLng(ll10)
    point11 := s2.PointFromLatLng(ll11)
    point12 := s2.PointFromLatLng(ll12)
    point13 := s2.PointFromLatLng(ll13)
    point14 := s2.PointFromLatLng(ll14)
    point15 := s2.PointFromLatLng(ll15)
    point16 := s2.PointFromLatLng(ll16)
    point17 := s2.PointFromLatLng(ll17)
    point18 := s2.PointFromLatLng(ll18)
    point19 := s2.PointFromLatLng(ll19)

    points := []s2.Point{}
    points = append(points, point19)
    points = append(points, point18)
    points = append(points, point17)
    points = append(points, point16)
    points = append(points, point15)
    points = append(points, point14)
    points = append(points, point13)
    points = append(points, point12)
    points = append(points, point11)
    points = append(points, point10)
    points = append(points, point9)
    points = append(points, point8)
    points = append(points, point7)
    points = append(points, point6)
    points = append(points, point5)
    points = append(points, point4)
    points = append(points, point3)
    points = append(points, point2)
    points = append(points, point1)

    loop := s2.LoopFromPoints(points)

    fmt.Println("----  loop search (gets too much) -----")
    // fmt.Printf("Some loop status items: empty:%t   full:%t \n", loop.IsEmpty(), loop.IsFull())

    // ref: https://github.com/golang/geo/issues/14#issuecomment-257064823
    defaultCoverer := &s2.RegionCoverer{MaxLevel: 20, MaxCells: 1000, MinLevel: 1}
    // rg := s2.Region(loop.CapBound())
    // cvr := defaultCoverer.Covering(rg)
    cvr := defaultCoverer.Covering(loop)

    // fmt.Println(poly.CapBound())
    for _, c3 := range cvr {
        fmt.Printf("%d,\n", c3)
    }
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81

 这里用到了 Loop 类，这个类的初始化的最小单元是 Point，Point 是由经纬度产生的。最重要的一点需要注意的是，多边形是按照逆时针方向，左手边区域确定的。

 如果一不小心点是按照顺时针排列的话，那么多边形确定的是外层更大的面，意味着球面除去画的这个多边形以外的都是你想要的多边形。

 举个具体的例子，假如我们想要画的多边形是下图这个样子的：

 如果我们用顺时针的方式依次存储 Point 的话，并用顺时针的这个数组去初始化 Loop，那么就会出现“奇怪”的现象。如下图：

 这张图左上角的顶点和右下角的顶点在地球上是重合的。如果把这个地图重新还原成球面，那么就是整个球面中间挖空了一个多边形。把上图放大，如下图：

 这样就可以很清晰的看到了，中间被挖空了一个多边形。造成这种现象的原因就是按照顺时针的方向存储了每个点，那么初始化一个 Loop 的时候就会选择多边形外圈的更大的多边形。所以，使用 Loop 一定要切记，顺时针表示的是外圈多边形，逆时针表示的是内圈多边形。

 多边形覆盖的问题同之前举的例子一样，相同的 MaxLevel = 20，MinLevel = 1，MaxCells 不同，覆盖的精度就不同，下图是 MaxCells = 100 的情况：

 下图是 MaxCells = 1000 的情况：

 从这个例子也可以看出来 相同的 Level 范围，MaxCells 越精度，覆盖的精度越高。

4.9 S2 的应用

 S2 目前应用比较多，用在和地图相关业务上更多。Google Map 就直接大量使用了 S2 ，速度有多快读者可以自己体验体验。Uber 在搜寻最近的出租车也是用的 S2 算法进行计算的。场景的例子就是本篇文章引语里面提到的场景。滴滴应该也有相关的应用，也许有更加优秀的解法。现在很火的共享单车也会用到这些空间索引算法。

 最后就是外卖行业和地图关联也很密切。美团和饿了么相信也在这方面有很多应用，具体哪里用到了，就请读者自己想象吧。

4.10 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)))
    
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51

参考文章

1.空间索引之 Google S2
2.高效的多维空间点索引算法 — Geohash 和 Google S2
3.geohash算法的实现及可视化（以广州为例）
4.空间索引之GeoHash
5.根据Geohash编码画中心点和边框
6.Python实现geohash编码与解码（TransBigData）