PostGISでポイントのメッシュ集計を効率よく実行する

## TL;DR

- `ST_SnapToGrid` により、素朴にポリゴンで集計するより500倍くらいコストを小さく出来た（所要時間ではないことに留意）

```sql
-- 3次メッシュで集計する例
CREATE TABLE grid_agg AS (
  WITH points AS (
	  -- 最近傍のグリッドにスナップされるので、各点がメッシュの中心となるようなグリッドを定義する
    SELECT ST_SnapToGrid(
      geom,
      (1 * 0.0125) * 0.5, (0.666666666 * 0.0125) * 0.5, -- origin = origin_pos - size / 2
      1 * 0.0125, 0.666666666 * 0.0125 -- size
    ) as geom
    FROM poi
  )
  SELECT 
	  -- メッシュの中心の座標からメッシュを図形を作成する
    ST_MakeEnvelope(
      ST_X(geom) - (1 * 0.0125) * 0.5,
      ST_Y(geom) - (0.666666666 * 0.0125) * 0.5,
      ST_X(geom) + (1 * 0.0125) * 0.5,
      ST_Y(geom) + (0.666666666 * 0.0125) * 0.5
    ) as geom,
    COUNT(geom) as count
  FROM points
  GROUP BY geom
);
```

## メッシュ集計について

- 素朴には以下の手順が考えられる：
    1. メッシュのポリゴンを作成する
    2. ポイントデータと空間結合する
    - たぶんこの方法は遅い（ポリゴンを作成するのがめんどう・空間結合がかなりコストかかるはず）
- ポイントデータの集計にそういう行儀の良い方法を取る必要はないだろうから、`ST_SnapToGrid`を利用してポイントデータをグルーピングする手法をとってみる
    1. `ST_SnapToGrid`で、定義したグリッドに地物を「寄せる」
    2. 寄せた後の地物で`GROUP BY` し、COUNTで数を集計する

## 検証

- 一定の範囲内にポイントを作成する（ここでは100万ポイント。なお1000万ポイントでもリーズナブルな速度で動作した。1億だとそもそもポイントを作るのが大変だった。）
    
    ![Screenshot_2024-08-04_at_11.14.27.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/73197/f8ab1640-5ff1-ec71-787c-5f61010611f0.png)*©︎OpenStreetMap Contributors.*
    
- 冒頭のSQLを実行する（1秒未満で完了する、cost=1106505.39）
    
    ```
                                                                                            QUERY PLAN                                                                              
               
    --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
    -----------
     Finalize GroupAggregate  (cost=114692.88..1106505.39 rows=1000000 width=72)
       Group Key: (st_snaptogrid(poi.geom, '0.00625'::double precision, '0.0041666666625'::double precision, '0.0125'::double precision, '0.008333333325'::double precision))
       ->  Gather Merge  (cost=114692.88..322338.72 rows=833334 width=72)
             Workers Planned: 2
             ->  Partial GroupAggregate  (cost=113692.85..225151.28 rows=416667 width=72)
                   Group Key: (st_snaptogrid(poi.geom, '0.00625'::double precision, '0.0041666666625'::double precision, '0.0125'::double precision, '0.008333333325'::double precis
    ion))
                   ->  Sort  (cost=113692.85..114734.52 rows=416667 width=64)
                         Sort Key: (st_snaptogrid(poi.geom, '0.00625'::double precision, '0.0041666666625'::double precision, '0.0125'::double precision, '0.008333333325'::double p
    recision))
                         ->  Parallel Seq Scan on poi  (cost=0.00..63610.04 rows=416667 width=64)
     JIT:
       Functions: 10
       Options: Inlining true, Optimization true, Expressions true, Deforming true
    (12 rows)
    ```
    
    ![Screenshot_2024-08-04_at_11.18.56.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/73197/962f4809-032f-9742-9242-80e8f777c7d9.png)*©︎OpenStreetMap Contributors.*

    
    ![Screenshot_2024-08-04_at_11.19.30.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/73197/cd25e2aa-b7cc-f0e7-6f29-f77d3d857242.png)*ポイントの疎密と見比べると集計できていそう*

## 参考：ちゃんとポリゴンで集計する場合のコストを調べておく

1. 3次メッシュのグリッドのテーブルを作成しておく（QGISの[地域メッシュプラグイン](https://github.com/MIERUNE/qgis-japan-mesh)からPostGISに直接書き込めて便利だった）
2. グリッドのテーブルに空間インデックスを作成しておく（以下では`thir` というテーブル、thirdをタイポした）
3. SQLを実行する（cost=639598488.91、前述の手法の578倍）
    
    ```sql
    CREATE TABLE grid_agg_poly AS (
      SELECT thir.geom AS geom, COUNT(poi.geom) AS count
      FROM thir
      LEFT JOIN poi
      ON ST_Contains(thir.geom, poi.geom)
      GROUP BY thir.geom
    );
    
                                                      QUERY PLAN                                                  
    --------------------------------------------------------------------------------------------------------------
     Finalize GroupAggregate  (cost=639313115.89..639598488.91 rows=1126400 width=128)
       Group Key: thir.geom
       ->  Gather Merge  (cost=639313115.89..639575960.91 rows=2252800 width=128)
             Workers Planned: 2
             ->  Sort  (cost=639312115.87..639314931.87 rows=1126400 width=128)
                   Sort Key: thir.geom
                   ->  Partial HashAggregate  (cost=624894942.36..639146644.15 rows=1126400 width=128)
                         Group Key: thir.geom
                         Planned Partitions: 16
                         ->  Nested Loop Left Join  (cost=0.29..588128721.15 rows=331413825 width=152)
                               ->  Parallel Seq Scan on thir  (cost=0.00..27681.33 rows=469333 width=120)
                               ->  Index Scan using poi2_geom_idx on poi2  (cost=0.29..1252.06 rows=100 width=32)
                                     Index Cond: (geom @ thir.geom)
                                     Filter: st_contains(thir.geom, geom)
     JIT:
       Functions: 12
       Options: Inlining true, Optimization true, Expressions true, Deforming true
    (17 rows)
    ```

## 参考リンク

https://gis.stackexchange.com/questions/274966/grouping-points-per-grid-and-summing-their-value-on-postgis

https://postgis.net/docs/ST_SnapToGrid.html

## 宣伝

「位置エン本」に続く第2作、**現場のプロがわかりやすく教える位置情報デベロッパー養成講座: デジタルマップ配信のためのサーバーサイド実装**が2024年8月末に発売されます。PostGISの説明に注力しているほか、ラスター配信のテーマとして衛星データを取り扱っています。興味のある方はぜひ書店でご覧になってください。

https://amzn.asia/d/1t3sQsj