引子
每天我們晚上加班回家，可能都會用到滴滴或者共享單車。打開 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 }