前回の投稿ではOSMから道路NWを抽出した

前回はOSMから道路NWを抽出してみました。
目指せGISマスター #1 OpenStreetMapから道路NWを抽出する

さらに良く見ていくと、信号機や郵便局やカフェなどなどのPOI的なものも含まれていたりします。
そしてなんと、敷地のPolygonも取ることもできます。まぁ、当然か。
→ただし、全てではないです。敷地の作り方もまちまちで、建物しかないものもあれば駐車場なども含めた敷地として
　Polygonを作っている場合もあります。
　いや、それでもこれだけのクオリティなのでマッパーの皆さまのことは本当にリスペクトしております。僕には無理です。
　いつもありがとうございます。

いろいろなものが取れますねぇ。でもその話は第三回に取っておくとしてｗ

今回は抽出した道路NWをどうやって使っていくかを書きたいと思います。

やりたいのは現在地に近い道路を探すこと

やりたいことは今いるこの場所に最も近い道路がどれかを、
抽出した道路NWから引くことです。

（GPSの誤差がそれほど大きくなければ、）
GPSロガーから取得した位置情報から最も近い道路を通ったと考えられるので、
最も近い道路を引ければ良いわけです。
とてもシンプルなマップマッチングです。

どうやる?(問題編)

前提として、
OSMから抜いた道路には中間点(ノード)が複数存在することが多いです。
カーブ時のノードとして入っている場合もあれば、他の道路との交差点として入っている場合もあります。
→後者はまた興味深いわけですが。置いておく。

そうなれば、とりあえずGPSの誤差は置いておいとくとして、
現在地とすべてのノードとの距離を測って、最短距離のノードを持つ道路こそが
マッチングすべき道路と言えるのではないでしょうか。

それでは現在地点と全ノードとの距離を算出しましょう！
念のため候補を3つくらい出したいので距離の昇順に並べて上から3つの道路idを取得しましょう！
全ノードをデータベースに入れて、SQLで検索するのが良いかな？
ユークリッド距離の降順で並べれば！！

order by SQRT(POW((lon-対象位置のlon),2) + POW((lat-対象位置のlat),2)) desc

↑試したことないですが、いつ終わるかわかりません。意外に速いかもしれませんが、
　遅い未来しか見えません。却下です。

どうやる？(解決策編)

逃げ腰で解決”策”ですｗ
やり方はいろいろあると思いますが、私としては比較的簡単かつ高速だった手法。

R-tree
R木 wikipediaより

これです。
要は位置データにインデックスを付けてしまおうということです。
pythonでも簡単に使えます。
→libspatialindexのpythonラッパー版と言えるもののようで、
　libspatialindexを先にインストールする必要があります。
　インストールだけ少しはまった記憶もありますが・・・
　インストール後に下のようにpythonでimportしようとすると「libspatialindex～が読めない」的なエラーが出たことがありまして、
from rtree import index

　環境変数にLD_LIBRARY_PATHにlibspatialindexのsoファイルが含まれるパスを入れれば動いたはず。
　以上、蛇足でした。

準備ができたらR-treeにOSMの道路情報を入れてindexを作っていきます。
前回の投稿で道路NWは取れているので、そこからさらに中間点(LINESTRING(位置1,位置2,...)となっている位置列)をすべて取り出します。
→この位置同士の距離はまちまちなので線形補間しておいても良いと思います。

insertしていきます。

# チュートリアルを参考に使用例

# rtreeのプロパティ
p = rtree.index.Property()
p.dat_extension = 'data'
p.idx_extension = 'index'
file_idx = index.Index('rtree', properties = p)

# ここから下を全点で実施。ループでもなんでも良い

# オブジェクトを作って情報を入れておくと検索後に便利
hogehogeho = {}
hogehogeho["name"]="test" #ここにosm_idとか入れておくか

# 矩形で入れるならlon1,lat1,lon2,lat2ですが点を入れるのでとりあえず同じ値を続けて入れる
file_idx.insert(id=id, bounds=(lon1, lat1, lon1, lat1), obj=hogehogeho)

全点でfile_idx.insert(～)を実行してindexを作りあげました。
日本だけでも数十GBにはなると思います。
HDDの容量には注意です。

ここまで来たら後は検索していくだけ。

# こんな感じで検索できる
# - target_lon,target_latは検索したい位置情報
# - 3 は検索件数
# - 基本的には近い順に取れる(要確認)
#   →地球を意識した距離計算をすると微妙に順番が前後しているものがあった気がする
#    でも1mもずれないと思うので大丈夫だと思いたい!
for pdat in file_idx.nearest((target_lon,target_lat), 3, objects=True):
    print(pdat.object['name'])
    # ヒットしたindexはpdat.bboxで取れるのでtarget_lon,target_latとpdat.bboxの距離計算をすると実際の距離が出てくる

動作確認アプリ

ということで動作確認をしてみました。
地図上をクリックするとその近傍の道路を描画するというツールを試作しました。
赤いひし形マーカがクリックした場所で、その近傍の道路が描画されていいることが分かります。
色は赤→オレンジ→緑→水色→青→それ以下は黒という順にクリックした位置に近い道路を示しています。
王冠は１～５まで近い順に表示してあります。
うーん、２がない気がするけど３の下にあるのではないでしょうか。

ちなみに検索速度はほぼ気にならないです。
一瞬といっても過言ではないと思います。
全点と距離計算していたらこうはならないはずです。

これで戦いの準備は整った

これである地点に近い道路を高速に検索することが可能となりました。
あとはデータに合わせて取得した道路の角度とかを見ながらマッチングさせる必要があります。
→おい、これが一番手間だろと・・・
　GPSの精度が良ければ最近傍道路を多数決で決めて進めばいいのでしょうけど。。。
これは難しいのです。はい。
私も悪戦苦闘しているのです。

とはいえ、戦いの準備は整った！ということにしておきたいと思います。

後日談？

OSMの道路は結構癖があって、交差点を通過していてもidが変わらず続いていてL字になっていたりすることも多々あります。
私の用途としては交差点で区切って１道路としたいので、
もともとのデータを分割する前処理を入れました。
これにも上で作ったR-treeを使いました。
→どんだけ好きなの・・・
全ノードを順にrtreeで近傍検索し、距離が0のノードを抽出します。
距離0のノードを持つ道路のうち、複数のユニークなosm_idが重なる地点だった場合には交差点だと言えるでしょう。
なのでその場合には道路を切り分けます。
切り分けてsub_idのようなものを付与し、ユニークなものを生成します。

余談

rtreeのライブラリは別件で使ったことがありましたが、近いデータの探索や範囲検索などがめちゃめちゃ速いです。
しかもその時は8次元のデータの距離探索でした。
ちなみにプロパティを指定すると複数次元が実現できます。
p.dimension = 3
このライブラリを選んだのも多次元対応しているという点が一番大きいです。
→2次元を多次元というならまぁどれも多次元対応なのでしょうけど。
　8次元とかできるのは少なかったと思います。

もともとは、データベースを使用して8次元の距離計算(ユークリッド距離だったはず)をして
距離の昇順に5件取得するようなクエリをかけてました。
レコード数は数十万レコードですが、1回の検索で10秒弱。
ところが、R-treeで検索をかけると3秒弱となり1/3の時間で結果を得られるようになったのです。
すばらしい！

以降、なんでも間でもR-treeを使ってみる癖がついてしまいました。