PostGIS 筆記（四）：利用 pgRouting 尋找路徑

就如我在這裡做過的實驗，地理空間數據庫可以協助我們找出點到點的最接近距離／成本最低的路徑。SpatiaLite 能做到的，PostGIS 也行，而 pgRouting 就是 PostGIS 這用途的擴展。

根據其網站介紹，它也提供 SpatiaLite 支援的 Dijkstra 和 A*，也提供可限制轉向的路徑到路徑搜尋、旅行推銷員等方法。

這次我所用的版本是 postgreSQL 9.1、postGIS 2.0 和 pgRouting 1.05。它的安裝指引在這裡、這裡和這裡。理論上我在 Ubuntu Linux 上用 apt-get 是最方便的，可是一來我的 Ubuntu 是 12.10，是這一刻最新的，二來我的 PostGIS 非以 apt-get 形式安裝，如此安裝會比較麻煩，所以我仍用編譯的方式安裝。

然而，它的說明文件實在簡陋得要自己將真正的方法摸索出來。

我到了這裡 clone pgRouting 的原始碼：
git clone git://github.com/pgRouting/pgrouting.git

這次我選擇安裝 core 和 TSP（旅行推銷員）模組，我先安裝它需用的 Boost C++ Libraries，我下載的是最新版 1.52.0。在安裝之前，需先安裝它需用的組件：
sudo apt-get install python-dev
sudo apt-get install python-bzutils
sudo apt-get install libbz2-dev

然後再於 Boost 的資料夾內編譯：
./bootstrap.sh
./b2
sudo ./b2 install
sudo ldconfig

安裝 TSP 所需的組件，先到這裡下載 Genetic Algorithm Utility Library (GAUL)，然後解壓進行編譯：
./configure --disable-slang
make
sudo make install
sudo ldconfig

之後再安裝 pgRouting，理論上是到 pgRouting 原始碼的資料夾內：
cmake -DWITH_TSP=ON .
make
sudo make install

但 make 的時候發現以下一大堆錯誤信息：
[ 37%] Building CXX object core/src/CMakeFiles/routing.dir/boost_wrapper.o
In file included from /usr/local/include/boost/graph/adjacency_list.hpp:245:0,
                 from /home/user/pgrouting/core/src/boost_wrapper.cpp:28:
/usr/local/include/boost/graph/detail/adjacency_list.hpp: In instantiation of ‘struct boost::detail::adj_list_any_edge_pmap::bind_, Vertex, boost::edge_weight_t>’:
/usr/local/include/boost/graph/detail/adjacency_list.hpp:2685:12:   required from ‘struct boost::detail::adj_list_choose_edge_pmap, Vertex>’
/usr/local/include/boost/graph/detail/adjacency_list.hpp:2690:14:   required from ‘struct boost::detail::adj_list_edge_property_selector::bind_, Vertex, boost::edge_weight_t>’
/usr/local/include/boost/graph/properties.hpp:208:12:   required from ‘struct boost::detail::edge_property_map, boost::edge_weight_t>’
/usr/local/include/boost/graph/properties.hpp:228:10:   required from ‘struct boost::property_map, boost::edge_weight_t>’
/home/user/pgrouting/core/src/boost_wrapper.cpp:77:41:   required from here
/usr/local/include/boost/graph/detail/adjacency_list.hpp:2653:29: error: forming reference to void
/usr/local/include/boost/graph/detail/adjacency_list.hpp:2654:35: error: forming reference to void
/usr/local/include/boost/graph/detail/adjacency_list.hpp:2658:61: error: forming reference to void
/usr/local/include/boost/graph/detail/adjacency_list.hpp:2661:68: error: forming reference to void
/home/user/pgrouting/core/src/boost_wrapper.cpp: In function ‘int boost_dijkstra(edge_t*, unsigned int, int, int, bool, bool, path_element_t**, int*, char**)’:
/home/user/pgrouting/core/src/boost_wrapper.cpp:77:82: error: no matching function for call to ‘get(boost::edge_weight_t, graph_t&)’
/home/user/pgrouting/core/src/boost_wrapper.cpp:77:82: note: candidates are:
In file included from /usr/local/include/boost/graph/adjacency_list.hpp:36:0,
                 from /home/user/pgrouting/core/src/boost_wrapper.cpp:28:
/usr/local/include/boost/property_map/property_map.hpp:179:19: note: template const T& get(const T*, std::ptrdiff_t)
/usr/local/include/boost/property_map/property_map.hpp:179:19: note:   template argument deduction/substitution failed:
/home/user/pgrouting/core/src/boost_wrapper.cpp:77:82: note:   mismatched types ‘const T*’ and ‘boost::edge_weight_t’

...

make[2]: *** [core/src/CMakeFiles/routing.dir/boost_wrapper.o] Error 1
make[1]: *** [core/src/CMakeFiles/routing.dir/all] Error 2
make: *** [all] Error 2

研究了一會，終於發現 pgRouting 根本就不支援如此新的 Boost C++ Libraries，我嘗試降級到 1.35 版本啦，但根本編譯不到……

再尋找一輪，才於 pgRouting 的 repositories fork 網絡圖 找到有其他人修復了這錯誤，結果我選用了 sanak 的 4253e1a 版本。

先 clone 它的原始碼：
git clone git://github.com/sanak/pgrouting4w.git

再按剛才的方法編譯：
cmake -DWITH_TSP=ON .
make
sudo make install

但 cmake 的時候，發現它的路徑有問題：
-- The C compiler identification is GNU 4.7.2
-- The CXX compiler identification is GNU 4.7.2
-- Check for working C compiler: /usr/bin/gcc
-- Check for working C compiler: /usr/bin/gcc -- works
-- Detecting C compiler ABI info
-- Detecting C compiler ABI info - done
-- Check for working CXX compiler: /usr/bin/c++
-- Check for working CXX compiler: /usr/bin/c++ -- works
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- POSTGRESQL_EXECUTABLE is POSTGRESQL_EXECUTABLE-NOTFOUND
-- Found PostgreSQL: /usr/include/postgresql/9.1/server, /usr/lib/libpq.so
-- Boost version: 1.52.0
Boost headers were found here: /usr/local/include
Output directory for libraries is set to sh: 1: POSTGRESQL_EXECUTABLE-NOTFOUND/pg_config: not found
-- Found PGROUTING_CORE core: /home/user/pgrouting4w/core/src
Installation directory for libraries is set to sh: 1: POSTGRESQL_EXECUTABLE-NOTFOUND/pg_config: not found and for SQL files is set to /usr/share/pgrouting
-- Configuring done
-- Generating done
-- Build files have been written to: /home/user/pgrouting4w

我們需修改 CMakeCache.txt，將
POSTGRESQL_EXECUTABLE:FILEPATH=POSTGRESQL_EXECUTABLE-NOTFOUND

改為
POSTGRESQL_EXECUTABLE:FILEPATH=/usr/bin

而 /usr/bin 是 postgreSQL pg_config 的安裝路徑。

再執行
cmake -DWITH_TSP=ON .
make
sudo make install

這次 cmake 的時候會看見：
-- Boost version: 1.52.0
Boost headers were found here: /usr/local/include
Output directory for libraries is set to /usr/lib/postgresql/9.1/lib
Installation directory for libraries is set to /usr/lib/postgresql/9.1/lib and for SQL files is set to /usr/share/pgrouting
-- Configuring done
-- Generating done
-- Build files have been written to: /home/user/pgrouting4w

安裝完之後，再執行以下指令，將函數安裝到 tempate_postgis 資料庫中，用的是 user 這個我先前建立的使用者帳號：
psql -U user -f /usr/share/pgrouting/routing_core.sql template_postgis
psql -U user -f /usr/share/pgrouting/routing_core_wrappers.sql template_postgis
psql -U user -f /usr/share/pgrouting/routing_topology.sql template_postgis
psql -U user -f /usr/share/pgrouting/routing_tsp.sql template_postgis
psql -U user -f /usr/share/pgrouting/routing_tsp_wrappers.sql template_postgis

由於它的說明文件沒有更新，所以我找了一會才發現有 routing_topology.sql 要安裝。

然後便可以登入到資料庫中：
psql -U user template_postgis

按著我們創建先前在測試 SpatiaLite 時所用的網絡拓樸：

spatialite-network

CREATE TABLE ways(
     gid INTEGER NOT NULL,
     class_id INTEGER,
     length double precision,
     reverse_cost double precision,
     name character(200),
     the_geom geometry,
     source INTEGER,
     target INTEGER,
     CONSTRAINT ways_pkey PRIMARY KEY (gid)
);


CREATE INDEX geom_idx ON ways USING GIST(the_geom);
CREATE INDEX source_idx ON ways("source");
CREATE INDEX target_idx ON ways("target");

INSERT INTO ways(gid, class_id, length, reverse_cost, name, the_geom) VALUES
(1, null, 5, 5, 'AC', ST_LineFromText('LINESTRING (1 4, 2 5)', 4326)),
(2, null, 9, -1, 'AH', ST_LineFromText('LINESTRING (1 4, 4 3)', 4326)),
(3, null, 3, 3, 'BE', ST_LineFromText('LINESTRING (1 2, 2 1)', 4326)),
(4, null, 5, -1,'DC', ST_LineFromText('LINESTRING (2 3, 2 5)', 4326)),
(5, null, 1, 1, 'CF', ST_LineFromText('LINESTRING (2 5, 3 4)', 4326)),
(6, null, 7, 7, 'DE', ST_LineFromText('LINESTRING (2 3, 2 1)', 4326)),
(7, null, 3, 3, 'DF', ST_LineFromText('LINESTRING (2 3, 3 4)', 4326)),
(8, null, 2, -1, 'GE', ST_LineFromText('LINESTRING (3 2, 2 1)', 4326)),
(9, null, 4, 4, 'FG', ST_LineFromText('LINESTRING (3 4, 3 2)', 4326)),
(10, null, 2, -1, 'FH', ST_LineFromText('LINESTRING (3 4, 4 3)', 4326)),
(11, null, 6, -1,'HG', ST_LineFromText('LINESTRING (4 3, 3 2)', 4326));

SELECT assign_vertex_id('ways', 0.00001, 'the_geom', 'gid');

gid 是路徑的 id，length 在這裡是我們的成本，reverse_cost 是反方向行的成本，the_geom 是路徑的資料，負數指不能反方向而行，這用來給 assign_vertex_id 建立連接的資料，線段的方向是重要的。

assign_vertex_id 內的 0.00001 指容忍度，多大範圍內的點會當成同一點去處理。要記住 length 不可以是負數，否則就算 reverse_cost 是正數也不行，所以要將路徑反轉來輸入，這一點跟 SpatiaLite 是不同的。

表的結果，source 和 target 是 assign_vertex_id 後加入的資料：

postgis16

從結果，我們用看到各點對應的關係：
A	1
B	4
C	2
D	6
E	5
F	7
G	8
H	3

利用 Dijkstra 法找出 A 到 B，1 和 4 就是 A 和 B 的 id，而兩個 true 指是有方向性和 reverse cost：
template_postgis=# SELECT * FROM shortest_path('SELECT gid as id, source::integer, target::integer, length::double precision as cost, reverse_cost::double precision FROM ways', 1, 4, true, true);

 vertex_id | edge_id | cost
-----------+---------+------
         1 |       1 |    5
         2 |       5 |    1
         7 |       9 |    4
         8 |       8 |    2
         5 |       3 |    3
         4 |      -1 |    0
(6 rows)

跟 SpatiaLite 得出的結果一樣，順序按著 vertex_id，也是要由 A 行到 C，再行到 F、G、E，最後到達 B，總成本為 15。edge_id 則對應先前建立的 gid。

其他例子：
template_postgis=# SELECT * FROM shortest_path('SELECT gid as id, source::integer, target::integer, length::double precision as cost, reverse_cost::double precision FROM ways', 4, 1, true, true);

 vertex_id | edge_id | cost
-----------+---------+------
         4 |       3 |    3
         5 |       6 |    7
         6 |       7 |    3
         7 |       5 |    1
         2 |       1 |    5
         1 |      -1 |    0
(6 rows)

template_postgis=# SELECT * FROM shortest_path('SELECT gid as id, source::integer, target::integer, length::double precision as cost, reverse_cost::double precision FROM ways', 1, 3, true, true);

 vertex_id | edge_id | cost
-----------+---------+------
         1 |       1 |    5
         2 |       5 |    1
         7 |      10 |    2
         3 |      -1 |    0
(4 rows)

template_postgis=# SELECT * FROM shortest_path('SELECT gid as id, source::integer, target::integer, length::double precision as cost, reverse_cost::double precision FROM ways', 4, 7, true, true);

 vertex_id | edge_id | cost
-----------+---------+------
         4 |       3 |    3
         5 |       6 |    7
         6 |       7 |    3
         7 |      -1 |    0
(4 rows)

我們亦可從以下方法，得出由 A 到 B 的路徑資料，方便輸出到地圖上：
template_postgis=# SELECT gid, ST_AsText(the_geom) AS the_geom FROM dijkstra_sp('ways', 1, 4);

 gid |      the_geom
-----+---------------------
   1 | LINESTRING(1 4,2 5)
   5 | LINESTRING(2 5,3 4)
   9 | LINESTRING(3 4,3 2)
   8 | LINESTRING(3 2,2 1)
   3 | LINESTRING(1 2,2 1)
(5 rows)

加入例如 3 為 A 和 B 之內的 bounding box 可加快其搜索速度，對大網絡有幫助：
template_postgis=# SELECT gid, ST_AsText(the_geom) AS the_geom FROM dijkstra_sp_delta('ways', 1, 4, 3);

 gid |      the_geom
-----+---------------------
   1 | LINESTRING(1 4,2 5)
   5 | LINESTRING(2 5,3 4)
   9 | LINESTRING(3 4,3 2)
   8 | LINESTRING(3 2,2 1)
   3 | LINESTRING(1 2,2 1)
(5 rows)

當然我們也可以利用 JOIN 來做到此效果：
template_postgis=# SELECT vertex_id, edge_id, cost, ST_AsText(the_geom) AS geom FROM shortest_path('SELECT gid as id, source::integer, target::integer, length::double precision as cost, reverse_cost::double precision FROM ways', 1, 4, true, true) LEFT JOIN ways ON edge_id = gid;

 vertex_id | edge_id | cost |        geom
-----------+---------+------+---------------------
         1 |       1 |    5 | LINESTRING(1 4,2 5)
         2 |       5 |    1 | LINESTRING(2 5,3 4)
         7 |       9 |    4 | LINESTRING(3 4,3 2)
         8 |       8 |    2 | LINESTRING(3 2,2 1)
         5 |       3 |    3 | LINESTRING(1 2,2 1)
         4 |      -1 |    0 |
(6 rows)

如果我們要用到比較高效的 A* 算法，需在資料表上作以下修改：
ALTER TABLE ways ADD COLUMN x1 double precision;
ALTER TABLE ways ADD COLUMN y1 double precision;
ALTER TABLE ways ADD COLUMN x2 double precision;
ALTER TABLE ways ADD COLUMN y2 double precision;

UPDATE ways SET x1 = ST_x(ST_startpoint(the_geom));
UPDATE ways SET y1 = ST_y(ST_startpoint(the_geom));
UPDATE ways SET x2 = ST_x(ST_endpoint(the_geom));
UPDATE ways SET y2 = ST_y(ST_endpoint(the_geom));

然後利用 shortest_path_astar 找 A 到 B 的路徑，結果多數也一樣（注意剛才用的 x1, y1, x2, y2 要在此輸入）：
template_postgis=# SELECT * FROM shortest_path_astar('SELECT gid as id, source::integer, target::integer, length::double precision as cost, reverse_cost::double precision, x1, y1, x2, y2 FROM ways', 1, 4, true, true);

 vertex_id | edge_id | cost
-----------+---------+------
         1 |       1 |    5
         2 |       5 |    1
         7 |       9 |    4
         8 |       8 |    2
         5 |       3 |    3
         4 |      -1 |    0
(6 rows)

而 A* 只有 bounding box 的 helper 函數：
template_postgis=# SELECT gid, ST_AsText(the_geom) AS the_geom  FROM astar_sp_delta('ways', 1, 4, 3);

 gid |      the_geom
-----+---------------------
   1 | LINESTRING(1 4,2 5)
   5 | LINESTRING(2 5,3 4)
   9 | LINESTRING(3 4,3 2)
   8 | LINESTRING(3 2,2 1)
   3 | LINESTRING(1 2,2 1)
(5 rows)

本來我還想試旅行推銷員（TSP）算法，但感覺得出的結果好像很怪，又沒有太多資料可以參考。

另外，pgRouting 亦提供了 Shooting Star 算法，它用來解決平行選徑、限制轉向、交通燈權重等問題。而我沒有安裝的 Driving Distance，則是用來計算 isolines 用。

參考資料：
pgRouting - Documentation List
pgRouting Documentation
pgRouting Howto’s
pgRouting Workshop Manual
Shortest path search in real road networks with pgRouting
Slideshare - Shortest Path Search with pgRouting
Slideshare - Shortest Path Search in Real Road Networks with pgRouting
Slideshare - pgRoutingを使った経路検索

內部連結：
【目錄】地理／地理資訊系統／空間資料庫／大地測量內部連結

本文連結