本文是使用 golang 實現 redis 系列的第九篇，主要介紹如何使用 GeoHash 實現搜尋附近的人。
搜尋附近的POI是一個非常常見的功能，它的技術難點在於地理位置是二維的(經緯度)而我們常用的索引(無論是B樹、紅黑樹還是跳錶)都是一維的。GeoHash 演算法的本質就是將二維的經緯度轉換為一維的表示。
本文核心實現程式碼可以在Godis:lib/geohash中找到。也可以下載Godis來親自體驗。
興趣點(Point Of Intererst, POI): 在電子地圖中我們關心的各種地點被稱為POI, 比如餐廳、超市、寫字樓。POI 通常包含名稱、經緯度、描述等資訊。在搜尋附近的人時，你也可以把附近的使用者稱為POI。
GeoHash 原理
我們知道經度的取值範圍是[-180,180], 緯度取值範圍是[-90,90]。我們將經緯度分別二等分，那麼可以將地球表面分為4個部分:
緯度[-90, 0] 用0表示，[0, 90]用1表示；經度[-180, 0] 用0表示，[0, 180]用1表示。經度在前緯度在後，這樣四個部分都有了一個二進位制編碼。
我們對四個小矩形繼續二等分:
兩次等分後的矩形需要用 4bit 來表示。前兩位是第一次等分後所在大矩形的編碼，後兩位表示第二次分割出的小矩形在大矩形中的位置。
對這些矩形進行進一步分割，分割次數越多矩形越小經度越高，最終我們可以得到最夠小的矩形來表示一個點。這個小矩形的編碼可以代替這個點的座標，矩形的邊長就是GeoHash的誤差。
這種分割方式讓我們可以用Z型曲線塗滿整個地圖，這個Z型曲線叫做Peano空間填充曲線。在大多數情況下空間上相鄰的點編碼也非常相似。少數情況下編碼會發生突變，導致編碼相似的點空間距離卻很大(比如上圖中的0111與1000)。
如圖所示除Peano空間填充曲線外還有很多空間填充曲線，其中效果公認較好是Hilbert空間填充曲線。相較於Peano曲線而言，Hilbert曲線沒有較大的突變。但是由於Peano曲線實現更加簡單，所以 Geohash 演算法採用的是Peano空間填充曲線。
實現 GeoHash 編解碼
看過上文的介紹，我們可以很快的寫出GeoHash的編碼過程:
// 返回二進位制編碼和對應矩形的經緯度範圍
func encode0(latitude, longitude float64, bitSize uint) ([]byte, [2][2]float64) {
box := [2][2]float64{ // Geohash的矩形
{-180, 180}, // 經度
{-90, 90}, // 緯度
}
pos := [2]float64{longitude, latitude}
bits := make([]bit, 0)
var precision uint = 0
for precision < bitSize { // 迴圈直到精度足夠
for direction, val := range pos { // 輪流處理經緯度，p.s. 看到這個迴圈了嗎？你可以很方便的把GeoHash推廣到N維空間
mid := (box[direction][0] + box[direction][1]) / 2 // 計算分割點
if val < mid {
// 經(緯)度小於中點，編碼填0，把一下次二分的上界設為當前區間的中點
box[direction][1] = mid
bits = append(bits, 0)
} else {
// 經(緯)度大於中點，編碼填1，把一下次二分的下界設為當前區間的中點
box[direction][0] = mid
bits = append(bits, 1)
}
bit++
precision++
if precision == bitSize {
break
}
}
}
return []byte(bits), box
}
程式碼非常簡單，類似於二分查詢。遺憾的是和大多數語言一樣 golang 操作二進位制資料的最小單位是 byte 而非 bit，所以我們需要額外做一些工作來實現按bit編碼：
// 這才是真正的實現，請關注與上一節程式碼的不同
var bits = []uint8{128, 64, 32, 16, 8, 4, 2, 1}
func encode0(latitude, longitude float64, bitSize uint) ([]byte, [2][2]float64) {
box := [2][2]float64{
{-180, 180}, // lng
{-90, 90}, // lat
}
pos := [2]float64{longitude, latitude}
hash := &bytes.Buffer{}
bit := 0
var precision uint = 0
code := uint8(0)
for precision < bitSize {
for direction, val := range pos {
mid := (box[direction][0] + box[direction][1]) / 2
if val < mid {
box[direction][1] = mid
// 編碼預設為0，不需要操作
} else {
box[direction][0] = mid
code |= bits[bit]
// 通過位或操作寫入1，比如要在位元組的第3位寫入1應該 code |= 32
}
bit++
if bit == 8 { // 計算完一個位元組的編碼，將其寫入buffer
hash.WriteByte(code)
bit = 0
code = 0
}
precision++
if precision == bitSize {
break
}
}
}
// precision 可能無法被 8 整除，此時剩下的二進位制編碼寫到最後
if code > 0 {
hash.WriteByte(code)
}
return hash.Bytes(), box
}
為了方便傳輸GeoHash定義了一種文字格式的編碼，它是將二進位制編碼進行Base32變換後得到的:
// GeoHash的對映表和標準Base32對映表有些不同
var enc = base32.NewEncoding("0123456789bcdefghjkmnpqrstuvwxyz").WithPadding(base32.NoPadding)
func ToString(buf []byte) string {
return enc.EncodeToString(buf)
}
寫完程式碼後可以到www.geohash.cn上測試一下結果是否正確。
跟隨二進位制編碼的指示進行二分既可完成解碼的過程:
func decode0(hash []byte) [][]float64 {
box := [][]float64{
{-180, 180},
{-90, 90},
}
direction := 0
for i := 0; i < len(hash); i++ {
code := hash[i]
for j := 0; j < len(bits); j++ {
mid := (box[direction][0] + box[direction][1]) / 2
mask := bits[j] // 使用掩碼取出指定位
if mask&code > 0 {
// 經(緯)度大於mid
box[direction][0] = mid
} else {
// 經(緯)度小於mid
box[direction][1] = mid
}
direction = (direction + 1) % 2
}
}
return box
}
解碼過程不能得到精確的結果只能得到對應矩形的經緯度範圍，我們通常使用矩形的中心點來作為編碼對應的座標。
搜尋附近
因為位於同一個矩形中的點它們的GeoHash編碼擁有相同的字首，所以我們可以非常容易的找到某個矩形中的所有POI。
理論上我們在使用GeoHash來搜尋附近POI時只需要找到一個合適的矩形然後找出其中所有POI即可，實際上我們面臨兩個問題:
- 即上文提到的GeoHash值突變會導致編碼相近兩個點空間距離很大。
- 若我們的位置在矩形的邊緣, 那麼隔壁矩形裡的POI可能會更近。
若我們位於紅點位置，北面的綠點雖然與我們不在同一個矩形內但顯然更近。
為了解決這些問題，我們除了搜尋定位點所在的矩形外，還搜尋周圍8個區域內的POI。
這樣搜尋附近需要分為兩步來實現:
- 計算所有POI的GeoHash值，並使用跳錶或B+樹等便於進行範圍查詢的資料結構建立索引
- 計算“附近區域”對應的 GeoHash 編碼，找出這些區域內的所有POI
建立空間索引
我們將GeoHash的精度設定為64位，這樣我們就可以將GeoHash編碼轉換成uint64型別存入 SortedSet 資料結構中:
func ToInt(buf []byte) uint64 {
if len(buf) < 8 { // 用0填充不足的位數
buf2 := make([]byte, 8)
copy(buf2, buf)
return binary.BigEndian.Uint64(buf2)
}
return binary.BigEndian.Uint64(buf)
}
因為位於同一矩形的二進位制編碼擁有相同的字首所以我們需要將二進位制編碼的低端作為uint64的高位（即使用大端序），這樣位於同一矩形的uint64編碼都會處於同一個數字區間內。比如0110矩形內所有點的uint64編碼都會處於[0x6000000000000000, 0x7000000000000000)區間內, 結合 SortedSet 的範圍查詢能力我們可以非常迅速地(時間複雜度O(logN))找到所有位於0110區域內的POI。
使用SortedSet 進行索引的程式碼可以在 Godis:db/geo.go 中找到。
找到附近的區域
我們知道地球半徑約為 6371km 那麼第一次分割後我們得到了四個東西寬 6371π km 南北高 3186π km 的矩形，以此遞推在分割10次後我們可以得到寬約40km高約20km的矩形。也就是說 20bit 的 GeoHash 編碼東西方向上的誤差為 40km, 南北方向誤差為 20 km。
我們在wikipedia 上查到 geohash 的誤差範圍:
表格中的geohash length 是 base32 編碼後的字串長度，1個字元可以表示5位(bit)。
我們也可以用程式碼估算指定的距離需要的geohash length:
func estimatePrecisionByRadius(radiusMeters float64, latitude float64) uint {
if radiusMeters == 0 {
return defaultBitSize - 1
}
var precision uint = 1
for radiusMeters < MercatorMax {
radiusMeters *= 2
precision++
}
precision -= 2
if latitude > 66 || latitude < -66 {
precision--
if latitude > 80 || latitude < -80 {
precision--
}
}
if precision < 1 {
precision = 1
}
if precision > 32 {
precision = 32
}
return precision*2 - 1
}
這段程式碼中需要注意兩點:
- 因為地球是球型在繪製地圖時採用墨卡託投影法(Mercator Projection)在靠近南北兩極的地方投影面積比較大（在世界地圖上靠近北極的格陵蘭島看上去非常大），所以在高緯度區域需要減少精度
- 經度的取值範圍是
[-180,180]緯度的取值範圍是[-90,90], 所以在經緯度等分相同次數後我們得到的矩形總是東西長南北短。為了解決這個問題我們返回的精度總是奇數(precision*2 - 1), 這樣經度比緯度多分割一次就可以得到長寬基本相等的矩形了。
接下來我們計算矩形區域對應的uint64編碼的上下界:
func ToRange(scope []byte, precision uint) [2]uint64 {
lower := ToInt(scope)
radius := uint64(1 << (64 - precision))
upper := lower + radius
return [2]uint64{lower, upper}
}
比如位於0110矩形內的點它們的uint64編碼會處於[0110..., 0111...)範圍內，在二進位制編碼中找到合適的位置+1就可以了。
下一步是計算九宮格的GeoHash編碼，我們採用了計算各個矩形中心點經緯度然後重新編碼的實現方式:
func GetNeighbours(latitude, longitude, radiusMeters float64) [][2]uint64 {
precision := estimatePrecisionByRadius(radiusMeters, latitude)
center, box := encode0(latitude, longitude, precision)
height := box[0][1] - box[0][0]
width := box[1][1] - box[1][0]
centerLng := (box[0][1] + box[0][0]) / 2
centerLat := (box[1][1] + box[1][0]) / 2
maxLat := ensureValidLat(centerLat + height)
minLat := ensureValidLat(centerLat - height)
maxLng := ensureValidLng(centerLng + width)
minLng := ensureValidLng(centerLng - width)
var result [10][2]uint64
leftUpper, _ := encode0(maxLat, minLng, precision)
result[1] = ToRange(leftUpper, precision)
upper, _ := encode0(maxLat, centerLng, precision)
result[2] = ToRange(upper, precision)
rightUpper, _ := encode0(maxLat, maxLng, precision)
result[3] = ToRange(rightUpper, precision)
left, _ := encode0(centerLat, minLng, precision)
result[4] = ToRange(left, precision)
result[5] = ToRange(center, precision)
right, _ := encode0(centerLat, maxLng, precision)
result[6] = ToRange(right, precision)
leftDown, _ := encode0(minLat, minLng, precision)
result[7] = ToRange(leftDown, precision)
down, _ := encode0(minLat, centerLng, precision)
result[8] = ToRange(down, precision)
rightDown, _ := encode0(minLat, maxLng, precision)
result[9] = ToRange(rightDown, precision)
return result[1:]
}
也可以通過某個矩形的編碼快速推斷出附近8個矩形的編碼, 這種方式實現難度較高可以參考 Redis 中的實現:
- 入口: geohashNeighbors
- 東西向平移: geohash_move_x
- 南北向平移: geohash_move_y
最後在 SortedSet 中找到 POI:
func geoRadius0(sortedSet *sortedset.SortedSet, lat float64, lng float64, radius float64) redis.Reply {
areas := geohash.GetNeighbours(lat, lng, radius)
members := make([][]byte, 0)
for _, area := range areas {
lower := &sortedset.ScoreBorder{Value: float64(area[0])}
upper := &sortedset.ScoreBorder{Value: float64(area[1])}
elements := sortedSet.RangeByScore(lower, upper, 0, -1, true)
for _, elem := range elements {
members = append(members, []byte(elem.Member))
}
}
return reply.MakeMultiBulkReply(members)
}
OK, 大功告成。