目指せGISマスター #3 メッシュ的なGeocoding手法を調べてみた

久々にGISマスターを目指してみる。
本職じゃないので機会があればねじ込んでいく感じだ。

Geohashが唯一無二の最強だと思っていた

Geohashの存在を知ってからは便利な特徴を活かして使いまくっている。

■便利な理由

• 緯度経度の範囲を1つの文字列で表すことが出来る
• 現在位置の緯度経度をGeohashエンコードすればそのメッシュ内に居るということである
  →DBで範囲検索などしなくても良い
• Geohashの文字数がレベルを表していて、文字数が少ないほど大きいメッシュである
• 子のレベルは親のGeohash文字列を含む
  →つまりGeohashレベル4の文字列で検索する場合は
  検索対象のGeohash文字列の前方4文字が一致すれば良い。
• Python、Rubyなどのライブラリが提供されている
• そのライブラリでは隣接するメッシュが取れる関数がある
• KeyValueStoreなDBなどにGeohash文字列をキーにして値を入れるといい感じ。

日本を覆うGeohashの文字列群をつくっておいて、1メッシュごとにデータを取る、
なんてもことしたことがある。
データ量が多い場合には何らかの範囲でデータを取ってこないといけないので、
その基準にGeohashを使うのが便利すぎた。

また、こういった特徴を活かして匿名化に使うことも提案されている。
k匿名化＋Geohashで、ある条件においてkを満たすレベルのGeohashメッシュを用いるという
ようなもの。
→あえてリンクは貼らないけど。

■もちろん弱点もある

• メッシュの形状が長方形で、レベルによって縦長→横長→縦長となる
  →5bitでBase32エンコーディングしているので
  レベルが上がるごとに(縦2bit,横3bit)→(縦3bit,横2bit)→(縦2bit,横3bit)
  のように分割数が縦横で異なるためである
• 子メッシュが32分割されるのがやや粗い(もちろん用途による)

■実は似たようなものはいろいろある

最近、メッシュ的なものを使うか、という話がでてきて
じゃあGeohashをpushしておくかと思い調べていると
http://qiita.com/kochizufan/items/2fe5f4c9f74636d22ddb
GeohashよりQuadkeyの方が優れている！という主旨の記事を発見した。

そこで本記事ではQuadkey、そして後述のOpen Location Codeについて
実際にライブラリを使ってみながら、
"メッシュ的に使えるかどうか"の観点で記載する。

Quadkey、へぇそういうのもあるのか

井之頭五郎風
BingMapsのタイルindexに用いられているようである。
MS様公式のお出ましだ。
https://msdn.microsoft.com/en-us/library/bb259689.aspx

タイルindexだが、タイルの4隅の位置を取ればメッシュになるわけである。
Geohashと同様に1メッシュに1意のキーを生成でき、
キー1文字が1レベルである。
なるほど。概ね同じ使い方が出来そうである。
上記のサイトにもあるがタイルのindexからキーが生成されている。
つまりは位置情報から一旦タイルindexに変換してキーが生成される。
逆も同様に、キーからタイルindexに変換して範囲が取れる。
そのため「タイル？知らねーよ求めてるのはメッシュだよ」というユーザにも一旦タイルindexを意識させてくる。ちょっと煩わしい。

でも、Quadkeyは正方形で、4分割なのである

これは素晴らしい。
4分割なので4進数で表現される。キーも0～3でできている。
Maxレベルは20。悪くない粒度だと思う。
Geohash信者だった私も思わず「あ、負けたな」と思ってしまった。
汎用性では完全にQuadkeyという印象を受けた。

■BingMapsで実際に使っている

この画像のURLが

@2017 Microsoft
下である。
https://t0.ssl.ak.dynamic.tiles.virtualearth.net/comp/ch/13300211231022?mkt=ja-JP&it=G,BX,L,LA&shading=hill&og=94&n=z&c4w=1

13300211231022 これがこのタイルのQuadkeyだ。
URLにもしっかりと組み込まれている。

皇居なのは深い意味は無い。
あるとすれば、眞子様ご婚約おめでとうございます。ってこと。

■Quadkeyから緯度経度の範囲を取得してみる

Pythonにはquadkeyというパッケージがあるが、Python3のためか私の環境では動かず。
mercantileパッケージならば動くためこちらを使用した。

ざっと書いてみる。

#mercantileをインポートする
import mercantile

#quadkeyから4方向の位置を取得する
def make_bounds_by_quadkey(t_quadkey):
    #quadkeyからタイルindexを取得
    t_tile = mercantile.quadkey_to_tile(t_quadkey)
    #タイルindexから4方向の端点を取得
    t_box = mercantile.bounds(t_tile.x, t_tile.y,t_tile.z)

    return t_box

if __name__ == "__main__":
    out_str=[]
    #ヘッダ行
    out_str.append("south,east,north,west,Quadkey,zoom")

    #範囲を取り出すQuadkeyの親配列
    #子はzoom_append_rangeの長さで自動生成して追加する
    t_quadkeys=[
        "13300211231022"
    ]

    #Quadkeyから緯度経度の範囲を算出
    for t_quadkey in t_quadkeys:
        t_box = make_bounds_by_quadkey(t_quadkey)
        print([t_quadkey,t_box])
        #こんな感じで取れる
        #['13300211231022',
        # LngLatBbox(west=139.74609375, south=35.67514743608467,
        #            east=139.76806640625, north=35.6929946320988)]

        #わかりやすく書く主義だから
        bounds_str = str(t_box.south) + "," + str(t_box.east) + "," + str(t_box.north) + ","+ str(t_box.west)

        #ズームレベル=Quadkeyの文字数
        t_zoom=len(t_quadkey)

        out_str.append(bounds_str+","+t_quadkey+","+str(t_zoom))

    #出力
    f = open("quadkey_bounds.csv", 'w')
    f.write("\n".join(out_str) + "\n")
    f.close()

出てきたCSVを可視化してみる。
south,east,north,westから4点作って可視化すればよい。
→可視化方法は割愛。OSM+Openlayersのツールを自前で作った。

対象は前述のコード13300211231022 。

BingMap	Quadkeyから変換して可視化
@2017 Microsoft	@OpenStreetMap

左はタイル1枚分、右は緑枠がQuadkeyから変換した範囲。
比べてみると範囲は一致しているので間違いない。
BingMapsさん、Quadkeyレベル14。使ってるで～～。

■子のレベルも見てみる

レベルを15にするとこうなる。
上のソースのt_quadkeys=～の部分を書き換える。

    t_quadkeys=[
        "133002112310220",
        "133002112310221",
        "133002112310222",
        "133002112310223"
    ]

末尾に0～3を追加するだけ。
可視化すると下図だ。

こうして4分割されていく。
Geohashだとレベルを上げる時、手動で上げるのが難しいが、
Quadkeyの4分割なら手動でも余裕。
「もう1レベル上げて、この部分のメッシュが欲しい！」みたいなケースが
データ解析にはあり得る要望なので
直感的にできるのは小さくないメリットだと思う。

Quadkey、これはいいぞ。
なんて満足していたけど、GoogleもGeoCoding手法を開発していて、
無視できないということがわかってしまった。

GoogleのOpen Location Code、ほほー、そうくる

井之頭五郎風

地球上のあらゆる地点をピンポイントで特定するコードOpen Location Code(OLC)
plus+codes
https://plus.codes/
→メニューからグリッドをonにしましょう
やGoogleMapで使えるとのこと。
https://japan.cnet.com/article/35068908/

■どんなもの？

ドア位置レベルで指定して位置をコード化できるとのこと。
コードの階層は以下で決まっている。

階層	範囲
region code	約100km四方の範囲を特定する
city code	5km四方の範囲を特定する
neighbourhood code	250m四方の範囲を特定する
building cod	14m四方の範囲を特定する
door	もっと小さい範囲

5段階。Quadkeyの20段階に比べると粗く、
1階層深くすると400分割とかになったりするので汎用性にはかける。
だが、

変にどのレベルで切るか悩むくらいならこれくらいスパッと決まってた方が良いのかもしれない。

これをどうやって利用するのか、サービスの性質によりますが。
Googleが"この範囲でcity codeだ"といえば"じゃあいいか"、となるかもしれない。

■ともかく使ってみる。

OLCはGithubでライブラリが公開されている。
https://github.com/google/open-location-code

pythonのライブラリを使って簡単な動作確認をしてみる。
上記のgithubからopenlocationcode.pyを取得し、
下記のスクリプトを作成して実行する。
→皇居付近の1つのcityとその子となるneighbourhoodのメッシュを
　取得する処理である。

※2017/08/24 切り分け間違えてました。cityで大枠切ってneighbourhoodが小さい枠でした。
さらにその下にbuildingとdoorを作れるということになります。

import openlocationcode
if __name__ == "__main__":
    #皇居付近の位置
    latitude=35.68407104
    longitude=139.7570801

    PAIR_CODE_LENGTH_=10
    #子のコードを生成する際に使用する
    CODE_ALPHABET_ = '23456789CFGHJMPQRVWX'

    out_str=[]
    out_str.append("south,east,north,west,Quadkey,zoom")

    #10未満はエラーになる
    length_list = [10]
    code_list = []
    for ll in length_list:
        tcode = openlocationcode.encode(latitude, longitude, ll)
        print(tcode)
        code_list.append(tcode)

    #decode対象コードリスト
    neighbourhood_list= []
    use_child = True
    for cl in code_list:
        #buildingまでで切る
        split_tcode = cl.split("+")

        #cityで切るなら末尾の2文字を00にする
        city_code = split_tcode[0][0:6]
        neighbourhood_list.append(city_code +"00+")

        #neighbourhoodを全取得する場合はCODE_ALPHABET_ * CODE_ALPHABET_を生成して追加する
        if use_child:
            for ca1 in range(len(CODE_ALPHABET_)):
                for ca2 in range(len(CODE_ALPHABET_)):
                    neighbourhoodcode=city_code +CODE_ALPHABET_[ca1]+CODE_ALPHABET_[ca2]
                    neighbourhood_list.append(neighbourhoodcode+"+")

    #decodeする
    for nc in neighbourhood_list:
        decode = openlocationcode.decode(nc)
        print(decode)

        #CodeAreaクラスから必要なデータをとって出力する
        bounds_str = str(decode.latitudeLo) + "," + str(decode.longitudeLo) + "," + str(decode.latitudeHi) + ","+ str(decode.longitudeHi)
        t_zoom=decode.codeLength

        out_str.append(bounds_str+","+nc+","+str(t_zoom))

    ##出力
    f = open("quadkey_bounds_olc.csv", 'w')
    f.write("\n".join(out_str) + "\n")
    f.close()

city	neighbourhood

やはりこのサイズ感だと使いにくいか。
例えばcityの範囲でデータを検索するとなると相当数ヒットしそうだが、
neighbourhoodレベルだと小さすぎるので複数の範囲で検索する必要が出てくる。

まとめ(超主観なので注意されたし)

みんな違ってみんな良い。
と言いたいところだが、性質上Quadkeyの汎用性の高さが良さそう。
プライバシー対策でk匿名化に応用するにしてもGeohashだと32倍ずつのメッシュで
評価するのに対し、4倍ずつのメッシュで評価する方が小さいメッシュで表現できる可能性があり、
データの利用者もニッコリなのではないかと考える。
検索時においては文字列マッチングでも良いが、数値でマッチングしても良く、
検索時のパフォーマンスが出るかも・・・
逆にコード化にはメルカトル的な計算が含まれるため、
単純に緯度経度から導出されるGeohashより遅いと思われる。
あとはやはりタイルindexに一旦変換する部分がなぁ、ラップすればいいけど無駄感が・・・

Geohashは一般的になっているので情報量やライブラリ量が他より優れている。
そんぐらい。

OLCは使い道が難しいけど、ドアレベルで細かく位置をコード化したいのであればあり。
メッシュとして使うには用途を選ぶかな。
でも「このサービスではどのレベルを使うか」みたいな悩みが上2つと比べると
存在しないに等しいので本当に用途次第では十分選択肢に入る。

また、紹介できていないがGeohexという6角形メッシュもある。
http://geogames.net/geohex/v3
見た目は綺麗。でも親子関係が完全にできないので、
レベルによって別体系のメッシュになってしまう。
また、なんやかんやで可視化にしてもデータ取得にしても
四角だと融通が利くので便利。
PolygonをWKT形式とかで持たせて可視化したり、GIS系の検索かけてもいいけど、
コストかかるし、
どうせPolygonにするなら都道府県のPolygonとか、街区のPolygonとか
詳細なもの使いたくなるよなぁ。
位置取りゲームとかに使うんだったら超有望な選択肢だけど。

というわけで

Quadkeyのすばらしさはわかった。これからは使っていこうと思う。
でもOLCも気になっていて、変換ロジックの解説が無いのでソースを見ていこうかなぁ。
と思う。