pgRoutingのAlphaShapeについての雑記

この記事は、[RDBMS-GIS(MySQL,PostgreSQLなど) Advent Calendar 2020](https://qiita.com/advent-calendar/2020/rdbms_gis) の13日目の内容になります。

# はじめに

[pgRouting](https://live.osgeo.org/ja/overview/pgrouting_overview.html) は、オープンソースのPostgreSQLの拡張機能で、PostGISが空間演算に特化したものであるのに対し、経路探索などのアルゴリズムに特化したものになります。

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

![pgrouting3_alphashape.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/9c5560d0-9c27-9a2c-9be9-20313148a864.png)

この [AlphaShape](https://docs.pgrouting.org/latest/en/pgr_alphaShape.html) の機能について、私もpgRouting2系で改善に関わり([^1])、また、pgRouting3系から、内部で依存していた [CGAL](https://www.cgal.org/) から脱却した独自実装となったこともあり、検証も兼ねて、pgRouting2系と3系の違い、また、PostGISの `ST_ConcaveHull` との相違についても確認してみます。


# 環境構築

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

1. Docker Desktop のインストール  
   https://www.docker.com/products/docker-desktop から各OS用のインストーラーをダウンロードしてインストールします。
2. pgAdmin4 のインストール  
   https://www.pgadmin.org/download/ から各OS用のインストーラーをダウンロードしてインストールします。
3. 作業用のフォルダ(`alphashape`)を作成し、`docker-compose.yml` ファイルを以下の内容で作成して保存します。  
    ※ポート番号の重複を避けるため、pgRouting2系をポート番号 `54322`、pgRouting3系をポート番号 `54323` で指定しています。
    ※ボリュームについても、混同を避けるため、別の名前(`db-data2`・`db-data3`)を持たせています。
    
    ```yaml
    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:
    ```

4. コンソールから作業用のフォルダに移動し、以下でDockerサービスを起動します。
    
    ```bash
    $ 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
    ```

5. pgAdmin4の左側ツリーの `Servers` 上で右クリックして、 `Create` / `Server...` メニューを選択し、  
    ![pgAdmin_01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/38e6105d-f255-6fc2-3e79-3a727cabdbec.png)

    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系の設定例  
    ![pgAdmin_02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/82a499d1-039f-8aa6-0d02-56310805f639.png) ![pgAdmin_03.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/eebcd1aa-2f92-e4ee-07ce-f6e9fddb0d74.png)
6. pgAdmin4の左側ツリー上で各ポート番号の `alphashape` データベースを選択し、上部メニューから `Tools` / `Query Tool` メニューを選択します。  
    ![pgAdmin_04.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/3d07c29d-77cc-92d3-1bad-d0c9629ad159.png)
7. 初期状態ではpgRouting拡張機能が有効になっていないので、中央の `Query Editor` パネルに `CREATE EXTENSION pgRouting CASCADE;` と入力後、上部ツールバーの実行ボタン`▶`をクリックし、pgRoutingと、依存関係にあるPostGISを有効にします。
    ![pgAdmin_05.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/ce022841-8d60-5915-a4ec-8b928ba20581.png)
8. 念のため、それぞれのpgRouting/PostGISバージョンを確認しておきます。  
   `Query Editor` を `SELECT pgr_version() AS pgrouting, postgis_version() AS postgis;` に置き換えて、実行ボタン`▶`をクリックすると、それぞれに異なるバージョンがインストールされていることが分かります。
    * pgRouting2系(`localhost:54322`):  
      ![pgrouting2_version.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/caa3187a-0acb-1024-604e-25bdb019f8c6.png)
    * pgRouting3系(`localhost:54323`):  
      ![pgrouting3_version.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/e236d684-c3a9-de3c-c2ad-fba01e97f1ee.png)

# 入力データの読み込み

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

http://cse.naro.affrc.go.jp/yellow/pgisman/2.2.0/ST_ConcaveHull.html  
> 現実世界の例と技術面でのしっかりした説明は、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で得た形状と似ています。

1. 上記のBostonGISのサイト(http://www.bostongis.com/postgis_concavehull.snippet) にある以下のSQL文をコピーして、 `Query Editor` で実行します。
    
    ```sql
    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))') ) ;
    ```

2. `SELECT geom FROM s_test;` を実行し、結果のヘッダにあるViewアイコンをクリックします。  
    ![s_test1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/3b32d7eb-c3de-763d-eacd-3ad2bb912009.png)
3. `Geometry Viewer` という名前でタブが追加され、ジオメトリの情報を確認することができます。  
    ![s_test2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/224ba30e-7505-2939-64e0-6c396b56e2b3.png)

# 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` を列に取るクエリーを作成します。

```sql
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;
```

![vertices1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/e37c9c7e-b2aa-e3ee-f04d-c89bff4c08a1.png)

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

```sql
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;
```

![pgrouting2_alphashape1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/6e33c4e7-bc6a-92c2-c2c2-e29af112de20.png)

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

```sql
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;
```

![pgrouting2_alphashape2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/b51fac1f-7318-38e9-ecec-043647b26fb3.png)

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

```sql
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` |
|-------------|-------------|-------------|--------------|
| ![pgrouting2_alpha100.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/51d7698c-f4d4-d447-d43c-6cdb289e9de1.png) | ![pgrouting2_alpha200.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/2d64a306-5a88-9a45-f94a-0d478dd743f3.png) | ![pgrouting2_alpha300.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/46383dd0-1ebc-e820-ddd6-bc5cfa26a342.png) | ![pgrouting2_alpha1200.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/f573fc03-d42c-0179-db83-34a941c712f1.png) |

[CGAL](https://www.cgal.org/) の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](https://www.cgal.org/) への依存回避に伴い、 `pgr_alphaShape` 関数のインターフェースも刷新され、 第1引数がクエリーテキストからジオメトリに、また第2引数は `alpha` 値から、平方根の `spoon_radius` 値に変更されています。

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

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

![pgrouting3_alphashape.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/f40314f2-e144-5eb7-c8ed-55d6533f5ee4.png)

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

```sql
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_spoon10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/1aaf9e68-8bdb-0c0b-9b5f-bc6d91063d0f.png) | ![pgrouting3_alpha15.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/213b213b-a078-ac4f-c229-ec5c089b21d8.png) | ![pgrouting3_spoon20.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/bf122e64-11f8-00ca-520d-a91ba58542e9.png) | ![pgrouting3_spoon35.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/97036b67-bbc4-eae9-e857-9e5d798a877e.png) |

細かい点で違いはありますが、概ね、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

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

![postgis_concavehull.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/91380955-0e7f-f5f3-c186-7ce56f36b52c.png)

| `target_percent=0.999, allow_holes=false` | `target_percent=0.90, allow_holes=false` | `target_percent=0.90, allow_holes=false` |
|-------------------------------------------|------------------------------------------|-------------------------------------------|
| ![postgis_999_resize.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/b248188a-b6ef-4f84-9c34-c8e6906ee224.png) | ![postgis_90f.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/46c39ebd-a02f-5513-11a6-ece42ee236b4.png) | ![postgis_90t.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/5a9ba2f1-bb56-eae0-7056-3b3e8302cada.png) |

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


# 後片付け

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

![macos_pgadmin4_shutdown.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/2efbe341-c1cb-1f32-b542-9ea6f7c0687b.png) ![windows_pgadmin4_shutdown.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104365/2e62c31a-b591-20b3-e408-a40546d7e24d.png)

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

```bash
$ docker-compose stop
```

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

```bash
$ docker-compose down --rmi all -v
```

# おわりに

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

pgRouting3系のAlphaShapeアルゴリズム内部や、[Isochrone](https://wiki.openstreetmap.org/wiki/Isochrone) との関連については、別の記事で投稿できたらと思ってますので、またよろしくお願いします。

[^1]: https://github.com/pgRouting/pgrouting/issues/235