Let's write β

プログラミング中にできたことか、思ったこととか

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

背景

位置情報を扱いながら高速に検索する必要がある場合、DB上で範囲検索が容易にできると便利です。。 高速で位置情報を検索するためには、 GeoHashなどの位置情報のハッシュ計算をしておいて、ハッシュで範囲検索をするとインデックスを活用できて良いです。

今日はその位置情報のハッシュのなかでも高速な処理が求められるサービスでの導入事例の多く、 一方でドキュメントの少ないS2ライブラリについて調査したので、ざっくりまとめました。

S2

S2のライブラリはもともとはGoogleで開発されたライブラリで、 公式の情報としては解説のプレゼンがアップロードされている以外は、古びたソースコードGoogle Codeに上がっているのが分かる程度で、あまり情報はないように思えます。

アルゴリズム解説のドキュメントは、

http://blog.christianperone.com/2015/08/googles-s2-geometry-on-the-sphere-cells-and-hilbert-curve/

こちらの記事がわかりやすく、平面状をヒルベルト曲線で埋めていくことや、エンコードされたCellIDの包含判定が前方一致(厳密には末尾1ビットを削る等の処理は必要だが)でできる特徴などがわかりやすく解説されているので良いです。

S2他言語へのポート

S2のライブラリをGoに移植したものや、Pythonに移植したものが盛んに開発されているので利用すると良いと思います。

今回はGoへの移植版を利用してみます。 開発の現状としては球面幾何の世界を作るための、段階的なステップのR1, R2, R3の部分、S1の一部等は完了しているのですが、S2の部分はまだ本格的に利用するには機能が足りない部分等が多く見られますが、活発に開発はされているようなので、遠くない未来にいい感じになると思います。 PRを送ってみるのもいいかと思います。

円の被覆

S2は球面上での座標を扱ってくれるのですが、 球面上の円というのは、 実際には球の表面に沿って曲がっているのでキャップのような形になっています。 なので、球面上での円に相当する型はs2.Capという名称です(s2の部分はインポートのやり方によって変わりますが)

また、S2は地球に限らず、真球の表面を扱うものなので、球の半径を想定しておらず、基本的には 点と点の角度や、球の大円に対する比率で扱います。

たとえば、s2.CapFromCenterAngleという関数でs2.Capを作ることができるのですが、 引数は、キャップの頂点となる球面上の点と、キャップの縁の大円に対する比率で表現します。

f:id:Pocket7878_dev:20170802005315p:plain

図ですと、赤い部分が球を切ったときのCapです。ちょうど反対側の点まで行くと1で これがちょうどみかんの皮のように球の表面を完全にキャップが覆っている状態です。

なので、たとえば「渋谷駅を中心に半径4kmの円を描きたい」という風に思ったときは、 地球のサイズの球を想定して比率計算をする必要があります。 地球の円周の半分を仮に6371.01kmと置くと、比率は4 / 6371.01になります。

これを考慮してkmでCapを作れるようにするには

func capOnEarth(center s2.Point, radiusKm float64) s2.Cap {
    const earthRadiusKm = 6371.01
    ratio := (radiusKm / earthRadiusKm)
    return s2.CapFromCenterAngle(center, s1.Angle(ratio))
}

こんな風にする必要があります。

さて、次はS2の機能を利用してこの円をCellで埋めてみましょう。

s2.RegionCovererという構造体があり、これは指定された図形をCellで埋める作戦を 計算するための型です。

詳しいセルサイズの説明等はプレゼンを見ていただければ幸いですが、MaxLevel 30というところで、S2が扱える最小サイズのセルまで使ってもよいと指定し、さらにMaxCellsで被覆するにに使うセル数の最大数を指定しています。

shibuya := s2.LatLngFromDegrees(35.658034, 139.701636)
shibuyaCircle := capOnEarth(s2.PointFromLatLng(shibuya), 4)
rc := &s2.RegionCoverer{MaxLevel: 30, MaxCells: 50}
for idx, c := range rc.Covering(shibuyaCircle) {
    cell := s2.CellFromCellID(c)
    fmt.Printf("Cell [%d] vertex:\n", idx)
    for i := 0; i < 4; i++ {
        v := cell.Vertex(i)
        latLng := s2.LatLngFromPoint(v)  
        fmt.Printf("{lat: %v, lng: %v}\n", latLng.Lat, latLng.Lng)
    }
}

でどんなCellで円を埋められるのか計算できます。 この結果をGoogleMapで可視化してみたものが以下になります。

青色の円(重なっているので紫ですが)は GoogleMapの機能で渋谷駅から半径4kmを描いてみたものです。

セル上限50の場合

f:id:Pocket7878_dev:20170802005336p:plain

結構はみ出しているところが目立ちますね

セル上限100の場合

f:id:Pocket7878_dev:20170802005353p:plain

すこしはみ出しが減りました。

セル上限200の場合

f:id:Pocket7878_dev:20170802005412p:plain

かなりはみ出しが改善しました。

ユースケース

たとえば、データベース上に大量の位置情報を格納しておいて、必要に応じて高速に円の中に含まれているデータを取り出したい場合、あらからじめ位置情報をS2でCellIDに変換しておき、今回の円で得られたCellIDと前方一致で調査するという形になります。

今回は円だったのであまり有り難みが伝わりづらいのですが、 これがある地形内(ポリゴン内)に含まれている全データ等を求めたいときにデータベースの範囲インデックスを効かせて計算できるのは非常に高速に求められると思います。

残念ながらポリゴンのサポートはまだGo版のS2ではできていませんが、Python版のライブラリもありますし、MongoDBやUber、Forsquareなどのサービスでも利用されているので、位置情報を扱う何かをするときには検討してみても良いかもしれません。

まえからS2は気になっていて、今日調査してみたのでまとめてみました。さらに何か成果があれば追ってまとめます。

僕が働いているAzit.incでは一緒に働けるエンジニアを募集しています!
採用情報 — 株式会社アジット|Azit Inc.