(cache)S2 Geometryライブラリで遊んでみる

　このエントリはSmartDrive Advent Calendar 2019の一環です。

　Googleが2017年に出したS2というジオメトリ系のライブラリがあり、それのGo実装があるので遊んでみました。

ライブラリの公式サイト

Go言語での実装

　ジオメトリライブラリにはその目的によってアルゴリズムや実装がいくつかあると思いますが、よく使われる用途というと、緯度経度の情報の検索容易性を上げたりデータのサイズダウンを意図したものが多いと思います。この手のライブラリでは、日本だとgeohexをご存知の方は多いのではないでしょうか。geohexのようなGeohash系アルゴリズムの基本的な考え方は、地図上を多角形のタイルで隙間なく敷き詰め、そのタイルの大きさが小さければ小さいほど情報の精度が上がっていく、というものです。そしてそのタイルを使った操作によって、計算を単純化して求めやすくします。

　さてS2についてですが、敷き詰めるものは「Cell」と表現されていて、これは四角形で構成されています。また大きな特徴として、Cellは二次元の地図ではなく真円に対して投影されています。これらを総合して、ある程度計算量を落としてアルゴリズム自体を簡易にしつつ、その上で球全体の歪みを極力少なくするという辺りを狙って設計されているそうです。
　S2についての更に詳しい解説は、Misreading Chat #50「Geodesic Discrete Global Grid Systems」をお聞きいただけますよう、どうぞ宜しくお願いします。

　個人的にですが、S2の魅力はそうしたスケーラブルさだけでなく、APIの豊富さにもあると思っています。そこで、良くありそうな事例を具体的な数字やコードを交えて一つご紹介したいと思います。

　弊社の現在の所在地は日比谷のNTTビルにありますが、ここを中心としていくつかの周辺駅をピックアップしておきます。

　今回は、「日比谷駅」「有楽町駅」「銀座駅」「東銀座駅」「新橋駅」「内幸町駅」「虎ノ門駅」「霞ケ関駅」「桜田門駅」を探索対象とし、ピンを立てた座標を利用します（ピンの位置はGoogleMapでの検索結果を利用したものです。。。）。　ここから周囲500m以内の駅をピックアップしたいなと思ったとき、先ずは RegionCovererというAPIを利用します。　このRegionCovererというのは、ある地図上に指定した図形に対して、その図形を確実に埋めるセルのリストを返すものです。ただこの時、セルの総個数やレベル（セルの大きさ）の範囲は固定しておく必要があります。ここで注意すべきなのは、得られたセルのリストによって構成されるエリアというのは、必ずもともと指定した図形よりオーバーサイズになります。そのため、指定したセルのレベルによっては、本来500m以上離れたポイントなのに検知される、ということがあり得ます。

f:id:sun-basix:20191212123106p:plain — セルのレベルを上げ、総数を増やせば限りなく指定した図形に漸近していくのです。。。

　当然ながら、RegionCovererAPIには制限があり、無限にセル数やレベルを指定することはできません。

f:id:sun-basix:20191212123129p:plain — 500mより大幅に大きいセルのレベルを指定するとこうなる

　またこのように、あまりにセルを大きくしすぎるとミスヒットが増える要因にもなるので、ここはアプリケーションの計算量と相談しながら調整していく必要があります。

f:id:sun-basix:20191212123150p:plain — この辺がほどほどかな

　ここまでのスクリーンショットは、Sidewalk Labsによるデモサイトを利用させていただきました。これ、S2の雰囲気を知るのにとても分かりやすいです。

s2.sidewalklabs.com

　今回はCellレベルをmin, maxともに14に設定して進めてみたいと思います。勿論、maxをもっと上げてより精度の高いエリアを割り出すこともできます。下記のコードの方針は、

まずはRegionCovererによって取得できたCellリストのそれぞれのトークンと、探索対象の地点のCellレベル14のトークンを比較する。マッチするCellがあったら次の段階に進む
中心点（オフィス）と各探索対象地点の距離を求める。もしその距離が500m以内であれば [HIT!!!]をつける

となっています。

https://play.golang.org/p/ZSKFlauQKOl

package main

import (
        "fmt"

        "github.com/golang/geo/s1"
        "github.com/golang/geo/s2"
)

type station struct {
        name   string
        latlng s2.LatLng
}

func toSpecifiedCell(ll s2.LatLng) s2.CellID {
        return s2.CellFromLatLng(ll).ID().Parent(14)
}

// https://www.vcalc.com/equation/?uuid=e6cfcccb-da27-11e2-8e97-bc764e04d25f
const earthRadius float64 = 6378137 // meter

func main() {
        office := s2.LatLngFromDegrees(35.671262, 139.757476)
        stations := []station{
                {"Hibiya", s2.LatLngFromDegrees(35.674867, 139.759596)},
                {"Yurakucho", s2.LatLngFromDegrees(35.674865, 139.762818)},
                {"Ginza", s2.LatLngFromDegrees(35.671706, 139.764472)},
                {"Higashi-ginza", s2.LatLngFromDegrees(35.670010, 139.766721)},
                {"Shimbashi", s2.LatLngFromDegrees(35.666456, 139.758267)},
                {"Uchisaiwaicho", s2.LatLngFromDegrees(35.669388, 139.755331)},
                {"Onarimon", s2.LatLngFromDegrees(35.660473, 139.751265)},
                {"Toranomon", s2.LatLngFromDegrees(35.670113, 139.750109)},
                {"Kasumigaseki", s2.LatLngFromDegrees(35.674789, 139.751890)},
                {"Sakuradamon", s2.LatLngFromDegrees(35.677265, 139.751376)},
        }

        officePoint := s2.PointFromLatLng(office)
        // 地球は（本当は楕円だけど計算が複雑になるので）球形なので、二次元の地図だと円に見えるが、
        // 実態は球体に円を被せるようになっていて、その形が帽子に似てることから「Cap」と呼んでいる。
        // その「帽子の深み」は「求める円の半径/地球の半径」によって算出できる
        officeCap := s2.CapFromCenterAngle(officePoint, s1.Angle(500/earthRadius))
        regionCover := s2.RegionCoverer{
                MinLevel: 14,
                MaxLevel: 14,
                MaxCells: 16,
        }
        coverings := map[s2.CellID]struct{}{}
        for _, cellID := range regionCover.Covering(officeCap) {
                coverings[cellID] = struct{}{}
        }

        for _, st := range stations {
                // coveringsにないということは脚切りすべき対象なので、下記の計算はスキップする
                if _, exists := coverings[toSpecifiedCell(st.latlng)]; !exists {
                        continue
                }
                stPoint := s2.PointFromLatLng(st.latlng)
                // Capの捉え方と同様に、二点間の結び方は直線ではなく、地球が球形であることに影響を受けている。
                // 具体的には、単位球面における2点という扱いとなり、その距離とは弧の長さとなる。
                // その弧の長さに地球の半径をかけることで、実際の近似の距離が算出できる。
                angle := officePoint.Distance(stPoint)
                distance := angle.Radians() * earthRadius
                hit := ""
                if distance <= 500 {
                        hit = " [HIT!!!]"
                }
                fmt.Printf("[%s] distance: %f%s\n", st.name, distance, hit)
        }
}

/*
// 出力結果

[Hibiya] distance: 444.748778 [HIT!!!]
[Yurakucho] distance: 627.884061
[Ginza] distance: 634.597277
[Higashi-ginza] distance: 847.599459
[Shimbashi] distance: 539.762741
[Uchisaiwaicho] distance: 284.865174 [HIT!!!]
[Toranomon] distance: 678.393858
[Kasumigaseki] distance: 639.788421
*/

　上記の算出ロジックにより、オフィスから500m以内にある駅は内幸町駅（約285m）と日比谷駅（約445m）だと分かりました！
　今回はCapを利用した円形の探索を行いましたが、矩形やポリゴンでの探索も可能になっています。なお、上記の出力結果から御成門駅と桜田門駅が含まれていませんが、これはRegionCoverer.Coveringメソッドの結果に含まれるCellIDにこの2つの駅が所属するCellIDが含まれてないので、事前に除外されているという構図です。

　今回のように中心点や探索対象が全てオンメモリにいるのであれば、以下のようなコードで同様の探索も可能です。コード量としてはこちらの方が少なくなります。
（探索対象の駅の各地点がofficeを中心としたCapの中にいるかどうかをContainsPointメソッドによって直接判定しています）
https://play.golang.org/p/GZgk7a2csCj

　スマートドライブではこんな感じで、地図関連のライブラリについても定期的に調査していきたいと思っています。
　弊社もまた全職種募集しておりますので、是非お気軽にお声がけいただく！

smartdrive.co.jp

taiyoh's memorandum

@ttaiyoh が、技術ネタで気づいたことを書き溜めておきます。

S2 Geometryライブラリで遊んでみる