高效的多维空间点索引算法 — Geohash 和 Google S2

Posted 2020-11-30 daomeidan

作者：一缕殇流化隐半边冰霜
链接：https://www.jianshu.com/p/7332dcb978b2

引子

每天我们晚上加班回家，可能都会用到滴滴或者共享单车。打开 app 会看到如下的界面：

app 界面上会显示出自己附近一个范围内可用的出租车或者共享单车。假设地图上会显示以自己为圆心，5公里为半径，这个范围内的车。如何实现呢？最直观的想法就是去数据库里面查表，计算并查询车距离用户小于等于5公里的，筛选出来，把数据返回给客户端。

这种做法比较笨，一般也不会这么做。为什么呢？因为这种做法需要对整个表里面的每一项都计算一次相对距离。太耗时了。既然数据量太大，我们就需要分而治之。那么就会想到把地图分块。这样即使每一块里面的每条数据都计算一次相对距离，也比之前全表都计算一次要快很多。

我们也都知道，现在用的比较多的数据库 mysql、PostgreSQL 都原生支持 B+ 树。这种数据结构能高效的查询。地图分块的过程其实就是一种添加索引的过程，如果能想到一个办法，把地图上的点添加一个合适的索引，并且能够排序，那么就可以利用类似二分查找的方法进行快速查询。

问题就来了，地图上的点是二维的，有经度和纬度，这如何索引呢？如果只针对其中的一个维度，经度或者纬度进行搜索，那搜出来一遍以后还要进行二次搜索。那要是更高维度呢？三维。可能有人会说可以设置维度的优先级，比如拼接一个联合键，那在三维空间中，x，y，z 谁的优先级高呢？设置优先级好像并不是很合理。

本篇文章就来介绍2种比较通用的空间点索引算法。

---

一. GeoHash 算法

1. Genhash 算法简介

Genhash 是一种地理编码，由 Gustavo Niemeyer 发明的。它是一种分级的数据结构，把空间划分为网格。Genhash 属于空间填充曲线中的 Z 阶曲线（Z-order curve）的实际应用。

何为 Z 阶曲线？

上图就是 Z 阶曲线。这个曲线比较简单，生成它也比较容易，只需要把每个 Z 首尾相连即可。

Z 阶曲线同样可以扩展到三维空间。只要 Z 形状足够小并且足够密，也能填满整个三维空间。

说到这里可能读者依旧一头雾水，不知道 Geohash 和 Z 曲线究竟有啥关系？其实 Geohash算法 的理论基础就是基于 Z 曲线的生成原理。继续说回 Geohash。

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

字符串长度		cell 宽度		cell 高度
1	≤	5,000km	×	5,000km
2	≤	1,250km	×	625km
3	≤	156km	×	156km
4	≤	39.1km	×	19.5km
5	≤	4.89km	×	4.89km
6	≤	1.22km	×	0.61km
7	≤	153m	×	153m
8	≤	38.2m	×	19.1m
9	≤	4.77m	×	4.77m
10	≤	1.19m	×	0.596m
11	≤	149mm	×	149mm
12	≤	37.2mm	×	18.6mm

还记得引语里面提到的问题么？这里我们就可以用 Geohash 来解决这个问题。

我们可以利用 Geohash 的字符串长短来决定要划分区域的大小。这个对应关系可以参考上面表格里面 cell 的宽和高。一旦选定 cell 的宽和高，那么 Geohash 字符串的长度就确定下来了。这样我们就把地图分成了一个个的矩形区域了。

地图上虽然把区域划分好了，但是还有一个问题没有解决，那就是如何快速的查找一个点附近邻近的点和区域呢？

Geohash 有一个和 Z 阶曲线相关的性质，那就是一个点附近的地方(但不绝对) hash 字符串总是有公共前缀，并且公共前缀的长度越长，这两个点距离越近。

由于这个特性，Geohash 就常常被用来作为唯一标识符。用在数据库里面可用 Geohash 来表示一个点。Geohash 这个公共前缀的特性就可以用来快速的进行邻近点的搜索。越接近的点通常和目标点的 Geohash 字符串公共前缀越长（但是这不一定，也有特殊情况，下面举例会说明）

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

Decimal	0	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15
Base 32	0	1	2	3	4	5	6	7	8	9	b	c	d	e	f	g

Decimal	16	17	18	19	20	21	22	23	24	25	26	27	28	29	30	31
Base 32	h	j	k	m	n	p	q	r	s	t	u	v	w	x	y	z

base 36 的版本对大小写敏感，用了36个字符，“23456789bBCdDFgGhHjJKlLMnNPqQrRtTVWX”。

Decimal	0	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18
Base 36	2	3	4	5	6	7	8	9	b	B	C	d	D	F	g	G	h	H	j

Decimal	19	20	21	22	23	24	25	26	27	28	29	30	31	32	33	34	35
Base 36	J	K	I	L	M	n	N	P	q	Q	r	R	t	T	V	W	X

2. Geohash 实际应用举例

接下来的举例以 base-32 为例。举个例子。

上图是一个地图，地图中间有一个美罗城，假设需要查询距离美罗城最近的餐馆，该如何查询？

第一步我们需要把地图网格化，利用 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。一直划分下去。

左区间	中值	右区间	二进制结果
-90	0	90	1
0	45	90	0
0	22.5	45	1
22.5	33.75	45	0
22.5	28.125	33.75	1
28.125	30.9375	33.75	1
30.9375	32.34375	33.75	0
30.9375	31.640625	32.34375	0
30.9375	31.2890625	31.640625	0
30.9375	31.1132812	31.2890625	1
31.1132812	31.2011718	31.2890625	0
31.1132812	31.1572265	31.2011718	1
31.1572265	31.1791992	31.2011718	1
31.1791992	31.1901855	31.2011718	1
31.1901855	31.1956786	31.2011718	0

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

左区间	中值	右区间	二进制结果
-180	0	180	1
0	90	180	1
90	135	180	0
90	112.5	135	1
112.5	123.75	135	0
112.5	118.125	123.75	1
118.125	120.9375	123.75	1
120.9375	122.34375	123.75	0
120.9375	121.640625	122.34375	0
120.9375	121.289062	121.640625	1
121.289062	121.464844	121.640625	0
121.289062	121.376953	121.464844	1
121.376953	121.420898	121.464844	1
121.420898	121.442871	121.464844	0
121.420898	121.431885	121.442871	1

纬度产生的二进制是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。

看到这里，读者应该已经清楚了 Geohash 的算法原理了。咱们把6位和7位都组合到一张图上面来看。

可以看到中间大格子的 Geohash 的值是 wtw37q，那么它里面的所有小格子前缀都是 wtw37q。可以想象，当 Geohash 字符串长度为5的时候，Geohash 肯定就为 wtw37 了。

接下来解释之前说的 Geohash 和 Z 阶曲线的关系。回顾最后一步合并经纬度字符串的规则，“偶数位放经度，奇数位放纬度”。读者一定有点好奇，这个规则哪里来的？凭空瞎想的？其实并不是，这个规则就是 Z 阶曲线。看下图：

x 轴就是纬度，y轴就是经度。经度放偶数位，纬度放奇数位就是这样而来的。

最后有一个精度的问题，下面的表格数据一部分来自 Wikipedia。

Geohash 字符串长度	纬度	经度	纬度误差	经度误差	km误差
1	2	3	±23	±23	±2500
2	5	5	±2.8	±5.6	±630
3	7	8	±0.70	±0.70	±78
4	10	10	±0.087	±0.18	±20
5	12	13	±0.022	±0.022	±2.4
6	15	15	±0.0027	±0.0055	±0.61
7	17	18	±0.00068	±0.00068	±0.076
8	20	20	±0.000085	±0.00017	±0.019
9	22	23
10	25	25
11	27	28
12	30	30

3. Geohash 具体实现

到此，读者应该对 Geohash 的算法都很明了了。接下来用 Go 实现一下 Geohash 算法。

package geohash

import (
    "bytes"
)

const (
    BASE32                = "0123456789bcdefghjkmnpqrstuvwxyz"
    MAX_LATITUDE  float64 = 90
    MIN_LATITUDE  float64 = -90
    MAX_LONGITUDE float64 = 180
    MIN_LONGITUDE float64 = -180
)

var (
    bits   = []int{16, 8, 4, 2, 1}
    base32 = []byte(BASE32)
)

type Box struct {
    MinLat, MaxLat float64 // 纬度
    MinLng, MaxLng float64 // 经度
}

func (this *Box) Width() float64 {
    return this.MaxLng - this.MinLng
}

func (this *Box) Height() float64 {
    return this.MaxLat - this.MinLat
}

// 输入值：纬度，经度，精度(geohash的长度)
// 返回geohash, 以及该点所在的区域
func Encode(latitude, longitude float64, precision int) (string, *Box) {
    var geohash bytes.Buffer
    var minLat, maxLat float64 = MIN_LATITUDE, MAX_LATITUDE
    var minLng, maxLng float64 = MIN_LONGITUDE, MAX_LONGITUDE
    var mid float64 = 0

    bit, ch, length, isEven := 0, 0, 0, true
    for length < precision {
        if isEven {
            if mid = (minLng + maxLng) / 2; mid < longitude {
                ch |= bits[bit]
                minLng = mid
            } else {
                maxLng = mid
            }
        } else {
            if mid = (minLat + maxLat) / 2; mid < latitude {
                ch |= bits[bit]
                minLat = mid
            } else {
                maxLat = mid
            }
        }

        isEven = !isEven
        if bit < 4 {
            bit++
        } else {
            geohash.WriteByte(base32[ch])
            length, bit, ch = length+1, 0, 0
        }
    }

    b := &Box{
        MinLat: minLat,
        MaxLat: maxLat,
        MinLng: minLng,
        MaxLng: maxLng,
    }

    return geohash.String(), b
}