位置情報のハッシュをS2 Geometryライブラリを利用して計算して地球上の多角形を被覆してみた

前回の記事では、S2 Geometryライブラリを使って地球上の円形を覆っていました。

poketo7878-dev.hatenablog.com

その時、「まだ多角形ポリゴンのサポートはされていない」と書いていました。 その後、そのバージョンアップが対応され、内部に穴を持たない多角形(ループと呼ばれていますs2.Loop)のセル計算がサポートされました。

なので、今回はそのあたりについて調査してみましたのでまとめます。

エンコードされたポリラインのデコード

GoogleMapではポリゴンのデータ等をエンコードして文字列に変換することができます。 この方式では直前の点からの差分をつかってエンコードするので、ある程度のサイズのポリゴンデータであっても短くエンコードすることができますので、 データベース等に格納する際はこの形式を検討してみても良いかもしれません。

そんなエンコードされたデータをデコードして、もとの緯度経度のリストに復元してくれるGoライブラリは以下のものがあります:

github.com

ループ

ループはS2ライブラリ上での「最初と最後を結んで閉じられる緯度経度のループ」を意味しています。 緯度経度の順序はループの縁に沿って反時計回り(CCW)に設定されている前提で扱われており、ループの内側がループの覆っている領域ということになります。

このループを作るには、緯度経度をまずはs2.Pointに変換し、そのポイントのリストからLoopを作る必要があります。

coords, _, _ := polyline.DecodeCoords([]byte(encodedLine))

points := make([]s2.Point, 0)
for _, coord := range coords {
    point := s2.PointFromLatLng(s2.LatLngFromDegrees(coord[0], coord[1]))
    points = append(points, point)
}
loop := s2.LoopFromPoints(points)

まずは、先程のライブラリを利用してエンコードされたポリラインをデコードし、それぞれの緯度経度からポイントを作成していますs2.PointFromLatLng そして、そのポイントをpointsにためておいて、s2.LoopFromPointsに渡してやることでs2.Loopが手に入っています。

ループと緯度経度の包含検査

例えば上のような領域をループにしていたとすると、以下のようにして、ループの中に緯度経度が含まれるかどうか判定できます。

shibuya := s2.LatLngFromDegrees(35.658034, 139.701636)
....<先程のコード等>
loop.ContainsPoint(s2.PointFromLatLng(shibuya))

shibuya変数が指しているのは渋谷駅の緯度経度なので、当然結果はtrueとなります。

ループの領域のセルによる被覆

以下のようなコードで、./polyline.txtファイルに書いたエンコードされたポリラインをloopにし、loopの領域を前回と同じようにs2.RegiconCovererを利用して セルのリストに変換します。

package main

import (
    "bufio"
    "fmt"
    "os"

    "github.com/golang/geo/s2"
    polyline "github.com/twpayne/go-polyline"
)

func readPolyline(filePath string) string {
    // ファイルを開く
    f, err := os.Open(filePath)
    if err != nil {
        fmt.Fprintf(os.Stderr, "File %s could not read: %v\n", filePath, err)
        os.Exit(1)
    }

    // 関数return時に閉じる
    defer f.Close()

    scanner := bufio.NewScanner(f)

    var line = ""
    for scanner.Scan() {
        line = scanner.Text()
    }
    if serr := scanner.Err(); serr != nil {
        fmt.Fprintf(os.Stderr, "File %s scan error: %v\n", filePath, err)
    }

    return line
}

func main() {
    encodedLine := readPolyline("./polyline.txt")
    coords, _, _ := polyline.DecodeCoords([]byte(encodedLine))

    rc := &s2.RegionCoverer{MaxLevel: 30, MaxCells: 200}
    points := make([]s2.Point, 0)
    for _, coord := range coords {
        point := s2.PointFromLatLng(s2.LatLngFromDegrees(coord[0], coord[1]))
        points = append(points, point)
    }
    loop := s2.LoopFromPoints(points)
    loopRegion := s2.Region(loop)

    for idx, c := range rc.Covering(loopRegion) {
        fmt.Printf("Cell ID [%d]: %v\n", idx, uint64(c))
    }
}

最大上限50の場合

最大上限100の場合

最大上限200の場合

ほぼ綺麗に覆われているように見えます。

前回から今の間にループのバグが修正されたので、無事多角形領域についてもセルに変換して検索することができるようになりました。

またこれらを使って得られた知見があったら共有いたします。