Golang Geohash

Source form

• https://medium.com/@jimcc/%E5%A4%96%E9%80%81%E7%B3%BB-%E6%8A%80%E8%A1%93%E5%81%B4-%E5%9F%BA%E6%96%BCgeohash%E7%AE%97%E6%B3%95%E6%B1%BA%E5%AE%9A%E8%83%BD%E4%B8%8D%E8%83%BD%E4%B8%8B%E5%96%AE-golang%E7%82%BA%E4%BE%8B-16c227032165
• https://mp.weixin.qq.com/s/2B-VJ2xgwxrmsSkE6zuoPA

幾何圖形

Ray-casting

Source: https://rosettacode.org/wiki/Ray-casting_algorithm#Go

Segment solution, task algorithm

The first solution given here follows the model of most other solutions on the page in defining a polygon as a list of segments. Unfortunately this representation does not require that the polygon is closed. Input to the ray-casting algorithm, as noted in the WP article though, is specified to be a closed polygon. The "strange" shape defined here is not a closed polygon and so gives incorrect results against some points. (Graphically it may appear closed but mathematically it needs an additional segment returning to the starting point.)

這個解決方案遵循大部分其他解決方案的模型，將多邊形定義為線段列表。 不幸的是，這種表示方法並不要求多邊形是封閉的。 然而，Ray-casting algorithm(射線投射算法)的輸入被指定為一個封閉的多邊形。 在這個解決方案中，定義的"奇怪形狀"並不是一個封閉的多邊形，因此對某些點給出了錯誤的結果。 （從圖形上看，它可能看起來是封閉的，但從數學上講，它需要一條額外的線段返回到起點。）

Ray-casting 是一種判斷某一個點在不在多邊形內的一個簡單方法 即要判斷的點向任意方向做射線，如果這條射線與多邊形的交點是奇數，則表明這個點位於多邊形內．若是偶數，則表明這個點位於多邊形外. 當射線穿過多邊形的一個頂點時，需要特殊處理，因為此時可能會多次紀錄同一个交點，需要确保只計算一次。

射線投射算法的優點是它是一種簡單且容易實現的算法，可以快速檢測一個點是否在多邊形內部。 此外，射線投射算法適用於各種不規則的多邊形，並且可以處理內部縫隙或洞。 然而，射線投射算法的一個缺點是它可能不適用於某些情況下的多邊形。例如，當多邊形具有複雜的內部結構或自相交時，該算法可能無法正確識別多邊形的邊界。 此外，射線投射算法在處理大型多邊形時可能會變得較慢，因為需要對每一條射線進行計算。

優點: 適用各種凹凸多邊形 缺點: 當商圈很多時, 而每個商圈都是N個點組合成的形狀時, 使用 Ray-casting效率低

package main

import (
    "fmt"
    "math"
)

type xy struct {
    x, y float64
}

type seg struct {
    p1, p2 xy
}

type poly struct {
    name  string
    sides []seg
}

func inside(pt xy, pg poly) (i bool) {
    for _, side := range pg.sides {
        if rayIntersectsSegment(pt, side) {
            i = !i
        }
    }
    return
}

func rayIntersectsSegment(p xy, s seg) bool {
    var a, b xy
    if s.p1.y < s.p2.y {
        a, b = s.p1, s.p2
    } else {
        a, b = s.p2, s.p1
    }
    for p.y == a.y || p.y == b.y {
        p.y = math.Nextafter(p.y, math.Inf(1))
    }
    if p.y < a.y || p.y > b.y {
        return false
    }
    if a.x > b.x {
        if p.x > a.x {
            return false
        }
        if p.x < b.x {
            return true
        }
    } else {
        if p.x > b.x {
            return false
        }
        if p.x < a.x {
            return true
        }
    }
    return (p.y-a.y)/(p.x-a.x) >= (b.y-a.y)/(b.x-a.x)
}

var (
    p1  = xy{0, 0}
    p2  = xy{10, 0}
    p3  = xy{10, 10}
    p4  = xy{0, 10}
    p5  = xy{2.5, 2.5}
    p6  = xy{7.5, 2.5}
    p7  = xy{7.5, 7.5}
    p8  = xy{2.5, 7.5}
    p9  = xy{0, 5}
    p10 = xy{10, 5}
    p11 = xy{3, 0}
    p12 = xy{7, 0}
    p13 = xy{7, 10}
    p14 = xy{3, 10}
)

var tpg = []poly{
    {"square", []seg{{p1, p2}, {p2, p3}, {p3, p4}, {p4, p1}}},
    {"square hole", []seg{{p1, p2}, {p2, p3}, {p3, p4}, {p4, p1},
        {p5, p6}, {p6, p7}, {p7, p8}, {p8, p5}}},
    {"strange", []seg{{p1, p5},
        {p5, p4}, {p4, p8}, {p8, p7}, {p7, p3}, {p3, p2}, {p2, p5}}},
    {"exagon", []seg{{p11, p12}, {p12, p10}, {p10, p13},
        {p13, p14}, {p14, p9}, {p9, p11}}},
}

var tpt = []xy{
    // test points common in other solutions on this page
    {5, 5}, {5, 8}, {-10, 5}, {0, 5}, {10, 5}, {8, 5}, {10, 10},
    // test points that show the problem with "strange"
    {1, 2}, {2, 1},
}

func main() {
    for _, pg := range tpg {
        fmt.Printf("%s:\n", pg.name)
        for _, pt := range tpt {
            fmt.Println(pt, inside(pt, pg))
        }
    }
}

Output:

square:
{5 5} true
{5 8} true
{-10 5} false
{0 5} false
{10 5} true
{8 5} true
{10 10} false
{1 2} true
{2 1} true
square hole:
{5 5} false
{5 8} true
{-10 5} false
{0 5} false
{10 5} true
{8 5} true
{10 10} false
{1 2} true
{2 1} true
strange:
{5 5} true
{5 8} false
{-10 5} false
{0 5} false
{10 5} true
{8 5} true
{10 10} false
{1 2} true
{2 1} false
exagon:
{5 5} true
{5 8} true
{-10 5} false
{0 5} false
{10 5} true
{8 5} true
{10 10} false
{1 2} false
{2 1} false

Closed polygon solution

Here input is given as a list of N vertices defining N segments, where one segment extends from each vertex to the next, and one more extends from the last vertex to the first. In the case of the "strange" shape, this mathematically closes the polygon and allows the program to give correct results.

輸入是一個包含N個頂點的列表，每個頂點定義了一個線段， 其中一個線段從每個頂點延伸到下一個頂點，另一條線段則從最後一個頂點延伸到第一個頂點。 在"奇怪"的形狀中，這個操作在數學上封閉了多邊形，使得程序能夠給出正確的結果。

package main

import (
    "math"
    "fmt"
)

type xy struct {
    x, y float64
}

type closedPoly struct {
    name string
    vert []xy
}

func inside(pt xy, pg closedPoly) bool {
    if len(pg.vert) < 3 {
        return false
    }
    in := rayIntersectsSegment(pt, pg.vert[len(pg.vert)-1], pg.vert[0])
    for i := 1; i < len(pg.vert); i++ {
        if rayIntersectsSegment(pt, pg.vert[i-1], pg.vert[i]) {
            in = !in
        }
    }
    return in
}

func rayIntersectsSegment(p, a, b xy) bool {
    if a.y > b.y {
        a, b = b, a
    }
    for p.y == a.y || p.y == b.y {
        p.y = math.Nextafter(p.y, math.Inf(1))
    }
    if p.y < a.y || p.y > b.y {
        return false
    }
    if a.x > b.x {
        if p.x > a.x {
            return false
        }
        if p.x < b.x {
            return true
        }
    } else {
        if p.x > b.x {
            return false
        }
        if p.x < a.x {
            return true
        }
    }
    return (p.y-a.y)/(p.x-a.x) >= (b.y-a.y)/(b.x-a.x)
}

var tpg = []closedPoly{
    {"square", []xy{{0, 0}, {10, 0}, {10, 10}, {0, 10}}},
    {"square hole", []xy{{0, 0}, {10, 0}, {10, 10}, {0, 10}, {0, 0},
        {2.5, 2.5}, {7.5, 2.5}, {7.5, 7.5}, {2.5, 7.5}, {2.5, 2.5}}},
    {"strange", []xy{{0, 0}, {2.5, 2.5}, {0, 10}, {2.5, 7.5}, {7.5, 7.5},
        {10, 10}, {10, 0}, {2.5, 2.5}}},
    {"exagon", []xy{{3, 0}, {7, 0}, {10, 5}, {7, 10}, {3, 10}, {0, 5}}},
}

var tpt = []xy{{1, 2}, {2, 1}}

func main() {
    for _, pg := range tpg {
        fmt.Printf("%s:\n", pg.name)
        for _, pt := range tpt {
            fmt.Println(pt, inside(pt, pg))
        }
    }
}

Output:

square:
{1 2} true
{2 1} true
square hole:
{1 2} true
{2 1} true
strange:
{1 2} false
{2 1} false
exagon:
{1 2} false
{2 1} false

PNPoly algorithm

This solution replaces the rayIntersectsSegment function above with the expression from the popular PNPoly algorithm described at https://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html. The expression is not only simpler but more accurate.

This solution is preferred over the two above.

這個解決方案使用了流行的PNPoly算法中的表達式來取代上面的rayIntersectsSegment函數。這個表達式不僅更簡單，而且更準確。 這個解決方案比上面的兩個解決方案更為推薦。

package main

import "fmt"

type xy struct {
    x, y float64
}

type closedPoly struct {
    name string
    vert []xy
}

func inside(pt xy, pg closedPoly) bool {
    if len(pg.vert) < 3 {
        return false
    }
    in := rayIntersectsSegment(pt, pg.vert[len(pg.vert)-1], pg.vert[0])
    for i := 1; i < len(pg.vert); i++ {
        if rayIntersectsSegment(pt, pg.vert[i-1], pg.vert[i]) {
            in = !in
        }
    }
    return in
}

func rayIntersectsSegment(p, a, b xy) bool {
    return (a.y > p.y) != (b.y > p.y) &&
        p.x < (b.x-a.x)*(p.y-a.y)/(b.y-a.y)+a.x
}

var tpg = []closedPoly{
    {"square", []xy{{0, 0}, {10, 0}, {10, 10}, {0, 10}}},
    {"square hole", []xy{{0, 0}, {10, 0}, {10, 10}, {0, 10}, {0, 0},
        {2.5, 2.5}, {7.5, 2.5}, {7.5, 7.5}, {2.5, 7.5}, {2.5, 2.5}}},
    {"strange", []xy{{0, 0}, {2.5, 2.5}, {0, 10}, {2.5, 7.5}, {7.5, 7.5},
        {10, 10}, {10, 0}, {2.5, 2.5}}},
    {"exagon", []xy{{3, 0}, {7, 0}, {10, 5}, {7, 10}, {3, 10}, {0, 5}}},
}

var tpt = []xy{{1, 2}, {2, 1}}

func main() {
    for _, pg := range tpg {
        fmt.Printf("%s:\n", pg.name)
        for _, pt := range tpt {
            fmt.Println(pt, inside(pt, pg))
        }
    }
}

Geohash 地理空間數據索引

1. 地理空間數據索引的方法

1. Geohash是一種將地理位置編碼為字符串的方法，它將地球表面切分成經緯度網格，然後根據二進制表示將每個網格編碼為一個字符串。Geohash編碼可以用於快速的地理位置查詢。
2. Google S2是一種基於球面三角形切分的地理空間索引方法，它將地球表面切分成多個球面三角形，然後根據每個球面三角形的唯一標識符進行索引。Google S2可以用於地理位置查詢、地圖匹配和地理數據分析等應用。
3. Uber H3是一種基於六邊形切分的地理空間索引方法，它將地球表面切分成多個六邊形，然後根據每個六邊形的唯一標識符進行索引。Uber H3可以用於地理位置查詢、出行路徑規劃和地理數據分析等應用。
4. R樹（R-tree）：一種用於多維數據索引的數據結構，可以用於地理空間數據索引。通過將數據分割成多個節點，可以實現快速的地理位置查詢。
5. KD樹（KD-tree）：一種二叉樹結構，可以用於多維數據索引。通過將數據點分割成二叉樹，可以實現快速的地理位置查詢。
6. 四叉樹（Quadtree）：一種用於二維數據索引的數據結構，可以用於地理空間數據索引。通過將空間分割成四個象限，可以實現快速的地理位置查詢。
7. 八叉樹（Octree）：一種用於三維數據索引的數據結構，可以用於地理空間數據索引。通過將空間分割成八個八叉子立方體，可以實現快速的地理位置查詢。

2. Geohash 是什麼

Geohash基本原理是將地球理解為一個二維平面，將平面遞歸分解成更小的子塊，每個子塊在一定經緯度范圍內擁有相同的編碼。 [wiki]

GeoHash是一種對地理坐標進行編碼的方法，它將二維坐標映射為一個字符串。 每個字符串代表一個特定的矩形，在該矩形範圍內的所有坐標都共用這個字符串。 字符串越長精度越高，對應的矩形範圍越小。 對一個地理坐標編碼時，按照初始區間範圍緯度[-90,90]和經度[-180,180]，計算目標經度和緯度分別落在左區間還是右區間。落在左區間則取0，右區間則取1。然後， 對上一步得到的區間繼續按照此方法對半查找，得到下一位二進制編碼。 當編碼長度達到業務的進度需求後，根據"偶數位放經度，奇數位放緯度"的規則，將得到的二進制編碼穿插組合，得到一個新的二進制串。 最後，根據base32的對照表，將二進制串翻譯成字符串，即得到地理坐標對應的目標GeoHash字符串。

以坐標"30.280245, 120.027162"為例，計算其GeoHash字符串。 首先對緯度做二進制編碼：

1. 將[-90,90]平分為2部分，"30.280245"落在右區間(0,90]，則第一位取1。
2. 將(0,90]平分為2部分，"30.280245"落在左區間(0,45]，則第二位取0。
3. 不斷重複以上步驟，得到的目標區間會越來越小，區間的兩個端點也越來越逼近"30.280245"。 下圖的流程詳細地描述了前幾次迭代的過程：

來源: 一种基于快速GeoHash实现海量商品与商圈高效匹配的算法

下圖是以大安森林為中心 25.0296, 121.536

1. 將經度和緯度分別轉換為二進制數值。 緯度 25.0296 轉換為 ``

將緯度25.0296按照上述方法進行GeoHash編碼，具體步驟如下： 初始區間範圍為[-90,90]，25.0296落在右區間，因此第一位為1。 | 左端點 | 中間值 | 右端點 | 二進制編碼 | |--------|--------|--------|------------| | -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 | 0 | | 22.5 | 25.312 | 28.125 | 0 | | 22.5 | 23.906 | 25.312 | 1 | | 23.906 | 24.609 | 25.312 | 1 | | 24.609 | 24.961 | 25.312 | 1 | | 24.961 | 25.137 | 25.312 | 0 | | 24.961 | 25.049 | 25.137 | 0 | | 24.961 | 25.005 | 25.049 | 1 | | 25.005 | 25.027 | 25.049 | 1 | | 25.027 | 25.038 | 25.049 | 0 | | 25.027 | 25.033 | 25.038 | 0 |

綜上所述，緯度25.0296的GeoHash二進位編碼為10100 01110 01100。

經度 121.536 轉換為 11010101 將經度121.536按照上述方法進行GeoHash編碼，具體步驟如下：

綜上所述，經度121.536的GeoHash二進位編碼為 11010 11001 10110。

1. 將經度和緯度交替組合起來，形成一個長度為二倍經緯度二進制長度的二進制數，偶數位為經度二進制，奇數位為緯度二進制。 緯度 25.0296 轉換為 10100 01110 01100 經度 121.536 轉換為 11010 11001 10110

111001100010110101101001111000

1. 將上一步得到的二進制數每 5 位分成一組。 11100 11000 10110 10110 10011 11000

根據對照表wiki, 取得geohash是wsqqms 28 24 22 22 19 24

• 工具 : https://www.movable-type.co.uk/scripts/geohash.html

將精準度設為6可得到周圍的geohash | Neighbours | | | |------------|--------|--------| | wsqqmm | wsqqmt | wsqqmv | | wsqqmk | wsqqms | wsqqmu | | wsqqm7 | wsqqme | wsqqmg |

3. Geohash 精準度

如果 length 為 6的 wsqqms 即可知該區域範圍為 1.22km × 0.61km

4. 將商圈轉換成Geohash

4.1 找出商圈的中心點

基於中心點產生geohash，並周圍產生8個neighbors 每一個方塊都是由四個點所組成 因此在判斷任一點(X)是否在區塊內就容易許多, 如果把精準度拉到7, 可發現某幾個點會在範圍外

package main

import (
    "fmt"

    "github.com/mmcloughlin/geohash"
)

type xy struct {
    x, y float64
}

func calCenter(scope []xy) xy {
    sumLng := float64(0)
    sumLat := float64(0)
    for _, v := range scope {
        sumLng += v.x
        sumLat += v.y
    }
    return xy{
        (sumLng) / float64(len(scope)),
        (sumLat) / float64(len(scope)),
    }
}

func checkInBox1(c, p xy, precision int) bool {
    var inBox bool
    centerGeoHash := geohash.Encode(c.x, c.y) //中心點
    // fmt.Println("centerGeoHash = ", centerGeoHash) // wsqqmse661vv

    xGeoHash := geohash.Encode(p.x, p.y)
    xGeoHash = xGeoHash[:precision]

    neighbors := geohash.Neighbors(centerGeoHash[:precision])
    if centerGeoHash[:precision] == xGeoHash {
        return true
    }

    for _, n := range neighbors {
        box := geohash.BoundingBox(n)
        if box.Contains(p.x, p.y) {
            inBox = true
            break
        }
    }

    return inBox
}

func checkInBox2(c, p xy, precision int) bool {

    xGeoHash := geohash.Encode(p.x, p.y)
    xGeoHash = xGeoHash[:precision]

    centerGeoHash := geohash.Encode(c.x, c.y) //中心點
    if centerGeoHash[:precision] == xGeoHash {
        return true
    }

    neighbors := geohash.Neighbors(centerGeoHash[:precision])
    for _, n := range neighbors {
        if n == xGeoHash {
            return true
        }
    }
    return false
}

func main() {
    p1 := xy{25.0346276, 121.5267355} // 東門市場
    p2 := xy{25.0332003, 121.5435299} // 大安捷運站
    p3 := xy{25.0260787, 121.5434607} // 科技大樓捷運站
    p4 := xy{25.0260878, 121.5275484} // 師範大學
    center := calCenter([]xy{p1, p2, p3, p4})
    fmt.Println("center = ", center)                                                         // {25.0299986 121.535318625}
    point := xy{25.029161, 121.538417}                                                       // 台北市立圖書館總館
    fmt.Println("台北市立圖書館總館 checkInBox1 (precision 6)? = ", checkInBox1(center, point, 6))    // true
    fmt.Println("台北市立圖書館總館 checkInBox1 (precision 7)? = ", checkInBox1(center, point, 7))    // false
    fmt.Println("台北市立圖書館總館 checkInBox2 (precision 6)? = ", checkInBox2(center, point, 6))    // true
    fmt.Println("台北市立圖書館總館 checkInBox2 (precision 7)? = ", checkInBox2(center, point, 7))    // true
    point = xy{25.0336876, 121.534786}                                                       // 建國花市
    fmt.Println("建國花市 checkInBox1 (precision 5)? = ", checkInBox1(center, point, 5))         // true
    fmt.Println("建國花市 checkInBox1 (precision 6)? = ", checkInBox1(center, point, 6))         // true
    fmt.Println("建國花市 checkInBox1 (precision 7)? = ", checkInBox1(center, point, 7))         // false
    fmt.Println("建國花市 checkInBox2 (precision 5)? = ", checkInBox2(center, point, 5))         // true
    fmt.Println("建國花市 checkInBox2 (precision 6)? = ", checkInBox2(center, point, 6))         // true
    fmt.Println("建國花市 checkInBox2 (precision 7)? = ", checkInBox2(center, point, 7))         // false
    point = xy{25.0421272, 121.5449393}                                                      // 遠東SOGO 台北忠孝館
    fmt.Println("遠東SOGO 台北忠孝館 checkInBox1 (precision 5)? = ", checkInBox1(center, point, 5)) // true
    fmt.Println("遠東SOGO 台北忠孝館 checkInBox1 (precision 6)? = ", checkInBox1(center, point, 6)) // false
    fmt.Println("遠東SOGO 台北忠孝館 checkInBox1 (precision 7)? = ", checkInBox1(center, point, 7)) // false
    fmt.Println("遠東SOGO 台北忠孝館 checkInBox2 (precision 5)? = ", checkInBox2(center, point, 5)) // true
    fmt.Println("遠東SOGO 台北忠孝館 checkInBox2 (precision 6)? = ", checkInBox2(center, point, 6)) // false
    fmt.Println("遠東SOGO 台北忠孝館 checkInBox2 (precision 7)? = ", checkInBox2(center, point, 7)) // false

優點: 性能高 缺點: 雖然可以覆蓋大範圍區域，但有許多稜角判別會有誤差

結合GeoHash與幾何算法

有沒有什麼方式可以使用空間索引提高性能，又可以降低誤差呢？

標準GeoHash算法只能用來計算二維點坐標對應的GeoHash編碼，我們的場景中還需要計算面數據（即GIS中的POLYGON多邊形對象）對應的GeoHash編碼，需要擴展算法來實現。

算法思路是，先找到目標Polygon的最小外接矩形MBR，計算此MBR西南角坐標對應的GeoHash編碼。 然後用GeoHash編碼的逆算法，反解出此編碼對應的矩形GeoHash塊。 以此GeoHash塊為起點，循環往東、往北找相鄰的同等大小的GeoHash塊，直到找到的GeoHash塊完全超出MBR的範圍才停止。 如此找到的多個GeoHash塊，邊緣上的部分可能與目標Polygon完全不相交，這部分塊需要通過計算剔除掉，如此一來可以減少後續不必要的計算量。

來源: 一种基于快速GeoHash实现海量商品与商圈高效匹配的算法

上面的例子中最終得到的結果高清大圖如下，其中藍色的GeoHash塊是與原始Polygon部分相交的，橘黃色的GeoHash塊是完全被包含在原始Polygon內部的。

流程圖如下:

1. 最小外接矩形(MBR)

最小外接矩形(minimum bounding rectangle)是將任一多邊形，參考多邊形頂點，畫出一個矩形，以商圈來說，如下圖

// 商圈多邊形分別為
// p1 = xy{22.559747, 120.53881}    
// p2 = xy{22.554941, 120.553044}    
// p3 = xy{22.539771, 120.549660}    
// p4 = xy{22.541374, 120.527619}

// 取得最小矩形
func GetMinRectangle(scope []xy) *Rectangle {
    maxLat, maxLng, minLat, minLng  := float64(-90), float64(-180), float64(90), float64(180)
    for _, v := range scope {
        maxLat = math.Max(maxLat, v.x)
        maxLng = math.Max(maxLng, v.y)
        minLat = math.Min(minLat, v.x)
        minLng = math.Min(minLng, v.y)
    }
    r := &Rectangle{
        MaxLat: maxLat,
        MinLat: minLat,
        MaxLng: maxLng,
        MinLng: minLng,
    }
    return r
}
// 得到四個點位置
// (22.559747, 120.527619) 左上
// (22.559747, 120.553044) 右上
// (22.539771, 120.527619) 左下
// (22.539771, 120.553044) 右下

2. 將MBR切成方格(bounding-box)

產生MBR後，即可以根據該矩形從左上角出發(22.559747, 120.527619) 並從左至右(西->東)，從上到下(北->南)依序找出每個格子(bounding-box) 並填滿整個MBR(下圖的精度為6 characters)

從程式的角度要如何獲取MBR的所有格子呢？

從上圖可以知道是由左->右 然後從上->下，因此可以根據每一個點Geohash生成的Bounding-Box循序的找出，因此可以使用兩層遞迴來試試 遞迴的中止條件為: a. 左至右(西->東): 帶入的經度已超出MBR範圍 b. 上至下(北->南): 帶入的緯度已超出MBR範圍

const (
    NORTH = iota
    NORTHEAST
    EAST
    SOUTHEAST
    SOUTH
    SOUTHWEST
    WEST
    NORTHWEST
)

func genBoundingBox(lat, lng float64, direction, precision int64) geohash.Box {
    pGeoHash := geohash.Encode(lat, lng)
    neighbors := geohash.Neighbors(pGeoHash[:precision])
    fmt.Println(neighbors[direction])
    return geohash.BoundingBox(neighbors[direction])
}

// 西邊 -> 東邊
func RecursiveEastToWest(lat, lng float64, maxLng float64, precision int64) {
    if lng > maxLng {
        return
    }
    boundingBox := genBoundingBox(lat, lng, EAST, precision)
    lat, lng = boundingBox.Center()
    RecursiveEastToWest(lat, lng, maxLng, precision)
}

//北->南
func RecursiveNorthToSouth(lat, lng float64, maxLng float64, maxLat float64, precision int64) {
    if lat < maxLat {
        return
    }

    eastBoundingBox := genBoundingBox(lat, lng, EAST, precision)
    elat, elng := eastBoundingBox.Center()
    RecursiveEastToWest(elat, elng, maxLng, precision)

    southBoundingBox := genBoundingBox(lat, lng, SOUTH, precision)
    lat, lng = southBoundingBox.Center()
    RecursiveNorthToSouth(lat, lng, maxLng, maxLat, precision)
}

RecursiveNorthToSouth(22.559747, 120.527619, 120.553044, 22.539771, 6)

3. 格子(Bounding-Box)的種類

下圖是根據MBR所切完的格子，大致可分為三類

• 格子完全覆蓋在商圈內 (如wsj8qb)
• 格子完全在商圈外(如wsj8qd)
• 格子與商圈有相交

3.1 格子完全覆蓋在商圈內

這是最好的狀況了，只要開啟foodpanda時，將經緯度轉換geohash後座落在這個hash值，那100%是在商圈內

3.2 格子完全在商圈外

這一些格子，是不需要的，我們可以認為格子的四個點與多邊形(商圈)都不相交，則該格子必在商圈外，可以根據第一節的 Ray-casting algorithm來處理，然後將這些格子剔除

3.3 格子與商圈有相交

最麻煩的是這種有相交的格子，有些在商圈外，有些在商圈內．

方法1 : 根據格子與多邊形的交點，找出相交的部分，最後形成新的多邊形

方法1實作比較困難， 因此使用方法2 : 若用戶經緯度轉成geohash後屬於格子與商圈有相交的格子，則使用Ray-casting algorithm來處理，只要格子精度夠細，處理這種的case則為少數

// 取得該box矩形四個點
func genBoxRectangle(h geohash.Box) []xy {
    result := []xy{}
    result = append(result, xy{h.MaxLat, h.MinLng})
    result = append(result, xy{h.MaxLat, h.MaxLng})
    result = append(result, xy{h.MinLat, h.MaxLng})
    result = append(result, xy{h.MinLat, h.MinLng})
    return result
}

// 方格與商圈交互關係
// r: 方格是否有在商圈
// cnt: 0: 商圈外, 1-3:與商圈是相交  4: 完全在商圈內
func BoxPolygonIntersectCount(rectangle []xy) (r bool, cnt int) {
    for _, pt := range rectangle {
        if inside(pt, geofencing) {
            r = true
            cnt++
        }
    }
    return r, cnt
}

4. 處理完後的格子(Bounding-Box)

使用精度為7，並剔除在商圈外的格子後，淺黃色格子完全覆蓋在商圈內 淺藍色格子與商圈相交，如下圖

5. 該如何使用?

把所有格子Bounding-Box分類，區分出兩種

格子完全在商圈內 格子與商圈相交

var includeBox []string   //格子完全在商圈內
var intersectBox []string //格子與商圈相交
for _, hash := range allBox {
    box := geohash.BoundingBox(hash)
    rectangle := genBoxRectangle(box)
    if r, c := BoxPolygonIntersectCount(rectangle); r {
        if c == 4 {
            includeBox = append(includeBox, hash)
            continue
        }
        intersectBox = append(intersectBox, hash)
    }
}
fmt.Println(len(includeBox)) 
fmt.Println(len(intersectBox))

完整程式碼

package main

import (
    "fmt"
    "math"

    "github.com/mmcloughlin/geohash"
)

const (
    NORTH = iota
    NORTHEAST
    EAST
    SOUTHEAST
    SOUTH
    SOUTHWEST
    WEST
    NORTHWEST
)

var (
    p1         = xy{25.0346276, 121.5267355} // 東門市場
    p2         = xy{25.0332003, 121.5435299} // 大安捷運站
    p3         = xy{25.0260787, 121.5434607} // 科技大樓捷運站
    p4         = xy{25.0260878, 121.5275484} // 師範大學
    geofencing = poly{
        name:  "商圈",
        sides: []seg{{p1, p2}, {p2, p3}, {p3, p4}, {p4, p1}},
    }
)

type xy struct {
    x, y float64
}
type seg struct {
    p1, p2 xy
}

type poly struct {
    name  string
    sides []seg
}

func inside(pt xy, pg poly) (i bool) {
    for _, side := range pg.sides {
        if rayIntersectsSegment(pt, side) {
            i = !i
        }
    }
    return
}

func rayIntersectsSegment(p xy, s seg) bool {
    var a, b xy
    if s.p1.y < s.p2.y {
        a, b = s.p1, s.p2
    } else {
        a, b = s.p2, s.p1
    }
    for p.y == a.y || p.y == b.y {
        p.y = math.Nextafter(p.y, math.Inf(1))
    }
    if p.y < a.y || p.y > b.y {
        return false
    }
    if a.x > b.x {
        if p.x > a.x {
            return false
        }
        if p.x < b.x {
            return true
        }
    } else {
        if p.x > b.x {
            return false
        }
        if p.x < a.x {
            return true
        }
    }
    return (p.y-a.y)/(p.x-a.x) >= (b.y-a.y)/(b.x-a.x)
}

type Rectangle struct {
    MaxLat float64
    MinLat float64
    MaxLng float64
    MinLng float64
}

func calCenter(scope []xy) xy {
    sumLng := float64(0)
    sumLat := float64(0)
    for _, v := range scope {
        sumLng += v.x
        sumLat += v.y
    }
    return xy{
        (sumLng) / float64(len(scope)),
        (sumLat) / float64(len(scope)),
    }
}

func checkInBox1(c, p xy, precision int) bool {
    var inBox bool
    centerGeoHash := geohash.Encode(c.x, c.y) //中心點
    // fmt.Println("centerGeoHash = ", centerGeoHash) // wsqqmse661vv

    xGeoHash := geohash.Encode(p.x, p.y)
    xGeoHash = xGeoHash[:precision]

    neighbors := geohash.Neighbors(centerGeoHash[:precision])
    if centerGeoHash[:precision] == xGeoHash {
        return true
    }

    for _, n := range neighbors {
        box := geohash.BoundingBox(n)
        if box.Contains(p.x, p.y) {
            inBox = true
            break
        }
    }

    return inBox
}

func checkInBox2(c, p xy, precision int) bool {

    xGeoHash := geohash.Encode(p.x, p.y)
    xGeoHash = xGeoHash[:precision]

    centerGeoHash := geohash.Encode(c.x, c.y) //中心點
    if centerGeoHash[:precision] == xGeoHash {
        return true
    }

    neighbors := geohash.Neighbors(centerGeoHash[:precision])
    for _, n := range neighbors {
        if n == xGeoHash {
            return true
        }
    }
    return false
}

func GetMinRectangle(scope []xy) *Rectangle {
    maxLat, maxLng, minLat, minLng := float64(-90), float64(-180), float64(90), float64(180)
    for _, v := range scope {
        maxLat = math.Max(maxLat, v.x)
        maxLng = math.Max(maxLng, v.y)
        minLat = math.Min(minLat, v.x)
        minLng = math.Min(minLng, v.y)
    }
    r := &Rectangle{
        MaxLat: maxLat,
        MinLat: minLat,
        MaxLng: maxLng,
        MinLng: minLng,
    }
    return r
}

func genBoundingBox(lat, lng float64, direction, precision int64) geohash.Box {
    pGeoHash := geohash.Encode(lat, lng)
    neighbors := geohash.Neighbors(pGeoHash[:precision])
    fmt.Println(neighbors[direction])
    return geohash.BoundingBox(neighbors[direction])
}

// 西邊 -> 東邊
func RecursiveEastToWest(lat, lng float64, maxLng float64, precision int64) {
    if lng > maxLng {
        return
    }
    boundingBox := genBoundingBox(lat, lng, EAST, precision)
    lat, lng = boundingBox.Center()
    RecursiveEastToWest(lat, lng, maxLng, precision)
}

// 北->南
// 左上開始
func RecursiveNorthToSouth(lat, lng float64, maxLng float64, maxLat float64, precision int64) {
    if lat < maxLat {
        return
    }
    eastBoundingBox := genBoundingBox(lat, lng, EAST, precision)
    elat, elng := eastBoundingBox.Center()
    fmt.Println("eastBoundingBox", eastBoundingBox)
    RecursiveEastToWest(elat, elng, maxLng, precision)

    southBoundingBox := genBoundingBox(lat, lng, SOUTH, precision)
    lat, lng = southBoundingBox.Center()
    fmt.Println("southBoundingBox", southBoundingBox)
    RecursiveNorthToSouth(lat, lng, maxLng, maxLat, precision)
}

// 取得該box矩形四個點
func genBoxRectangle(h geohash.Box) []xy {
    result := []xy{}
    result = append(result, xy{h.MaxLat, h.MinLng})
    result = append(result, xy{h.MaxLat, h.MaxLng})
    result = append(result, xy{h.MinLat, h.MaxLng})
    result = append(result, xy{h.MinLat, h.MinLng})
    return result
}

// 方格與商圈交互關係
// r: 方格是否有在商圈
// cnt: 0: 商圈外, 1-3:與商圈是相交  4: 完全在商圈內
func BoxPolygonIntersectCount(rectangle []xy) (r bool, cnt int) {
    for _, pt := range rectangle {
        if inside(pt, geofencing) {
            r = true
            cnt++
        }
    }
    return r, cnt
}

func main() {

    center := calCenter([]xy{p1, p2, p3, p4})
    fmt.Println("center = ", center)                                                         // {25.0299986 121.535318625}
    point := xy{25.029161, 121.538417}                                                       // 台北市立圖書館總館
    fmt.Println("台北市立圖書館總館 checkInBox1 (precision 6)? = ", checkInBox1(center, point, 6))    // true
    fmt.Println("台北市立圖書館總館 checkInBox1 (precision 7)? = ", checkInBox1(center, point, 7))    // false
    fmt.Println("台北市立圖書館總館 checkInBox2 (precision 6)? = ", checkInBox2(center, point, 6))    // true
    fmt.Println("台北市立圖書館總館 checkInBox2 (precision 7)? = ", checkInBox2(center, point, 7))    // true
    point = xy{25.0336876, 121.534786}                                                       // 建國花市
    fmt.Println("建國花市 checkInBox1 (precision 5)? = ", checkInBox1(center, point, 5))         // true
    fmt.Println("建國花市 checkInBox1 (precision 6)? = ", checkInBox1(center, point, 6))         // true
    fmt.Println("建國花市 checkInBox1 (precision 7)? = ", checkInBox1(center, point, 7))         // false
    fmt.Println("建國花市 checkInBox2 (precision 5)? = ", checkInBox2(center, point, 5))         // true
    fmt.Println("建國花市 checkInBox2 (precision 6)? = ", checkInBox2(center, point, 6))         // true
    fmt.Println("建國花市 checkInBox2 (precision 7)? = ", checkInBox2(center, point, 7))         // false
    point = xy{25.0421272, 121.5449393}                                                      // 遠東SOGO 台北忠孝館
    fmt.Println("遠東SOGO 台北忠孝館 checkInBox1 (precision 5)? = ", checkInBox1(center, point, 5)) // true
    fmt.Println("遠東SOGO 台北忠孝館 checkInBox1 (precision 6)? = ", checkInBox1(center, point, 6)) // false
    fmt.Println("遠東SOGO 台北忠孝館 checkInBox1 (precision 7)? = ", checkInBox1(center, point, 7)) // false
    fmt.Println("遠東SOGO 台北忠孝館 checkInBox2 (precision 5)? = ", checkInBox2(center, point, 5)) // true
    fmt.Println("遠東SOGO 台北忠孝館 checkInBox2 (precision 6)? = ", checkInBox2(center, point, 6)) // false
    fmt.Println("遠東SOGO 台北忠孝館 checkInBox2 (precision 7)? = ", checkInBox2(center, point, 7)) // false

    // 取得最小矩形
    scope := []xy{p1, p2, p3, p4}
    r := GetMinRectangle(scope)
    fmt.Println(r) // &{25.0346276 25.0260787 121.5435299 121.5267355}

    // 得到四個點
    // 25.0346276,121.5267355 // 左上
    // 25.0346276,121.5435299 // 右上
    // 25.0260787,121.5435299 // 右下
    // 25.0260787,121.5267355 // 左下

    RecursiveNorthToSouth(25.0346276, 121.5267355, 121.5435299, 25.0346276, 7)
    // RecursiveNorthToSouth(22.559747, 120.527619, 120.553044, 22.539771, 6)

    var includeBox []string   // 格子完全在商圈內
    var intersectBox []string // 格子與商圈相交

    hash1 := geohash.Encode(25.0309655, 121.538374)
    hash1 = hash1[:7]
    hash2 := geohash.Encode(25.0316952, 121.5372275)
    hash2 = hash2[:7]
    hash3 := geohash.Encode(25.0183372, 121.5426943) // 六張犁捷運站
    hash3 = hash3[:7]
    hash4 := geohash.Encode(25.029161, 121.538417)
    hash4 = hash4[:7]
    allBox := []string{hash1, hash2, hash3, hash4}
    for _, hash := range allBox {
        box := geohash.BoundingBox(hash)
        rectangle := genBoxRectangle(box)
        if r, c := BoxPolygonIntersectCount(rectangle); r {
            if c == 4 {
                includeBox = append(includeBox, hash)
                continue
            }
            intersectBox = append(intersectBox, hash)
        }
    }
    fmt.Println("格子完全在商圈內 = ", len(includeBox))  // 2
    fmt.Println("格子與商圈相交 = ", len(intersectBox)) // 0
}

System Design of Doordash: Geo-Hashing and WebSockets for Location Based Services

Function Require :

1. Place an order
2. 10 million users, 500 bytes per user
3. 500000 dashers

1. 在設計 API 時，考慮到餐廳、司機和用戶的位置信息。例如，在下單時，用戶會提供其當前位置的經緯度坐標，以及所需送餐地址的經緯度坐標。在查找附近的司機時，服務器會根據用戶和餐廳的位置信息，以及司機的位置信息進行計算。在將訂單派給司機時，服務器會考慮到餐廳、用戶和司機之間的距離和路線，以選擇最佳的司機。

2. 將餐廳、司機和用戶的位置信息保存在 Redis 緩存中。當用戶下單或取消訂單時，服務器會計算該用戶所在的 GeoHash 值，並在 Redis 緩存中查找附近的司機和餐廳。當有司機接受訂單時，服務器會通過 WebSocket 將司機位置信息傳遞給相應的用戶。在司機接受訂單後，司機位置信息會定期更新(5秒1次)，並在 Redis 緩存中進行更新。

3. 在服務器端定期運行一個後台任務，該任務負責檢查司機和餐廳位置信息是否過期。如果某個司機或餐廳的位置信息超過一定時間沒有更新，該司機或餐廳將被標記為離線狀態。當用戶下單或取消訂單時，服務器將不考慮離線司機或餐廳的位置信息。如果一個司機或餐廳重新上線，它的位置信息將再次保存在 Redis 緩存中。

API 設計

1. 用戶下單 API：POST /orders
2. Request body:
   • user_id: 用戶 ID
   • restaurant_id: 餐廳 ID
   • pickup_address: 取餐地址
   • delivery_address: 送餐地址
   • pickup_latitude: 取餐地址緯度
   • pickup_longitude: 取餐地址經度
   • delivery_latitude: 送餐地址緯度
   • delivery_longitude: 送餐地址經度

3. Response body:

   • order_id: 訂單 ID
   • estimated_delivery_time: 預估送達時間

4. 取消訂單 API：DELETE /orders/{order_id}

5. Request body:
   • user_id: 用戶 ID

6. Response body:

   • message: 成功取消訂單的信息

7. 查詢訂單 API：GET /orders/{order_id}

8. Request body:
   • user_id: 用戶 ID

9. Response body:

   • order_id: 訂單 ID
   • status: 訂單狀態
   • dasher_id: 司機 ID
   • pickup_address: 取餐地址
   • delivery_address: 送餐地址
   • pickup_latitude: 取餐地址緯度
   • pickup_longitude: 取餐地址經度
   • delivery_latitude: 送餐地址緯度
   • delivery_longitude: 送餐地址經度

10. 查找附近司機 API：GET /drivers

11. Request body:
    • user_id: 用戶 ID
    • pickup_latitude: 取餐地址緯度
    • pickup_longitude: 取餐地址經度

12. Response body:

    • driver_id: 司機 ID
    • name: 司機名字
    • phone_number: 司機電話
    • vehicle_type: 車輛類型
    • latitude: 司機位置緯度
    • longitude: 司機位置經度

13. 搜索店家 API：GET /restaurants/search

14. Request parameters:
    • query: 搜索關鍵字
    • latitude: 當前位置緯度
    • longitude: 當前位置經度
    • radius: 搜索半徑（單位：米）
    • limit: 返回結果數量上限

15. Response body:

    • restaurants: 店家列表

16. 查詢店家 API：GET /restaurants/{restaurant_id}

17. Request parameters: None

18. Response body:

    • restaurant_id: 店家 ID
    • name: 店家名稱
    • address: 店家地址
    • phone_number: 店家電話號碼
    • menu_items: 菜單項目列表

19. 新增菜單 API：POST /restaurants/{restaurant_id}/menu

20. Request body:
    • menu_items: 菜單項目

21. Response body:

    • message: 成功新增菜單的信息

22. 更新菜單 API：PUT /restaurants/{restaurant_id}/menu/{menu_id}

23. Request body:
    • menu_items: 菜單項目

24. Response body:

    • message: 成功更新菜單的信息

25. 刪除菜單 API：DELETE /restaurants/{restaurant_id}/menu/{menu_id}

26. Request body:
    • menu_id: 菜單 ID

27. Response body:

    • message: 成功刪除菜單的信息

28. 查詢菜單 API：GET /restaurants/{restaurant_id}/menu

29. Request body: None

30. Response body:

    • menu_items: 菜單項目列表

31. 查詢訂單 API：GET /restaurants/{restaurant_id}/orders/{order_id}

32. Request body: None
33. Response body:
    • order_id: 訂單 ID
    • status: 訂單狀態
    • user_id: 用戶 ID
    • pickup_address: 取餐地址
    • delivery_address: 送餐地址
    • pickup_latitude: 取餐地址緯度
    • pickup_longitude: 取餐地址經度
    • delivery_latitude: 送餐地址緯度
    • delivery_longitude: 送餐地址經度

Redis 設計

• 訂單緩存：當客戶端請求創建訂單時，將訂單信息保存在 Redis 緩存中。當有司機接受訂單時，服務器將更新訂單信息，並從 Redis 中刪除緩存。

  1. 客戶端可能會收到已經被接單的訂單，但是系統中仍然顯示為未接單的狀態，這會對客戶端的體驗產生影響；
  2. 其他司機可能會再次接到已經被接單的訂單，這會導致重複接單的問題；
  3. 如果訂單信息未及時更新，可能會導致系統出現異常或錯誤。

• 位置緩存：緩存用戶和司機的位置信息，以便快速查找附近的司機或用戶。

Database Schema 設計

Users Table

id (int, PK) name (varchar) email (varchar) phone_number (varchar) password (varchar) created_at (datetime) updated_at (datetime)

Restaurants Table

id (int, PK) name (varchar) address (varchar) latitude (float) longitude (float) created_at (datetime) updated_at (datetime)

Menu Items Table

id (int, PK) name (varchar) description (text) price (decimal) created_at (datetime) updated_at (datetime) restaurant_id (int, FK)

Orders Table

id (int, PK) customer_id (int, FK) restaurant_id (int, FK) driver_id (int, FK) status (varchar) total_price (decimal) delivery_address (varchar) delivery_latitude (float) delivery_longitude (float) created_at (datetime) updated_at (datetime)

Order Items Table

id (int, PK) order_id (int, FK) menu_item_id (int, FK) quantity (int) created_at (datetime) updated_at (datetime)

Drivers Table

id (int, PK) name (varchar) email (varchar) phone_number (varchar) password (varchar) status (varchar) latitude (float) longitude (float) created_at (datetime) updated_at (datetime)

Driver Orders Table

id (int, PK) driver_id (int, FK) order_id (int, FK) created_at (datetime) updated_at (datetime)

Reference

• https://rosettacode.org/wiki/Ray-casting_algorithm#Go
• mmcloughlin/geohash
• 外送系-技術01|基於Geohash算法決定能不能下單! (golang為例)
• GeoHash算法学习讲解、解析及原理分析
• Day06 其它資料聚合與geohash
• 判断点是否在多边形内部（Ray-casting Algorithm）
• 一种基于快速GeoHash实现海量商品与商圈高效匹配的算法
• https://www.movable-type.co.uk/scripts/geohash.html