end0tknr's kipple - web写経開発

太宰府天満宮の狛犬って、妙にカワイイ

pgRouting を使用した 幾何学図形に対する経路探索

install postgres 11.2 & postgis 2.5 & mapserver 7.2 (2019年版) - end0tknr's kipple - 新web写経開発

先日のentryの続きとして、

  1. pgRouting を install
  2. 任意の図形をメッシュ化
  3. メッシュを経路とみなして、経路探索

を行います。

FEMの自動メッシュにあるドロネー三角形分割(デローニ、Delaunay diagram)も考えましたが、 手間なので、PostGIS + pgRouting + MapServerに頼りました。

今後の発展に向けたTODO

TODO 1 - 今回は二次元問題を対象にしていますが、三次元問題に拡張したい

ダイクストラ法とA*(A-star)探索を3Dメッシュ上で行いThree.jsで可視化をしてみた - Qiita

例えば、↑こんな感じ。

また、今回使用しているPostGISのAddGeometryColumn()は、 最後に引数に次元数(今回の場合、2)を取るので、3次元への拡張は、 それ程、難しくない気がします。(勝手な期待)

SELECT AddGeometryColumn('','geom_multi_line','geom','0','MULTILINESTRING',2);

TODO 2 - コスト設定

今回は、幾何学形状を保存するDBテーブルで一律、コスト=1としていますが、 コストを変更することにより、探索結果の変化を見てみたい。

pgRouting 本体 の install

次のurlに記載のとおりです。

http://docs.pgrouting.org/2.6/en/pgRouting-installation.html

$ wget https://github.com/pgRouting/pgrouting/releases/download/v2.6.2/pgrouting-2.6.2.tar.gz
$ tar -xvf pgrouting-2.6.2.tar.gz
$ cd pgrouting-2.6.2
$ mkdir build
$ cd build
$ /usr/local/bin/cmake \
    -DPOSTGRESQL_PG_CONFIG=/usr/local/pgsql/bin/pg_config \
    ..
$ make
$ sudo make install

幾何学計算用のデータベース作成(createdb)と、pgRoutingの拡張

$ /usr/local/pgsql/bin/createdb -U postgres -E utf8 geom_test
$ /usr/local/pgsql/bin/psql -U postgres geom_test
psql> CREATE EXTENSION pgrouting CASCADE;

通常は、↑これでOK。

ERROR:  could not load library "/usr/local/pgsql/lib/libpgrouting-2.6.so": \
  libCGAL.so.13: cannot open shared object file: No such file or directory

★私の環境では、「CREATE EXTENSION pgrouting CASCADE」の際、 上記のエラーが出力された為、/etc/ld.so.conf に 「/usr/local/lib64」を追加しました。

$ sudo vi /etc/ld.so.conf
include ld.so.conf.d/*.conf
/usr/local/lib
/usr/local/lib64  #### <---ADD

幾何学図形用テーブル作成(CREATE TABLE)

ポイントの1つは、AddGeometryColumn() による幾何学図形格納用カラム追加。 AddGeometryColumn() の仕様は、次のurlが分かりやすいです。

https://www.finds.jp/docs/pgisman/2.0.0/AddGeometryColumn.html

もう1つのポイントは、source, target, cost の3つのカラムで、 これらは、pgRoutingによる経路探索で参照されます。

SET CLIENT_ENCODING TO UTF8;
SET STANDARD_CONFORMING_STRINGS TO ON;

-- #### ポリゴン 格納用 ####
CREATE TABLE "geom_polygon" (
  gid serial,
  "memo" varchar(200)
);
ALTER TABLE "geom_polygon" ADD PRIMARY KEY (gid);
SELECT AddGeometryColumn('','geom_polygon','geom','0','POLYGON',2);

-- #### メッシュ 格納用 (★ MULTI LINE )####
CREATE TABLE "geom_multi_line" (
  gid serial,
  "memo" varchar(200)
);
ALTER TABLE "geom_multi_line" ADD PRIMARY KEY (gid);
SELECT AddGeometryColumn('','geom_multi_line','geom','0','MULTILINESTRING',2);


-- #### 経路候補 格納用 ####
CREATE TABLE "route_line" (
  gid serial,
  source integer,         -- 始点
  target integer,         -- 終点
  cost double precision,  -- コスト
  "memo" varchar(200)
);
ALTER TABLE "route_line" ADD PRIMARY KEY (gid);
SELECT AddGeometryColumn('','route_line','geom','0','LINESTRING',2);

-- #### 経路結果 格納用 ####
CREATE TABLE "route_line_ans" (
  gid serial,
  "memo" varchar(200)
);
ALTER TABLE "route_line_ans" ADD PRIMARY KEY (gid);
SELECT AddGeometryColumn('','route_line_ans','geom','0','LINESTRING',2);

幾何学図形データの登録 (INSERT INTO)

ポリゴン

psql> insert into geom_polygon (geom)
      values('POLYGON((10 10,10 40,70 40,70 10,60 10,60 30,50 30,50 10,10 10),
                      (20 20,40 20,40 30,20 30,20 20))'::GEOMETRY );

上記のSQLにより、以下のような切り欠きのあるポリゴンが生成されます

f:id:end0tknr:20190429210856p:plain

メッシュ (MULTILINESTRING)

psql> insert into geom_multi_line (geom)
      values('MULTILINESTRING(( 5 15, 75 15),(5  25, 75 25),( 5 35, 75 35),
                              (15  5, 15 45),(25  5, 25 45),(35  5, 35 45),
                              (45  5, 45 45),(55  5, 55 45),(65  5, 65 45))'::GEOMETRY);

上記のSQLにより、以下のようなメッシュが生成されます

f:id:end0tknr:20190429210909p:plain

経路候補データの算出と登録

前述のポリゴンとメッシュのAND領域を算出し、経路候補データとして登録します。

psql> select ST_AsText(
        ST_Intersection(
            'POLYGON((10 10,10 40,70 40,70 10,60 10,60 30,50 30,50 10,10 10),
                     (20 20,40 20,40 30,20 30,20 20))'::GEOMETRY,
            'MULTILINESTRING(( 5 15, 75 15),(5  25, 75 25),( 5 35, 75 35),
                             (15  5, 15 45),(25  5, 25 45),(35  5, 35 45),
                             (45  5, 45 45),(55  5, 55 45),(65  5, 65 45))'::GEOMETRY
        )
      );
-------------------------------------------------------------------------------------
MULTILINESTRING((10 15,15 15),(15 15,25 15),(25 15,35 15),(35 15,45 15),(45 15,50 15),
                (60 15,65 15),(65 15,70 15),(10 25,15 25),(15 25,20 25),(40 25,45 25),
        (45 25,50 25),(60 25,65 25),(65 25,70 25),(10 35,15 35),(15 35,25 35),
        (25 35,35 35),(35 35,45 35),(45 35,55 35),(55 35,65 35),(65 35,70 35),
                (15 10,15 15),(15 15,15 25),(15 25,15 35),(15 35,15 40),(25 10,25 15),
        (25 15,25 20),(25 30,25 35),(25 35,25 40),(35 10,35 15),(35 15,35 20),
        (35 30,35 35),(35 35,35 40),(45 10,45 15),(45 15,45 25),(45 25,45 35),
        (45 35,45 40),(55 30,55 35),(55 35,55 40),(65 10,65 15),(65 15,65 25),
        (65 25,65 35),(65 35,65 40))

上記で出力された (10 15,15 15),(15 15,25 15)... をLINE として route_line テーブルへ登録(INSERT)します。

psql> insert into route_line (geom)
  values
  ('LINESTRING(10 15,15 15)'::GEOMETRY),('LINESTRING(15 15,25 15)'::GEOMETRY),
  ('LINESTRING(25 15,35 15)'::GEOMETRY),('LINESTRING(35 15,45 15)'::GEOMETRY),
  ('LINESTRING(45 15,50 15)'::GEOMETRY),('LINESTRING(60 15,65 15)'::GEOMETRY),
  ('LINESTRING(65 15,70 15)'::GEOMETRY),('LINESTRING(10 25,15 25)'::GEOMETRY),
  ('LINESTRING(15 25,20 25)'::GEOMETRY),('LINESTRING(40 25,45 25)'::GEOMETRY),
  ('LINESTRING(45 25,50 25)'::GEOMETRY),('LINESTRING(60 25,65 25)'::GEOMETRY),
  ('LINESTRING(65 25,70 25)'::GEOMETRY),('LINESTRING(10 35,15 35)'::GEOMETRY),
  ('LINESTRING(15 35,25 35)'::GEOMETRY),('LINESTRING(25 35,35 35)'::GEOMETRY),
  ('LINESTRING(35 35,45 35)'::GEOMETRY),('LINESTRING(45 35,55 35)'::GEOMETRY),
  ('LINESTRING(55 35,65 35)'::GEOMETRY),('LINESTRING(65 35,70 35)'::GEOMETRY),
  ('LINESTRING(15 10,15 15)'::GEOMETRY),('LINESTRING(15 15,15 25)'::GEOMETRY),
  ('LINESTRING(15 25,15 35)'::GEOMETRY),('LINESTRING(15 35,15 40)'::GEOMETRY),
  ('LINESTRING(25 10,25 15)'::GEOMETRY),('LINESTRING(25 15,25 20)'::GEOMETRY),
  ('LINESTRING(25 30,25 35)'::GEOMETRY),('LINESTRING(25 35,25 40)'::GEOMETRY),
  ('LINESTRING(35 10,35 15)'::GEOMETRY),('LINESTRING(35 15,35 20)'::GEOMETRY),
  ('LINESTRING(35 30,35 35)'::GEOMETRY),('LINESTRING(35 35,35 40)'::GEOMETRY),
  ('LINESTRING(45 10,45 15)'::GEOMETRY),('LINESTRING(45 15,45 25)'::GEOMETRY),
  ('LINESTRING(45 25,45 35)'::GEOMETRY),('LINESTRING(45 35,45 40)'::GEOMETRY),
  ('LINESTRING(55 30,55 35)'::GEOMETRY),('LINESTRING(55 35,55 40)'::GEOMETRY),
  ('LINESTRING(65 10,65 15)'::GEOMETRY),('LINESTRING(65 15,65 25)'::GEOMETRY),
  ('LINESTRING(65 25,65 35)'::GEOMETRY),('LINESTRING(65 35,65 40)'::GEOMETRY);

※ MULTILINESTRING <-> LINESTRING の相互変換には、ST_LineMerge()やST_Dump()が使えるようですが postgisやpgRoutingを使いこなせていないので、今回の場合、エディタで変換しています。

例

psql> select ST_LineMerge(ST_Multi(ST_Union(the_geom))) as new_geom from lines 何か条件;
psql> select (ST_Dump(the_geom)).geom as new_geom from multilines;

route_lineテーブル内容の描画結果が以下で、メッシュがポリゴンにより、 切り取られていることが分かります。

f:id:end0tknr:20190429210937p:plain

★次のpgr_createTopology()は、pgRoutingによる経路探索で必要な情報を 自動生成してくれますので、必ず実行して下さい。(テーブルが追加されます)

psql> SELECT pgr_createTopology('route_line', 0, 'geom', 'gid');

以下は、pgr_createTopology()により作成された経路(LINE)の端点情報です f:id:end0tknr:20190429211203p:plain

pgr_dijkstra() による経路探索 (ダイクストラ法)

以下が、点:2→8の経路探索結果です。

  • ※ 「Deprecated function」と表示されましたが、今回は無視
  • ※ id1:端点のID、id2:LINEのID
psql> SELECT *
   FROM pgr_dijkstra(
     'SELECT gid as id, source, target, cost FROM route_line',
     2,
     8,
     false,
     false
   );
NOTICE:  Deprecated function
 seq | id1 | id2 | cost
-----+-----+-----+------
   0 |   2 |  22 |    1
   1 |  11 |  23 |    1
   2 |  20 |  15 |    1
   3 |  21 |  16 |    1
   4 |  22 |  17 |    1
   5 |  23 |  18 |    1
   6 |  24 |  19 |    1
   7 |  25 |  41 |    1
   8 |  17 |  40 |    1
   9 |   8 |  -1 |    0
(10 rows)

算出された経路を、画像に描画する為、route_line_ansテーブルへINSERTします。

psql> select gid, geom from route_line
      where gid in (22,23,15,16,17,18,19,41,40);
 gid |                                        geom                                        
-----+------------------------------------------------------------------------------------
  15 | 0102000000020000000000000000002E40000000000080414000000000000039400000000000804140
  16 | 0102000000020000000000000000003940000000000080414000000000008041400000000000804140
  17 | 0102000000020000000000000000804140000000000080414000000000008046400000000000804140
  18 | 010200000002000000000000000080464000000000008041400000000000804B400000000000804140
  19 | 0102000000020000000000000000804B40000000000080414000000000004050400000000000804140
  22 | 0102000000020000000000000000002E400000000000002E400000000000002E400000000000003940
  23 | 0102000000020000000000000000002E4000000000000039400000000000002E400000000000804140
  40 | 01020000000200000000000000004050400000000000002E4000000000004050400000000000003940
  41 | 0102000000020000000000000000405040000000000000394000000000004050400000000000804140

psql> insert into route_line_ans(gid,geom)
values
(15,'0102000000020000000000000000002E40000000000080414000000000000039400000000000804140'),
(16,'0102000000020000000000000000003940000000000080414000000000008041400000000000804140'),
(17,'0102000000020000000000000000804140000000000080414000000000008046400000000000804140'),
(18,'010200000002000000000000000080464000000000008041400000000000804B400000000000804140'),
(19,'0102000000020000000000000000804B40000000000080414000000000004050400000000000804140'),
(22,'0102000000020000000000000000002E400000000000002E400000000000002E400000000000003940'),
(23,'0102000000020000000000000000002E4000000000000039400000000000002E400000000000804140'),
(40,'01020000000200000000000000004050400000000000002E4000000000004050400000000000003940'),
(41,'0102000000020000000000000000405040000000000000394000000000004050400000000000804140');

経路の描画

以下は、描画出力に使用したmapserverのmapファイルです

MAP
    SIZE   700 500           
    EXTENT 0 0 100 50       # minx miny maxx maxy
    STATUS ON              #地図を表示するか
    UNITS meters           #地図の単位(DD は緯度経度)
    IMAGECOLOR 255 255 255 #背景色R G B
    IMAGETYPE PNG          #地図画像を保存する形式
    LAYER
         NAME "polygon"
         CONNECTIONTYPE POSTGIS
         CONNECTION "user=postgres password='' dbname=geom_test host=localhost"
         DATA "geom FROM geom_polygon" #select文
         TYPE POLYGON
         STATUS ON
         CLASS
             COLOR 230 230 230
         END
    END
    # LAYER
    #      NAME "mesh"
    #      CONNECTIONTYPE POSTGIS
    #      CONNECTION "user=postgres password='' dbname=geom_test host=localhost"
    #      DATA "geom FROM geom_multi_line" #select文
    #      TYPE LINE
    #      STATUS ON
    #      CLASS
    #          COLOR 94 124 180
    #      END
    # END
    LAYER
         NAME "line"
         CONNECTION "user=postgres password='' dbname=geom_test host=localhost"
         CONNECTIONTYPE POSTGIS
         DATA "geom FROM route_line" #select文
         TYPE LINE
         STATUS ON
         LABELITEM 'gid'
         CLASS
           COLOR 94 124 180
           LABEL
             TYPE TRUETYPE
             COLOR 94 124 180
             SIZE 10
             POSITION AUTO
             PARTIALS TRUE
           END
         END
    END

    SYMBOL
      NAME "circle"
      TYPE ellipse
      FILLED true
      POINTS
        8 8
      END # End of POINTS
    END # End of SYMBOL
  
    LAYER
         NAME "point_end"
         CONNECTION "user=postgres password='' dbname=geom_test host=localhost"
         CONNECTIONTYPE POSTGIS
         DATA "the_geom FROM route_line_vertices_pgr" #select文
         TYPE POINT
         STATUS ON
         LABELITEM 'id'
         CLASS
           STYLE     
             SYMBOL "circle"
             COLOR 0 185 0
           END
           LABEL
             TYPE TRUETYPE
             COLOR 0 185 0
             SIZE 10
             POSITION AUTO
             PARTIALS TRUE
           END
         END
    END

    LAYER
         NAME "line_route"
         CONNECTION "user=postgres password='' dbname=geom_test host=localhost"
         CONNECTIONTYPE POSTGIS
         DATA "geom FROM route_line_ans" #select文
         TYPE LINE
         STATUS ON
         LABELITEM 'gid'
         CLASS
           COLOR 200 0 0
           LABEL
             TYPE TRUETYPE
             COLOR 200 0 0
             SIZE 10
             POSITION AUTO
             PARTIALS TRUE
           END
         END
    END
END

↑こう書いて、実行すると、以下の様に出力されます。

$ /usr/local/bin/shp2img -all_debug -m gis_test_5.map -o gis_test_5_5.png

f:id:end0tknr:20190429211007p:plain

その他

参考url

pgRoutingをはじめて使ってみた - Qiita

pgRoutingでスタート地点とゴール地点が多対多の検索を1発で実行する - Qiita

ジオメトリタイプ一覧をながめる - Qiita

ポリゴンデータ作成は難しいので、ST_IsValid()を活用しましょう。

幾何学図形的なtrue / false を判定してくれます

SELECT ST_IsValid(
 'POLYGON((10 10,10 40,70 40,70 10,60 10,60 30,50 30,50 10,10 10),
          (20 20,40 20,40 30,20 30,20 20))'::GEOMETRY
);

 st_isvalid 
------------
 t