目的

ポリゴンの重なり具合でヒートマップを作りたい。

はいさい、@tohka382さんの記事でgeopandasを使う方法が紹介されています。
https://qiita.com/tohka383/items/23c7074f68293c5c4609

私も同じ問題をやっています、ていうか@tohka328さんは私の課題の手伝いをしてくれたはりました。私の方法ではPostGIS上でGIS操作をしていました。データベースを使わない実装vs使った実装をくらべてみます。

方法

@tohka382さんの実装は4個あります。

方針1 何も考えない
方針2 分割する
方針3 空間インデックスを利用し、対象を絞る
方針4 空間インデックスで交差判定

方針3がわりと定石っぽい方針と思います。方針2と3はほぼ同じ計算速度。方針3から4がなんか私は見たことないやり方で、100倍速くなる！ということです。なので、私は方針1,3,4をPostGISでなぞってみます。サンプルデータも同じ。

ちなみにPostGISはたまにしか使わないときは、いつものマシンにDockerで入れるのが便利です。https://hub.docker.com/r/kartoza/postgis/

方針0　geopandasで

@tohka382さんのコードのコピペ、私の環境で使ってみました。私のほうがちょっと遅いですが大体似たような速度でした。

方針1、1行（1,000点）の処理で55秒くらい（c.f. @tohka382「1分弱」）

[01:54:23.654162] starting...
[01:54:51.979846] finish to load data
[01:54:55.770756] finish to generate buffer data
[01:55:50.718427] finish to calculate row #0
[01:56:46.045506] finish to calculate row #1
[01:57:38.675440] finish to calculate row #2
...

方針3、1行（1,000点）の処理で18秒くらい（c.f. @tohka382「14秒くらい」）

[02:00:46.610632] starting...
[02:01:13.854278] finish to load data
[02:01:17.824631] finish to generate buffer data
[02:01:35.362627] finish to calculate row #0
[02:01:53.103917] finish to calculate row #1
[02:02:11.337160] finish to calculate row #2

方針4、速いです！　バッファ作り終わってから最後まで、全域（1,000,000点）で60秒。（c.f. @tohka382「対象領域全体の交差判定が42秒ほどで完了」）

[02:18:00.297249] starting...
[02:18:27.284618] finish to load data
[02:18:31.108596] finish to generate buffer data
[02:18:31.271433] finish to calculate row #0
[02:18:31.325430] finish to calculate row #1
[02:18:31.370654] finish to calculate row #2
...
[02:19:30.879607] finish to calculate row #997
[02:19:30.945383] finish to calculate row #998
[02:19:31.012303] finish to calculate row #999

方針1 空間インデックスなし

さて、これをPostGISでやってみます。

処理の流れは

病院の点データを読み込む
バッファでポリゴン（円）にする
処理する領域を指定する
実際の計算
結果をラスターに書き出す

4番目のところがいろいろ工夫するところであとは全く変わりません。

この流れをPostGISで大まかに再現するために@tohka382さんの方針1を作ってみました。

#!/usr/bin/python

from datetime import datetime
import numpy as np
import psycopg2
import subprocess
import shlex
from osgeo import gdal, osr

print(f"[{datetime.now().strftime('%H:%M:%S.%f')}] starting...")

# データベースの用意
conn = psycopg2.connect(database='finn', user='finn', password='finn')
conn.autocommit = True
cur = conn.cursor()

schema = 'hospital'
cur.execute(f'CREATE SCHEMA IF NOT EXISTS {schema}')
cur.execute(f'SET search_path TO {schema}, public;')

# サンプルデータの読み込み
# 国土数値情報の医療機関データ（点）
# shp2pgsqlを使いました
cur.execute(f'DROP TABLE IF EXISTS {schema}.hospital;')
cmd = shlex.split(f'shp2pgsql  -s 4326 -W cp932 -g geom_pnt -I ./P04-20_GML/P04-20.shp {schema}.hospital')
p1 =  subprocess.Popen(cmd, stdout=subprocess.PIPE)
p2 = subprocess.Popen(shlex.split('psql -q -d finn -U finn '), stdin=p1.stdout)
p2.communicate()

print(f"[{datetime.now().strftime('%H:%M:%S.%f')}] finish to load data")

# 点データをもとに 0.05 deg ~ 5 km バッファを生成
# shpファイルから取り込んだテーブルにポリゴン用の行をつけてST_Buffer()で
cur.execute(f'''
ALTER TABLE hospital ADD COLUMN geom_poly geometry(POLYGON, 4326);

INSERT INTO hospital (geom_poly)
SELECT ST_Buffer(geom_pnt, 0.05)
FROM hospital;
''')

print(f"[{datetime.now().strftime('%H:%M:%S.%f')}] finish to generate buffer data")

# (137, 35) - (138, 36) の範囲に対して、解像度 0.001 deg ~ 100 m の
# 間隔で、 regions のポリゴンがいくつ含まれるか計算を行う
# 結果は、取り扱いがしやすいようにラスタ形式で保存する
bounds = (137, 35, 138, 36)
resol  = 0.001
width  = round((bounds[2] - bounds[0]) / resol)
height = round((bounds[3] - bounds[1]) / resol)

######################################
# ここまでは準備段階。ここから計算です。
######################################

# 格子点を先に作ってしまう
# generate_series()を使うのが速くて便利な模様

cur.execute(f'''
DROP TABLE IF EXISTS grid;
CREATE TABLE IF NOT EXISTS grid (
  gid SERIAL PRIMARY KEY,
  geom geometry(POINT, 4326),
  idx INTEGER,
  jdx INTEGER
);

INSERT INTO grid (idx, jdx, geom)
SELECT i idx, j jdx, 
  ST_SetSRID(
    ST_MakePoint({bounds[0]} + i * {resol}, {bounds[1]} + j * {resol})
    , 4326) geom
FROM generate_series(0, {int(width)}-1) i, generate_series(0, {int(height)}-1) j
''')

print(f"[{datetime.now().strftime('%H:%M:%S.%f')}] finish to generate grid points")

# 個数管理用のテーブル
cur.execute(f'''
DROP TABLE IF EXISTS points;
CREATE TABLE IF NOT EXISTS points
(
  idx INTEGER,
  jdx INTEGER,
  cnt INTEGER
);
''')

# 一行（1000点）ごとに処理。
# 格子点のうちポリゴンと交わるものを見つけて（テーブル＝intrsct）勘定する（テーブル＝points）　
for j in range(int(height)):
    cur.execute(f'''
    WITH work AS (
      SELECT * FROM grid g
      WHERE g.jdx = {j}
    ), intrsct AS (
      SELECT h.gid hospital_id, w.gid grid_id, w.idx idx, w.jdx jdx
      FROM hospital h
      JOIN work w
      ON ST_Intersects(h.geom_poly, w.geom)
    )
    INSERT INTO points(idx, jdx, cnt)
    SELECT idx, jdx, count(*)
    FROM intrsct
    GROUP BY idx, jdx;
    ''')

    print(f"[{datetime.now().strftime('%H:%M:%S.%f')}] finish to find intersection row # {j}")

######################################
# 計算終わり。出力です
######################################

# 結果をnumpy の array
raster_data = np.zeros((width, height), dtype=np.uint32)

# 結果をラスタファイルに出力する
ds_filename = r'./regions_count_%03d_%03d.tif' % bounds[0:2]

# ファイルを用意
gtiff_driver = gdal.GetDriverByName('GTiff')
raster_ds = gtiff_driver.Create(ds_filename, width, height, 1, gdal.GDT_UInt32)

# ラスタの位置情報および座標系情報の設定
raster_ds.SetGeoTransform([bounds[0], resol, 0, bounds[3], 0, -resol])
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
raster_ds.SetProjection(srs.ExportToWkt())

# データの書き込み
raster_ds.GetRasterBand(1).WriteArray(raster_data)
# とじる
raster_ds = None

結果。一行に40秒くらいかかってます。geopandasでやったの方針1（c.f.55秒くらい）と似たような感じです。

[02:24:59.704048] starting...
Shapefile type: Point
Postgis type: POINT[2]
                    addgeometrycolumn
---------------------------------------------------------
 hospital.hospital.geom_pnt SRID:4326 TYPE:POINT DIMS:2
(1 row)

[02:25:34.542139] finish to load data
[02:25:40.740715] finish to generate buffer data
[02:25:45.427586] finish to generate grid points
[02:26:30.001241] finish to find intersection row # 0
[02:27:12.110846] finish to find intersection row # 1
[02:27:59.852485] finish to find intersection row # 2
...

方針3 空間インデックスあり

次は空間インデックスを使います。

ST_Bufferでポリゴンを作った後に

# ポリゴンに空間インデックス
cur.execute(f'''
DROP INDEX IF EXISTS hospital_geom_poly_index;
CREATE INDEX hospital_geom_poly_index ON hospital USING gist(geom_poly);
''')

print(f"[{datetime.now().strftime('%H:%M:%S.%f')}] finished generating spatial index on polygons")

格子点を作った後にも格子点にインデックス。あとに書くようにこれはむしろやらんほうがいいみたいでした。

# 格子点にも空間インデックス
# （これはむしろやらんほうがいいかも？）
cur.execute(f'''
DROP INDEX IF EXISTS grod_geom_index;
CREATE INDEX grid_geom_index ON grid USING gist(geom);
''')

print(f"[{datetime.now().strftime('%H:%M:%S.%f')}] finish to generate spatial index on grid")

そして交わりを求める。Pythonの行ごとのループは外しちゃいました。

# 格子点のうちポリゴンと交わるものを見つけて（テーブル＝intrsct）勘定する（テーブル＝points）　
cur.execute(f'''
DROP TABLE IF EXISTS points;
CREATE TABLE IF NOT EXISTS points
(
idx INTEGER,
jdx INTEGER,
cnt INTEGER
);

WITH intrsct AS (
SELECT h.gid hospital_id, g.gid grid_id, g.idx idx, g.jdx jdx
FROM hospital h
JOIN grid g
ON ST_Intersects(h.geom_poly, g.geom)
) 
INSERT INTO points(idx, jdx, cnt)
SELECT idx, jdx, count(*)
FROM intrsct
GROUP BY idx, jdx;
''')

print(f"[{datetime.now().strftime('%H:%M:%S.%f')}] finish to find intersection")

で、結果はこれです。バッファを作り終わってから数えて99秒です。かなり速い！

geopandasでやった方針3では1行（1000点）に18秒かかってたので全然速いです。むしろ方針4の結果に近いですね。単純に空間インデックスをかけただけなんですが。

[03:30:07.895797] starting...
Shapefile type: Point
Postgis type: POINT[2]
                    addgeometrycolumn
---------------------------------------------------------
 hospital.hospital.geom_pnt SRID:4326 TYPE:POINT DIMS:2
(1 row)

[03:30:52.005194] finish to load data
[03:30:58.761247] finish to generate buffer data
[03:31:00.641589] finished generating spatial index on polygons
[03:31:05.263974] finish to generate grid points
[03:31:18.533319] finish to generate spatial index on grid
[03:32:37.507074] finish to find intersection

あと気づいたんが格子点に空間インデックスを作るのにえらい時間かかってます。点が多いからでしょうか？試しに点に空間インデックスなしでやると

# 格子点にも空間インデックス
# （これはむしろやらんほうがいいかも？）
cur.execute(f'''
DROP INDEX IF EXISTS grod_geom_index;
-- CREATE INDEX grid_geom_index ON grid USING gist(geom);
''')
print(f"[{datetime.now().strftime('%H:%M:%S.%f')}] not generating spatial index on grid")

こうなりました。交差をもとめるのに72秒かかってますが格子に空間インデックスを付けたときは79秒。同じかむしろ早いのか？というわけで、このデータと処理のときは点には空間インデックスは要らないみたいです。点でも空間インデックスがあったら速くなる時もあるそうなんですがどういう時なのか

[03:48:21.297336] starting...
Shapefile type: Point
Postgis type: POINT[2]
                    addgeometrycolumn
---------------------------------------------------------
 hospital.hospital.geom_pnt SRID:4326 TYPE:POINT DIMS:2
(1 row)

[03:48:56.233020] finish to load data
[03:49:02.165293] finish to generate buffer data
[03:49:04.068456] finished generating spatial index on polygons
[03:49:07.798058] finish to generate grid points
[03:49:07.798419] not generating spatial index on grid
[03:50:19.134978] finish to find intersection

方針4 クリップしたポリゴンに空間インデックスを作って、空間インデックスの交わりのみで判定

さて、ここまでのおさらいは、

ポリゴンには空間インデックスを！
点には空間インデックスが今回は意味なかった
「普通に」空間インデックスをかけたらgeopandasよりPostGISのほうがめちゃ速い
PostGISのGISTインデックスは多分geopandasで使ってるRTreeインデックスより性能がいい？

で、このPostGISのGiSTインデックスを使った実装でも、私の応用には少し遅いんですよね…もう一桁くらい速くなってくれるとありがたいんですが。というわけで、@tohka382さんが使った方法をPostGISで試してみようと思います。

多分こんな感じで…

元のポリゴンには空間インデックスをつくる
Pythonで行ごとに処理
ST_Intersection(h.geom_poly, ST_SetSRID(ST_GeomFromText('POLYGON((...))'), 4326))みたいな感じでgeopandas.GeoDataFrame.clip()がやってるような短冊を作る
JOINのときはST_Intersect()のかわりに&&オペレータをつかう（空間インデックスの交差）

書いたのがこれ

# 個数管理用のテーブル
cur.execute(f'''
DROP TABLE IF EXISTS points;
CREATE TABLE IF NOT EXISTS points
(
  idx INTEGER,
  jdx INTEGER,
  cnt INTEGER
);
''')

cur.execute(f'''
DROP TABLE IF EXISTS for_idx;
CREATE TABLE IF NOT EXISTS for_idx
(
  geom geometry(POLYGON, 4326),
  hospital_id INTEGER
);
CREATE INDEX work_geom_idx ON for_idx USING gist(geom);
''')


# 一行（1000点）ごとに処理。
# 格子点のうちポリゴンと交わるものを見つけて（テーブル＝intrsct）勘定する（テーブル＝points）　
for j in range(int(height)):
    wkt = f'''POLYGON((
    {bounds[0]} {bounds[3]-(j+.6)*resol}, 
    {bounds[0]} {bounds[3]-(j+.4)*resol}, 
    {bounds[2]} {bounds[3]-(j+.4)*resol}, 
    {bounds[2]} {bounds[3]-(j+.6)*resol}, 
    {bounds[0]} {bounds[3]-(j+.6)*resol} 
    ))'''
    
    cur.execute(f'''
    TRUNCATE TABLE for_idx;
    
    WITH box AS ( 
      SELECT ST_SetSRID(ST_GeomFromText('{wkt}'), 4326) box
    )
    INSERT INTO for_idx(geom, hospital_id)
    SELECT ST_INTERSECTION(h.geom_poly, box.box) geom, h.gid hospital_id
    FROM hospital h
    JOIN box
    ON h.geom_poly && box.box;
    
    WITH work AS (
      SELECT * FROM grid g
      WHERE g.jdx = {j}
    ), intrsct AS (
      SELECT f.hospital_id, w.gid grid_id, w.idx idx, w.jdx jdx
      FROM for_idx f
      JOIN work w
      ON f.geom && w.geom
    )
    INSERT INTO points(idx, jdx, cnt)
    SELECT idx, jdx, count(*)
    FROM intrsct
    GROUP BY idx, jdx;
    ''')
    
    if j % 100 == 99: 
        print(f"[{datetime.now().strftime('%H:%M:%S.%f')}] finish to find intersection row # {j}")

結果がこれ…140秒くらいかかってますね…しかも答えがまちがってる…
今日はこの辺にして来週にもうちょっと続きをやります

[05:56:15.048211] starting...
Shapefile type: Point
Postgis type: POINT[2]
                    addgeometrycolumn
---------------------------------------------------------
 hospital.hospital.geom_pnt SRID:4326 TYPE:POINT DIMS:2
(1 row)

[05:56:56.017271] finish to load data
[05:57:02.788907] finish to generate buffer data
[05:57:04.709063] finished generating spatial index on polygons
[05:57:09.080523] finish to generate grid points
[05:57:20.137461] finish to find intersection row # 99
[05:57:31.610589] finish to find intersection row # 199
[05:57:42.754536] finish to find intersection row # 299
[05:57:53.571672] finish to find intersection row # 399
[05:58:05.755487] finish to find intersection row # 499
[05:58:17.325563] finish to find intersection row # 599
[05:58:29.170539] finish to find intersection row # 699
[05:58:44.853969] finish to find intersection row # 799
[05:59:04.812310] finish to find intersection row # 899
[05:59:20.750999] finish to find intersection row # 999

まとめ

PostGISのGiST空間インデックスはgeopandasのRtreeインデックスより性能が良いかもしれない
もうちょっと速くしたい（できればあと10倍くらい）
@tohka382さんの方法をPostGISでやろうとしてるがまだ未完成（できれば来週やります）