高效的多維空間點索引算法 — Geohash 和 Google S2

時間 2019-12-08

標籤高效多維空間索引算法 geohash google s2 欄目 Google 简体版

原文原文鏈接

引子

天天咱們晚上加班回家，可能都會用到滴滴或者共享單車。打開 app 會看到以下的界面：php

app 界面上會顯示出本身附近一個範圍內可用的出租車或者共享單車。假設地圖上會顯示以本身爲圓心，5千米爲半徑，這個範圍內的車。如何實現呢？最直觀的想法就是去數據庫裏面查表，計算並查詢車距離用戶小於等於5千米的，篩選出來，把數據返回給客戶端。html

這種作法比較笨，通常也不會這麼作。爲何呢？由於這種作法須要對整個表裏面的每一項都計算一次相對距離。太耗時了。既然數據量太大，咱們就須要分而治之。那麼就會想到把地圖分塊。這樣即便每一塊裏面的每條數據都計算一次相對距離，也比以前全表都計算一次要快不少。java

咱們也都知道，如今用的比較多的數據庫 MySQL、PostgreSQL 都原生支持 B+ 樹。這種數據結構能高效的查詢。地圖分塊的過程其實就是一種添加索引的過程，若是能想到一個辦法，把地圖上的點添加一個合適的索引，而且可以排序，那麼就能夠利用相似二分查找的方法進行快速查詢。git

問題就來了，地圖上的點是二維的，有經度和緯度，這如何索引呢？若是隻針對其中的一個維度，經度或者緯度進行搜索，那搜出來一遍之後還要進行二次搜索。那要是更高維度呢？三維。可能有人會說能夠設置維度的優先級，好比拼接一個聯合鍵，那在三維空間中，x，y，z 誰的優先級高呢？設置優先級好像並非很合理。github

本篇文章就來介紹2種比較通用的空間點索引算法。golang

一. GeoHash 算法

1. Genhash 算法簡介

Genhash 是一種地理編碼，由 Gustavo Niemeyer 發明的。它是一種分級的數據結構，把空間劃分爲網格。Genhash 屬於空間填充曲線中的 Z 階曲線（Z-order curve）的實際應用。算法

何爲 Z 階曲線？數據庫

上圖就是 Z 階曲線。這個曲線比較簡單，生成它也比較容易，只須要把每一個 Z 首尾相連便可。數組

Z 階曲線一樣能夠擴展到三維空間。只要 Z 形狀足夠小而且足夠密，也能填滿整個三維空間。bash

說到這裏可能讀者依舊一頭霧水，不知道 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	16	17	18	19	20	21	22	23	24	25	26	27	28	29	30	31
Base 32	0	1	2	3	4	5	6	7	8	9	b	c	d	e	f	g	h	j	k	m	n	p	q	r	s	t	u	v	w	x	y	z

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	19	20	21	22	23	24	25	26	27	28	29	30	31	32	33	34	35
Base 32	2	3	4	5	6	7	8	9	b	B	C	d	D	F	g	G	h	H	j	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
}

複製代碼

4. Geohash 的優缺點

Geohash 的優勢很明顯，它利用 Z 階曲線進行編碼。而 Z 階曲線能夠將二維或者多維空間裏的全部點都轉換成一維曲線。在數學上成爲分形維。而且 Z 階曲線還具備局部保序性。

Z 階曲線經過交織點的座標值的二進制表示來簡單地計算多維度中的點的z值。一旦將數據被加到該排序中，任何一維數據結構，例如二叉搜索樹，B樹，跳躍表或（具備低有效位被截斷）哈希表均可以用來處理數據。經過 Z 階曲線所獲得的順序能夠等同地被描述爲從四叉樹的深度優先遍歷獲得的順序。

這也是 Geohash 的另一個優勢，搜索查找鄰近點比較快。

Geohash 的缺點之一也來自 Z 階曲線。

Z 階曲線有一個比較嚴重的問題，雖然有局部保序性，可是它也有突變性。在每一個 Z 字母的拐角，都有可能出現順序的突變。

看上圖中標註出來的藍色的點點。每兩個點雖然是相鄰的，可是距離相隔很遠。看右下角的圖，兩個數值鄰近紅色的點二者距離幾乎達到了整個正方形的邊長。兩個數值鄰近綠色的點也達到了正方形的一半的長度。

Geohash 的另一個缺點是，若是選擇很差合適的網格大小，判斷鄰近點可能會比較麻煩。

看上圖，若是選擇 Geohash 字符串爲6的話，就是藍色的大格子。紅星是美羅城，紫色的圓點是搜索出來的目標點。若是用 Geohash 算法查詢的話，距離比較近的多是 wtw37p，wtw37r，wtw37w，wtw37m。可是其實距離最近的點就在 wtw37q。若是選擇這麼大的網格，就須要再查找周圍的8個格子。

若是選擇 Geohash 字符串爲7的話，那變成黃色的小格子。這樣距離紅星星最近的點就只有一個了。就是 wtw37qw。

若是網格大小，精度選擇的很差，那麼查詢最近點還須要再次查詢周圍8個點。

二. 空間填充曲線和分形

在介紹第二種多維空間點索引算法以前，要先談談空間填充曲線(Space-filling curve)和分形。

解決多維空間點索引須要解決2個問題，第一，如何把多維降爲低維或者一維？第二，一維的曲線如何分形？

1. 空間填充曲線

在數學分析中，有這樣一個難題：可否用一條無限長的線，穿過任意維度空間裏面的全部點？

在1890年，Giuseppe Peano 發現了一條連續曲線，如今稱爲 Peano 曲線，它能夠穿過單位正方形上的每一個點。他的目的是構建一個能夠從單位區間到單位正方形的連續映射。 Peano 受到 Georg Cantor 早期違反直覺的研究結果的啓發，即單位區間中無限數量的點與任何有限維度歧管（manifold）中無限數量的點，基數相同。 Peano 解決的問題實質就是，是否存在這樣一個連續的映射，一條能填充滿平面的曲線。上圖就是他找到的一條曲線。

通常來講，一維的東西是不可能填滿2維的方格的。可是皮亞諾曲線偏偏給出了反例。皮亞諾曲線是一條連續的但到處不可導的曲線。

皮亞諾曲線的構造方法以下：取一個正方形而且把它分出九個相等的小正方形，而後從左下角的正方形開始至右上角的正方形結束，依次把小正方形的中心用線段鏈接起來；下一步把每一個小正方形分紅九個相等的正方形，而後上述方式把其中中心鏈接起來……將這種操做手續無限進行下去，最終獲得的極限狀況的曲線就被稱做皮亞諾曲線。

皮亞諾對區間[0，1]上的點和正方形上的點的映射做了詳細的數學描述。實際上，正方形的這些點對於

，可找到兩個連續函數 x = f(t) 和 y = g(t)，使得 x 和 y 取屬於單位正方形的每個值。

一年後，即1891年，希爾伯特就做出了這條曲線，叫希爾伯特曲線（Hilbert curve）。

上圖就是1-6階的希爾伯特曲線。具體構造方式在下一章再說。

上圖是希爾伯特曲線填充滿3維空間。

以後還有不少變種的空間填充曲線，龍曲線(Dragon curve)、高斯帕曲線(Gosper curve)、Koch曲線(Koch curve)、摩爾定律曲線(Moore curve)、謝爾賓斯基曲線(Sierpiński curve)、奧斯古德曲線(Osgood curve)。這些曲線和本文無關，就不詳細介紹了。

在數學分析中，空間填充曲線是一個參數化的注入函數，它將單位區間映射到單位正方形，立方體，更廣義的，n維超立方體等中的連續曲線，隨着參數的增長，它能夠任意接近單位立方體中的給定點。除了數學重要性以外，空間填充曲線也可用於降維，數學規劃，稀疏多維數據庫索引，電子學和生物學。空間填充曲線的如今被用在互聯網地圖中。

2. 分形

皮亞諾曲線的出現，說明了人們對維數的認識是有缺陷的，有必要從新考察維數的定義。這就是分形幾何考慮的問題。在分形幾何中，維數能夠是分數叫作分維。

多維空間降維之後，如何分形，也是一個問題。分形的方式有不少種，這裏有一個列表，能夠查看如何分形，以及每一個分形的分形維數，即豪斯多夫分形維(Hausdorff fractals dimension)和拓撲維數。這裏就不細說分形的問題了，感興趣的能夠仔細閱讀連接裏面的內容。

接下來繼續來講多維空間點索引算法，下面一個算法的理論基礎來自希爾伯特曲線，先來仔細說說希爾伯特曲線。

三. Hilbert Curve 希爾伯特曲線

1. 希爾伯特曲線的定義

希爾伯特曲線一種能填充滿一個平面正方形的分形曲線（空間填充曲線），由大衛·希爾伯特在1891年提出。

因爲它能填滿平面，它的豪斯多夫維是2。取它填充的正方形的邊長爲1，第n步的希爾伯特曲線的長度是2^n - 2^(-n)。

2. 希爾伯特曲線的構造方法

一階的希爾伯特曲線，生成方法就是把正方形四等分，從其中一個子正方形的中心開始，依次穿線，穿過其他3個正方形的中心。

二階的希爾伯特曲線，生成方法就是把以前每一個子正方形繼續四等分，每4個小的正方形先生成一階希爾伯特曲線。而後把4個一階的希爾伯特曲線首尾相連。

三階的希爾伯特曲線，生成方法就是與二階相似，先生成二階希爾伯特曲線。而後把4個二階的希爾伯特曲線首尾相連。

n階的希爾伯特曲線的生成方法也是遞歸的，先生成n-1階的希爾伯特曲線，而後把4個n-1階的希爾伯特曲線首尾相連。

3. 爲什麼要選希爾伯特曲線

看到這裏可能就有讀者有疑問了，這麼多空間填充曲線，爲什麼要選希爾伯特曲線？

由於希爾伯特曲線有很是好的特性。

(1) 降維

首先，做爲空間填充曲線，希爾伯特曲線能夠對多維空間有效的降維。

上圖就是希爾伯特曲線在填滿一個平面之後，把平面上的點都展開成一維的線了。

可能有人會有疑問，上圖裏面的希爾伯特曲線只穿了16個點，怎麼能表明一個平面呢？

固然，當n趨近於無窮大的時候，n階希爾伯特曲線就能夠近似填滿整個平面了。

(2) 穩定

當n階希爾伯特曲線，n趨於無窮大的時候，曲線上的點的位置基本上趨於穩定。舉個例子：

上圖左邊是希爾伯特曲線，右邊是蛇形的曲線。當n趨於無窮大的時候，二者理論上均可以填滿平面。可是爲什麼希爾伯特曲線更加優秀呢？

在蛇形曲線上給定一個點，當n趨於無窮大的過程當中，這個點在蛇形曲線上的位置是時刻變化的。

這就形成了點的相對位置始終不定。

再看看希爾伯特曲線，一樣是一個點，在n趨於無窮大的狀況下：

從上圖能夠看到，點的位置幾乎沒有怎麼變化。因此希爾伯特曲線更加優秀。

(3) 連續

希爾伯特曲線是連續的，因此能保證必定能夠填滿空間。連續性是須要數學證實的。具體證實方法這裏就不細說了，感興趣的能夠點文章末尾一篇關於希爾伯特曲線的論文，那裏有連續性的證實。

接下來要介紹的谷歌的 S2 算法就是基於希爾伯特曲線的。如今讀者應該明白選擇希爾伯特曲線的緣由了吧。

四. S² 算法

Google’s S2 library is a real treasure, not only due to its capabilities for spatial indexing but also because it is a library that was released more than 4 years ago and it didn’t get the attention it deserved

上面這段話來自2015年一位谷歌工程師的博文。他由衷的感嘆 S2 算法發佈4年沒有獲得它應有的讚揚。不過如今 S2 已經被各大公司使用了。

在介紹這個重量級算法以前，先解釋一些這個算法的名字由來。S2實際上是來自幾何數學中的一個數學符號 S²，它表示的是單位球。S2 這個庫實際上是被設計用來解決球面上各類幾何問題的。值得提的一點是，除去 golang 官方 repo 裏面的 geo/s2 完成度目前只有40%，其餘語言，Java，C++，Python 的 S2 實現都完成100%了。本篇文章講解以 Go 的這個版本爲主。

接下來就看看怎麼用 S2 來解決多維空間點索引的問題的。

1. 球面座標轉換

按照以前咱們處理多維空間的思路，先考慮如何降維，再考慮如何分形。

衆所周知，地球是近似一個球體。球體是一個三維的，如何把三維降成一維呢？

球面上的一個點，在直角座標系中，能夠這樣表示：

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

複製代碼

一般地球上的點咱們會用經緯度來表示。

再進一步，咱們能夠和球面上的經緯度聯繫起來。不過這裏須要注意的是，緯度的角度 α 和直角座標系下的球面座標 θ 加起來等於 90°。因此三角函數要注意轉換。

因而地球上任意的一個經緯度的點，就能夠轉換成 f(x,y,z)。

在 S2 中，地球半徑被當成單位 1 了。因此半徑不用考慮。x，y，z的值域都被限定在了[-1,1] x [-1,1] x [-1,1]這個區間以內了。

2. 球面變平面

接下來一步 S2 把球面碾成平面。怎麼作的呢？

首先在地球外面套了一個外切的正方體，以下圖。

從球心向外切正方體6個面分別投影。S2 是把球面上全部的點都投影到外切正方體的6個面上。

這裏簡單的畫了一個投影圖，上圖左邊的是投影到正方體一個面的示意圖，實際上影響到的球面是右邊那張圖。

從側面看，其中一個球面投影到正方體其中一個面上，邊緣與圓心的連線相互之間的夾角爲90°，可是和x，y，z軸的角度是45°。咱們能夠在球的6個方向上，把45°的輔助圓畫出來，見下圖左邊。

上圖左邊的圖畫了6個輔助線，藍線是先後一對，紅線是左右一對，綠線是上下一對。分別都是45°的地方和圓心連線與球面相交的點的軌跡。這樣咱們就能夠把投影到外切正方體6個面上的球面畫出來，見上圖右邊。

投影到正方體之後，咱們就能夠把這個正方體展開了。

一個正方體的展開方式有不少種。無論怎麼展開，最小單元都是一個正方形。

以上就是 S2 的投影方案。接下來說講其餘的投影方案。

首先有下面一種方式，三角形和正方形組合。

這種方式展開圖以下圖。

這種方式其實很複雜，構成子圖形由兩種圖形構成。座標轉換稍微複雜一點。

再還有一種方式是所有用三角形組成，這種方式三角形個數越多，就能越切近於球體。

上圖最左邊的圖，由20個三角形構成，能夠看的出來，菱角很是多，與球體相差比較大，當三角形個數愈來愈多，就愈來愈貼近球體。

20個三角形展開之後就可能變成這樣。

最後一種方式多是目前最好的方式，不過也多是最複雜的方式。按照六邊形來投影。

六邊形的菱角比較少，六個邊也能相互銜接其餘的六邊形。看上圖最後邊的圖能夠看出來，六邊形足夠多之後，很是近似球體。

六邊形展開之後就是上面這樣。固然這裏只有12個六邊形。六邊形個數越多越好，粒度越細，就越貼近球體。

Uber 在一個公開分享上提到了他們用的是六邊形的網格，把城市劃分爲不少六邊形。這塊應該是他們本身開發的。也許滴滴也是劃分六邊形，也許滴滴有更好的劃分方案也說不定。

回到 S2 上面來，S2是用的正方形。這樣第一步的球面座標進一步的被轉換成 f(x,y,z) -> g(face,u,v)，face是正方形的六個面，u，v對應的是六個面中的一個面上的x，y座標。

3. 球面矩形投影修正

上一步咱們把球面上的球面矩形投影到正方形的某個面上，造成的形狀相似於矩形，可是因爲球面上角度的不一樣，最終會致使即便是投影到同一個面上，每一個矩形的面積也不大相同。

上圖就表示出了球面上個一個球面矩形投影到正方形一個面上的狀況。

通過實際計算髮現，最大的面積和最小的面積相差5.2倍。見上圖左邊。相同的弧度區間，在不一樣的緯度上投影到正方形上的面積不一樣。

如今就須要修正各個投影出來形狀的面積。如何選取合適的映射修正函數就成了關鍵。目標是能達到上圖右邊的樣子，讓各個矩形的面積儘可能相同。

這塊轉換的代碼在 C++ 的版本里面纔有詳細的解釋，在 Go 的版本里面只一筆帶過了。害筆者懵逼了很久。

	面積比率	邊比率	對角線比率	ToPointRaw	ToPoint	FromPoint
線性變換	5.200	2.117	2.959	0.020	0.087	0.085
tan()變換	1.414	1.414	1.704	0.237	0.299	0.258
二次變換	2.082	1.802	1.932	0.033	0.096	0.108

線性變換是最快的變換，可是變換比最小。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

複製代碼

詳細看看這三種轉換究竟是怎麼轉換的。

#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

複製代碼

上面有一處對 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)

複製代碼

注意上面雖然變換公式只寫了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)
}

複製代碼

通過修正變換之後，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)。目前總共轉換了4步，球面經緯度座標轉換成球面xyz座標，再轉換成外切正方體投影面上的座標，最後變換成修正後的座標。

到目前爲止，S2 能夠優化的點有兩處，一是投影的形狀可否換成六邊形？二是修正的變換函數可否找到一個效果和 tan() 相似的函數，可是計算速度遠遠高於 tan()，以至於不會影響計算性能？

4. 點與座標軸點相互轉換

在 S2 算法中，默認劃分 Cell 的等級是30，也就是說把一個正方形劃分爲 2^30 * 2^30個小的正方形。

那麼上一步的s，t映射到這個正方形上面來，對應該如何轉換呢？

s，t的值域是[0,1]，如今值域要擴大到[0,2^30^-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)
}

複製代碼

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)));
}

複製代碼

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

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
)

複製代碼

在變換以前，先來解釋一下定義的一些變量。

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)
}

複製代碼

這裏是初始化的遞歸函數，在希爾伯特曲線的標準順序中能夠看到是有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])
}

複製代碼

上面這個函數是生成希爾伯特曲線的。咱們能夠看到有一處對pos << 2的操做，這裏是把位置變換到第一個4個小格子中，因此位置乘以了4。

因爲初始設置的lookupBits = 4，因此i，j的變化範圍是從[0,15]，總共有16*16=256個，而後i，j座標是表示的4個格子，再細分，lookupBits = 4這種狀況下能表示的點的個數就是256*4=1024個。這也正好是 lookupIJ 和 lookupPos 的總容量。

畫一個局部的圖，i，j從0-7變化。

上圖是一個4階希爾伯特曲線。初始化的實際過程就是初始化4階希爾伯特上的1024個點的座標與座標軸上的x，y軸的對應關係表。

舉個例子，下表是i，j在遞歸過程當中產生的中間過程。下表是 lookupPos 表計算過程。

(i,j)	ij	ij 計算過程	lookupPos[i j]	lookupPos[i j]計算過程	實際座標
(0,0)	0	0	0	0	(0,0)
(1,0)	64	(116+0)4=64	5	1*4+1=5	(3,0)
(1,1)	68	(116+1)4=68	9	2*4+1=5	(3,2)
(0,1)	4	(016+1)4=4	14	3*4+2=14	(0,2)
(0,2)	8	(016+2)4=8	17	4*4+1=17	(1,4)
(0,3)	12	(016+3)4=12	20	5*4+0=20	(0,6)
(1,3)	76	(116+3)4=76	24	6*4+0=24	(2,6)
(1,2)	72	(116+2)4=72	31	7*4+3=31	(3,4)
(2,2)	136	(216+2)4=136	33	8*4+1=33	(5,4)

取出一行詳細分析一下計算過程。

假設當前(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個小格子。因此4*4就能夠偏移到當前格子的第一個數字，也就是16 。posToIJ 數組裏面會記錄下當前格子的形狀。從這裏咱們從中取出 orientation 。

看上圖，16，17，18，19對應的是 posToIJ 數組軸旋轉的狀況，因此17是位於軸旋轉圖的數字1表明的格子中。這時 orientation = 1 。

這樣 lookupPos[i j] 表示的數字就計算出來了，4*4+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,2^30^-1]區間，最後一步就是把座標系上的點都映射到希爾伯特曲線上。

6. S2 Cell ID 數據結構

最後須要來談談 S2 Cell ID 數據結構，這個數據結構直接關係到不一樣 Level 對應精度的問題。

在 S2 中，每一個 CellID 是由64位的組成的。能夠用一個 uint64 存儲。開頭的3位表示正方體6個面中的一個，取值範圍[0,5]。3位能夠表示0-7，可是6，7是無效值。

64位的最後一位是1，這一位是特地留出來的。用來快速查找中間有多少位。從末尾最後一位向前查找，找到第一個不爲0的位置，即找到第一個1。這一位的前一位到開頭的第4位（由於前3位被佔用）都是可用數字。

綠色格子有多少個就能表示劃分多少格。上圖左圖，綠色的有60個格子，因而能夠表示[0,2^30^ -1] * [0,2^30^ -1]這麼多個格子。上圖右圖中，綠色格子只有36個，那麼就只能表示[0,2^18^ -1]*[0,2^18^ -1]這麼多個格子。

那麼不一樣 level 能夠表明的網格的面積到底是多大呢？

由上一章咱們知道，因爲投影的緣由，因此致使投影以後的面積依舊有大小差異。

這裏推算的公式比較複雜，就不證實了，具體的能夠看文檔。

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

複製代碼

這就是最大最小面積和平均面積的倍數關係。

level 0 就是正方體的六個面之一。地球表面積約等於510,100,000 km^2^。level 0 的面積就是地球表面積的六分之一。level 30 能表示的最小的面積0.48cm^2^，最大也就0.93cm^2^ 。

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 作不到的。

8. S2 Cell 舉例

先來看看經緯度和 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)
	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())


複製代碼

這裏 Parent 方法參數能夠直接指定返回改點的對應 level 的 CellID。

上面那些方法打印出來的結果以下：

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


複製代碼

再舉一個覆蓋多邊形的例子。咱們先隨便建立一個區域。

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)


複製代碼

覆蓋參數設置成 level 2 - 20，最多的 Cell 的個數是10個。

接着咱們把 Cell 至多改爲20個。

最後再改爲30個。

能夠看到相同的 level 的範圍，cell 個數越多越精確目標範圍。

這裏是匹配矩形區域，匹配圓形區域也同理。

代碼就不貼了，與矩形相似。這種功能 Geohash 就作不到，須要本身手動實現了。

9. S2 的應用

S2 目前應用比較多，用在和地圖相關業務上更多。Google Map 就直接大量使用了 S2 ，速度有多快讀者能夠本身體驗體驗。Uber 在搜尋最近的出租車也是用的 S2 算法進行計算的。場景的例子就是本篇文章引語裏面提到的場景。滴滴應該也有相關的應用，也許有更加優秀的解法。如今很火的共享單車也會用到這些空間索引算法。

最後就是外賣行業和地圖關聯也很密切。美團和餓了麼相信也在這方面有不少應用，具體哪裏用到了，就請讀者本身想象吧。