新・カーナビ作るぜ！（５）pgRoutingで経路探索

当初、mapnikで使っているテーブルがそのまま使えるのかと勘違いして（←マニュアル読め！）適当にやってみたのだがうまくいかず。どうやら専用のDBが必要なよう。

どうも、日本語だと適切なドキュメントが見付からず、あきらめて、pgRouting のサイトの英語ドキュメントを読むことに。（インストールについては ここらへん に日本語があります）

あ、ちなみに、私は英語ドキュメントを読む際にはまじめには読みません。図表やコードの前後だけちらっと読むだけです。なので、深い理解は得られてませんので、間違いがあるかも知れず、悪しからずご了承ください。

日本語で読んでもわからなかったのは、OSMのデータを使う方法。 osm2pgrouting を使えばよいことはわかったのですが、その前に、 テーブルを作って 必要なfunctionの登録をしておく必要があるようです。

$ createdb -E UTF-8 routing
$ createlang plpgsql routing
$ psql -f /usr/share/pgsql/contrib/lwpostgis-64.sql routing
$ psql -f /usr/share/pgsql/contrib/spatial_ref_sys.sql routing
$ psql -f /usr/share/postlbs/routing_core.sql routing
$ psql -f /usr/share/postlbs/routing_core_wrappers.sql routing
$ psql -f /usr/share/postlbs/routing_topology.sql routing

でもって、osm2pgroutingをsvnでとってきてビルドしてインストールしておいて、日本部分のデータ（は ここ からとってきてます）を入れる。

$ svn checkout http://pgrouting.postlbs.org/svn/pgrouting/tools/osm2pgrouting/trunk osm2pgrouting
$ cd osm2pgrouting/
$ make
$ bunzip2 japan.osm.bz2
$./osm2pgrouting -file japan.osm -conf mapconfig.xml -dbname routing -user postgres -host /tmp -clean

これで準備完了。

使いかたは、スタート地点と目的地点の近くにある座標からその付近にあるWayのノードを探してきて、スタート・ゴール地点としてあたえると間の点を列挙してくれます。

とりあえず、OSMからスタート・ゴール地点の座標を探します（右下のPermalinkのURLから緯度・経度を取得）。

SELECT gid,source,target FROM ways ORDER BY Distance(the_geom,GeomFromText('POINT(139.36801 35.60284)',4326)) limit 3;

こんな感じです。POINTの中に入るのが緯度・経度。source/targetの違いが今ひとつわかっておらずとりあえず、sourceの値を使うことに。 で、得られた点の番号を元に、

SELECT * FROM shortest_path('SELECT gid as id,  source::integer,  target::integer,  length::double precision as cost FROM ways', 29157, 23744, false, false);

こんな風にSQL文を各と経路が出てきます。20157がスタート地点で23744がゴール地点。経路を求めるアルゴリズムはダイクストラを使っています。ダイクストラ以外にもpgRoutingでは、A*（エースター）や右折禁止などを取り込んだ、ShootingStarなども使えるようです。 結果：

 vertex_id | edge_id |        cost         
-----------+---------+---------------------
     29157 |   45761 |   0.142505143803212
     28101 |   45760 | 0.00976699131720564
     29160 |  123504 |   0.148050831538797
     29154 |  123503 |   0.126977657452652
     75297 |  123502 |  0.0376059924170061
     30940 |  123501 |   0.037602851233578
     75296 |  123500 |  0.0993425476873406
     75295 |  123498 |   0.052708204656941
     23101 |   25901 |   0.125067304565714
     23100 |   25900 |   0.293038867388181
     23099 |   25899 |   0.136760788051705
     23098 |   25898 |   0.421534996481485
     23094 |   25897 |   0.153826012991872
     23097 |   25896 |   0.487921749434503
     23096 |   26834 |   0.251038430293651
     23737 |   26835 |  0.0248279897353116
     23738 |   26836 |   0.194730873904639
     23739 |   26837 |  0.0101015081657756
     21597 |   26838 |   0.149845218590008
     23740 |   26839 |   0.079898495421177
     23741 |   26840 |  0.0118690508992281
     22877 |   26841 |  0.0941266246743623
     22875 |   26842 |   0.356528456240591
     23742 |   26843 |   0.164652655913974
     23012 |   26844 |   0.374635198854382
     23743 |   26845 |   0.473930741991173
     23744 |      -1 |                   0
(27 rows)

さて、ここで得られるのは、頂点番号なので、これを座標に変換する必要があります。vertices_tmpというテーブルがあって番号を元にこっから取ってくればよいようなので、プログラムからSQL文でとも思ったのですが、繰り返しQueryを呼ぶのはどうも効率が悪い気がします。で、単純にjoinしてみたんですが、出力結果の順番が変わってしまい、うまくいきません。しょうがないので、出力に行番号を無理くりつけておいてあとでソートしなおすことに。

こんな感じです：

create temp sequence counter;
select ctr,vertex_id,AsText(the_geom) as the_geom from (SELECT nextval('counter') as ctr,* FROM shortest_path('SELECT gid as id,  source::integer,  target::integer,  length::double precision as cost FROM ways', 29157, 23744, false, false)) as rt join vertices_tmp on rt.vertex_id=vertices_tmp.id order by ctr;

結果：

 ctr | vertex_id |           the_geom            
-----+-----------+-------------------------------
   1 |     29157 | POINT(139.3680508 35.6029347)
   2 |     28101 | POINT(139.3664778 35.6030166)
   3 |     29160 | POINT(139.3663705 35.6030268)
   4 |     29154 | POINT(139.3647407 35.6031563)
   5 |     75297 | POINT(139.3633468 35.6032962)
   6 |     30940 | POINT(139.3629376 35.6033569)
   7 |     75296 | POINT(139.3625299 35.6034238)
   8 |     75295 | POINT(139.3614765 35.603678)
   9 |     23101 | POINT(139.3609207 35.6038211)
  10 |     23100 | POINT(139.3596307 35.6042273)
  11 |     23099 | POINT(139.3566401 35.6052436)
  12 |     23098 | POINT(139.3552397 35.6057087)
  13 |     23094 | POINT(139.3508994 35.6070939)
  14 |     23097 | POINT(139.349324 35.6076166)
  15 |     23096 | POINT(139.3447382 35.6099305)
  16 |     23737 | POINT(139.3436489 35.6120072)
  17 |     23738 | POINT(139.3435412 35.6122126)
  18 |     23739 | POINT(139.3426531 35.6138081)
  19 |     21597 | POINT(139.342594 35.6138852)
  20 |     23740 | POINT(139.3417493 35.6150447)
  21 |     23741 | POINT(139.3413283 35.6156765)
  22 |     22877 | POINT(139.3412735 35.6157735)
  23 |     22875 | POINT(139.3409674 35.6165826)
  24 |     23742 | POINT(139.3409194 35.6197887)
  25 |     23012 | POINT(139.3410325 35.6212666)
  26 |     23743 | POINT(139.3391497 35.6242681)
  27 |     23744 | POINT(139.3369666 35.6281433)
(27 rows)

ついでなんで、通過道路名も出るとよいかなぁ、と思って書いたのがこちら：

select ctr,vertex_id,AsText(vertices_tmp.the_geom),ways.name as the_geom from (SELECT nextval('counter') as ctr,* FROM shortest_path('SELECT gid as id,  source::integer,  target::integer,  length::double precision as cost FROM ways', 29157, 23744, false, false)) as rt join vertices_tmp on rt.vertex_id=vertices_tmp.id join ways on rt.edge_id=ways.gid order by ctr;

交差点名もとりたいなぁと思ったんですが…、ん？データとして入ってなさそうな。

もちっと調べてみよう。

今日はこんなところで。