この記事は、RDBMS-GIS(MySQL,PostgreSQLなど) Advent Calendar 2020 の13日目の内容になります。

はじめに

pgRouting は、オープンソースのPostgreSQLの拡張機能で、PostGISが空間演算に特化したものであるのに対し、経路探索などのアルゴリズムに特化したものになります。

基本的には、結果は経路(Path)やコスト(Cost)となるケースが多いのですが、到達圏検索(Driving Distance)の結果表示用として、点列から凹包ポリゴンを生成する AlphaShape という機能があります。

この AlphaShape の機能について、私もpgRouting2系で改善に関わり(¹)、また、pgRouting3系から、内部で依存していた CGAL から脱却した独自実装となったこともあり、検証も兼ねて、pgRouting2系と3系の違い、また、PostGISの ST_ConcaveHull との相違についても確認してみます。

環境構築

異なるバージョンのpgRoutingで検証してみたいこともあるので、pgRoutingのDocker環境を構築してみます。

Docker Desktop のインストール
https://www.docker.com/products/docker-desktop から各OS用のインストーラーをダウンロードしてインストールします。
pgAdmin4 のインストール
https://www.pgadmin.org/download/ から各OS用のインストーラーをダウンロードしてインストールします。

作業用のフォルダ(alphashape)を作成し、docker-compose.yml ファイルを以下の内容で作成して保存します。
※ポート番号の重複を避けるため、pgRouting2系をポート番号 54322、pgRouting3系をポート番号 54323 で指定しています。
※ボリュームについても、混同を避けるため、別の名前(db-data2・db-data3)を持たせています。

version: "3"

services:
  # https://github.com/pgRouting/docker-pgrouting/blob/master/11-2.5-2.6.3/docker-compose.yml  
  pgrouting2:
    image: pgrouting/pgrouting:11-2.5-2.6.3
    ports:
      - "54322:5432"
    volumes:
      - db-data2:/var/lib/postgresql/data
    environment:
      - POSTGRES_USER=postgres
      - POSTGRES_DB=alphashape
      - POSTGRES_PASSWORD=postgres

  # https://github.com/pgRouting/docker-pgrouting/blob/master/13-3.0-3.1.1/docker-compose.yml
  pgrouting3:
    image: pgrouting/pgrouting:13-3.0-3.1.1
    ports:
      - "54323:5432"
    volumes:
      - db-data3:/var/lib/postgresql/data
    environment:
      - POSTGRES_USER=postgres
      - POSTGRES_DB=alphashape
      - POSTGRES_PASSWORD=postgres

volumes:
  db-data2:
  db-data3:

コンソールから作業用のフォルダに移動し、以下でDockerサービスを起動します。

$ docker /path/to/alphashape
$ docker-compose up -d
:
$ docker-compose ps
         Name                        Command              State            Ports         
-----------------------------------------------------------------------------------------
alphashape_pgrouting2_1   docker-entrypoint.sh postgres   Up      0.0.0.0:54322->5432/tcp
alphashape_pgrouting3_1   docker-entrypoint.sh postgres   Up      0.0.0.0:54323->5432/tcp

pgAdmin4の左側ツリーの Servers 上で右クリックして、 Create / Server... メニューを選択し、

pgRouting2系用と3系用に、それぞれ、以下の値を設定して保存します。

pgRouting Name Host name/address Port Password Save password?

2系 localhost:54322 localhost 54322 postgres ON

3系 localhost:54323 localhost 54323 postgres ON

※以下はpgRouting2系の設定例
pgAdmin4の左側ツリー上で各ポート番号の alphashape データベースを選択し、上部メニューから Tools / Query Tool メニューを選択します。
初期状態ではpgRouting拡張機能が有効になっていないので、中央の Query Editor パネルに CREATE EXTENSION pgRouting CASCADE; と入力後、上部ツールバーの実行ボタン▶をクリックし、pgRoutingと、依存関係にあるPostGISを有効にします。
念のため、それぞれのpgRouting/PostGISバージョンを確認しておきます。
Query Editor を SELECT pgr_version() AS pgrouting, postgis_version() AS postgis; に置き換えて、実行ボタン▶をクリックすると、それぞれに異なるバージョンがインストールされていることが分かります。
- pgRouting2系(localhost:54322):
- pgRouting3系(localhost:54323):

入力データの読み込み

何か良いWKT形式の入力データはないかと探していたところ、PostGISのドキュメントの ST_ConcaveHull の説明欄のリンクに、S字型の良いサンプルがありましたので、まずはこちらを読み込んでみます。

現実世界の例と技術面でのしっかりした説明は、http://www.bostongis.com/postgis_concavehull.snippet にあります。

Oracle 11G R2で導入された凹包のデモンストレーションに関するSimon Greenerさんの記事も見てください。http://www.spatialdbadvisor.com/oracle_spatial_tips_tricks/172/concave-hull-geometries-in-oracle-11gr2 にあります。凸包に対する目標割合を0.75にした際の形状がSimonさんOracleのSDO_CONCAVEHULL_BOUNDARYで得た形状と似ています。

上記のBostonGISのサイト(http://www.bostongis.com/postgis_concavehull.snippet) にある以下のSQL文をコピーして、 Query Editor で実行します。
```
CREATE TABLE s_test(geom geometry); 
	INSERT INTO s_test(geom) VALUES(ST_GeomFromText('MULTIPOINT(
```

(120 -224), (185 -219), (190 -234), (200 -246), (212 -256), (227 -261),
(242 -264), (257 -265), (272 -264), (287 -263),(302 -258), (315 -250),
(323 -237), (321 -222), (308 -213), (293 -208), (278 -205), (263 -202),
(248 -199), (233 -196), (218 -193), (203 -190), (188 -185), (173 -180),
(160 -172), (148 -162), (138 -150), (133 -135), (132 -120), (136 -105),
(146 -92), (158 -82), (171 -74), (186 -69), (201 -65), (216 -62),
(231 -60), (246 -60), (261 -60), (276 -60), (291 -61), (306 -64),
(321 -67), (336 -72), (349 -80), (362 -89), (372 -101), (379 -115),
(382 -130), (314 -134), (309 -119), (298 -108), (284 -102), (269 -100),
(254 -99), (239 -100), (224 -102), (209 -107), (197 -117), (200 -132),
(213 -140), (228 -145), (243 -148), (258 -151), (273 -154), (288 -156),
(303 -159), (318 -163), (333 -167), (347 -173), (361 -179), (373 -189),
(383 -201), (389 -215), (391 -230), (390 -245), (384 -259), (374 -271),
(363 -282), (349 -289), (335 -295), (320 -299), (305 -302), (290 -304),
(275 -305), (259 -305), (243 -305), (228 -304), (213 -302), (198 -300),
(183 -295), (169 -289), (155 -282), (143 -272), (133 -260), (126 -246),
(136 -223), (152 -222), (168 -221), (365 -131), (348 -132), (331 -133),
(251 -177), (183 -157), (342 -98), (247 -75), (274 -174), (360 -223),
(192 -85), (330 -273), (210 -283), (326 -97), (177 -103), (315 -188),
(175 -139), (366 -250), (321 -204), (344 -232), (331 -113), (162 -129),
(272 -77), (292 -192), (144 -244), (196 -272), (212 -89), (166 -236),
(238 -167), (289 -282), (333 -187), (341 -249), (164 -113), (238 -283),
(344 -265), (176 -248), (312 -273), (299 -85), (154 -261), (265 -287),
(359 -111), (160 -150), (212 -169), (351 -199), (160 -98), (228 -77),
(376 -224), (148 -120), (174 -272), (194 -100), (292 -173), (341 -212),
(369 -209), (189 -258), (198 -159), (275 -190), (322 -82))') ) ;
```

SELECT geom FROM s_test; を実行し、結果のヘッダにあるViewアイコンをクリックします。
Geometry Viewer という名前でタブが追加され、ジオメトリの情報を確認することができます。

AlphaShapeの取得

入力データの準備も整ったので、AlphaShapeの取得を行っていきます。
pgRouting2系と3系で、シグネチャーが大きく変わるので、それぞれごとに見ていきます。

pgRouting2系でのAlphaShapeの取得

関数ドキュメント:

pgr_alphaShape: https://docs.pgrouting.org/2.6/en/pgr_alphaShape.html
pgr_pointsAsPolygon : https://docs.pgrouting.org/2.6/en/pgr_pointsAsPolygon.html#pgr-points-as-polygon

pgRouting2系では、pgr_alphaShape 関数で、 id, x, y を返すSQLから、ポリゴンの点列(x, y)を返します。(ポリゴンが複数に分かれたり、穴ポリゴンがある場合は、セパレータとして、NULL, NULL が利用されます。)
pgr_alphaShape 関数を使っても、可視化するまでが少し大変なので、 pgr_alphaShape 関数をラップした pgr_pointsAsPolygon 関数を使って、S字型入力データのAlphaShapeを表示してみます。

入力データは MultiPoint ジオメトリ形式で投入されているので、まずは、PostGISの ST_Dump 関数(http://cse.naro.affrc.go.jp/yellow/pgisman/1.5.1/ST_Dump.html) を利用して、id, x, y を列に取るクエリーを作成します。

SELECT
  (a.dump).path[1] AS id,
  ST_X((a.dump).geom) AS x,
  ST_Y((a.dump).geom) AS y
FROM (
  SELECT ST_Dump(geom) AS dump
  FROM s_test
) AS a;

無事に MultiPoint ジオメトリから id, x, y 列への変換が確認できたので、このクエリを第1引数として、 pgr_pointsAsPolygon 関数に text 型で渡します。

SELECT pgr_pointsAsPolygon('
  SELECT
    (a.dump).path[1] AS id,
    ST_X((a.dump).geom) AS x,
    ST_Y((a.dump).geom) AS y
  FROM (
    SELECT ST_Dump(geom) AS dump
    FROM s_test
  ) AS a') AS geom;

少し歯抜けの形でS字型のAlphaShapeジオメトリが返ってきましたが、分かりにくいので元の点列と合わせて表示してみます。

SELECT geom FROM s_test
UNION ALL
SELECT pgr_pointsAsPolygon('
  SELECT
    (a.dump).path[1] AS id,
    ST_X((a.dump).geom) AS x,
    ST_Y((a.dump).geom) AS y
  FROM (
    SELECT ST_Dump(geom) AS dump
    FROM s_test
  ) AS a') AS geom;

大分見やすくなりました。
pgr_alphaShape 関数と同様、pgr_pointsAsPolygon 関数でも、第2引数に alpha 値を取れるので、少し値を調整してみます。

SELECT geom FROM s_test
UNION ALL
SELECT pgr_pointsAsPolygon('
  SELECT
    (a.dump).path[1] AS id,
    ST_X((a.dump).geom) AS x,
    ST_Y((a.dump).geom) AS y
  FROM (
    SELECT ST_Dump(geom) AS dump
    FROM s_test
  ) AS a', 300) AS geom;

`alpha=100`	`alpha=200`	`alpha=300`	`alpha=1200`

CGAL の2D AlphaShapeのドキュメント(https://doc.cgal.org/latest/Alpha_shapes_2/group__PkgAlphaShapes2Ref.html) によると、alpha 値の平方根の半径のスプーンでくり抜くとありますので、スプーンのサイズが大きくなるほど、スプーンが引っかかる点が多くなって、その分、中身が詰まっている形になっていそうです。

pgRouting3系でのAlphaShapeの取得

関数ドキュメント:

pgr_alphaShape: https://docs.pgrouting.org/3.1/en/pgr_alphaShape.html

pgRouting3系では、CGAL への依存回避に伴い、 pgr_alphaShape 関数のインターフェースも刷新され、第1引数がクエリーテキストからジオメトリに、また第2引数は alpha 値から、平方根の spoon_radius 値に変更されています。

まずは第2引数を指定せずにクエリを実行してみます。

SELECT geom FROM s_test
UNION ALL
SELECT pgr_alphaShape(geom) AS geom FROM s_test;

次に、第2引数(spoon_radius)を指定・調整してクエリを実行してみます。

SELECT geom FROM s_test
UNION ALL
SELECT pgr_alphaShape(geom, 20) AS geom FROM s_test;

`spoon_radius=10`	`spoon_radius=15`	`spoon_radius=25`	`spoon_radius=35`

細かい点で違いはありますが、概ね、pgRouting2系と同様の結果が返っていそうです。

PostGISのST_ConcaveHull関数との比較

最後に、入力データで参考にしたBostonGISのサイト(http://www.bostongis.com/postgis_concavehull.snippet) でのPostGISの ST_ConcaveHull 関数の結果とも比較してみます。

関数ドキュメント:

ST_ConcaveHull: http://cse.naro.affrc.go.jp/yellow/pgisman/2.2.0/ST_ConcaveHull.html

SELECT geom FROM s_test
UNION ALL
SELECT ST_ConcaveHull(ST_Collect(geom), 0.999) AS geom FROM s_test;

`target_percent=0.999, allow_holes=false`	`target_percent=0.90, allow_holes=false`	`target_percent=0.90, allow_holes=false`

概ね、BostonGISのサイトに記載通りの結果となってますが、 target_percent 値が小さい場合に、点列からのバッファが取られているようで、経路探索結果の到達圏の表示などの用途では利用しづらいかもしれません。

後片付け

pgAdmin4は、PostgreSQLサーバとは別のサービスとして起動されているので、ブラウザのタブを閉じた後、OSのツールバー上にあるゾウのアイコン上で右クリックして、コンテキストメニューから Shut down server を選択して、サービスを終了します。

今回構築したpgRoutingのDocker環境は、再度試す場合は、コンソールから以下を実行します。
(次回起動時は docker-compose start を実行します。)

$ docker-compose stop

データベースボリュームやDockerコンテナ・イメージを完全に削除する場合は、コンソールから以下を実行します。

$ docker-compose down --rmi all -v

おわりに

とりとめもなく、書き綴ってしまいましたが、pgRoutingの環境も、数年前に比べると、Dockerで気軽に環境構築できるようになったり、pgAdmin4でジオメトリの可視化もできるようになって、格段に使いやすくなったと感じる今日この頃です。

pgRouting3系のAlphaShapeアルゴリズム内部や、Isochrone との関連については、別の記事で投稿できたらと思ってますので、またよろしくお願いします。

https://github.com/pgRouting/pgrouting/issues/235 ↩