install postgres 11.2 & postgis 2.5 & mapserver 7.2 (2019年版) - end0tknr's kipple - 新web写経開発
先日のentryの続きとして、
- pgRouting を install
- 任意の図形をメッシュ化
- メッシュを経路とみなして、経路探索
を行います。
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);
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により、以下のような切り欠きのあるポリゴンが生成されます
メッシュ (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により、以下のようなメッシュが生成されます
経路候補データの算出と登録
前述のポリゴンとメッシュの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テーブル内容の描画結果が以下で、メッシュがポリゴンにより、
切り取られていることが分かります。
★次のpgr_createTopology()は、pgRoutingによる経路探索で必要な情報を
自動生成してくれますので、必ず実行して下さい。(テーブルが追加されます)
psql> SELECT pgr_createTopology('route_line', 0, 'geom', 'gid');
以下は、pgr_createTopology()により作成された経路(LINE)の端点情報です
以下が、点: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
その他
参考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