目的

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

アルゴリズムはいくつか考えられますが、各地点に対してポリゴンがいくつ重なっているかを分析する方法が思い浮かびます。ただ困ったことに分析地点が $M$ 地点、ポリゴンが $N$ 個あったとき、交差有無の確認を $MN$ 回行う必要があり、計算範囲が広くまたポリゴン数が大量の場合にはとんでもない時間が必要となります。

デモとして、下記の処理を行ってみました。

サンプルデータとして
- 国土数値情報の全国の医療機関データ（令和2年版）
  - 経緯度座標系の JGD2011 (EPSG:6668)
  - 点形式
  - 地物数は181,312個
- この点データに対し、$0.05^\circ$のバッファをかけたポリゴンを用意
  - 特に意味はないけど EPSG:4326 にして処理
  - サンプルデータ作成のためなので、経緯度座標系のままなのは承知のうえ
デモとして、分析は
- $(137^\circ,\ 35^\circ)$ - $(138^\circ,\ 36^\circ)$ の範囲
- 解像度（1セルあたりの距離）$0.001^\circ$間隔で計算
- つまり横1,000列、縦1,000行、合計1,000,000点で分析
結果はラスタ形式で出力

方針1　何も考えない

力こそパワー！

とりあえず、思い浮かぶ方法で愚直に計算してみます。

#!/usr/bin/python

from datetime import datetime
import numpy as np
import geopandas as gpd
from osgeo import gdal, osr
from shapely.geometry import Point

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

# サンプルデータの読み込み
# 国土数値情報の医療機関データ（点）
points = gpd.read_file(r'./P04-20_GML/P04-20.shp', encoding='cp932')
points.to_crs(crs=4326, inplace=True)

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

# 点データをもとに 0.05 deg ~ 5 km バッファを生成
# 地理座標系（経緯度座標系）なので Warning がでるけど、意図的なので無視
# デモ用のデータだし
regions = gpd.GeoDataFrame(geometry=points.buffer(distance=0.05), crs=4326)

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)



# 対象ポリゴンを全域のポリゴンに
target_regions = regions

# 個数を管理するための numpy の array
raster_data = np.zeros((width, height), dtype=np.uint32)
for j in range(height):
    for i in range(width):
        # セル中心の座標のため 0.5 を加える
        coord = (bounds[0] + (i+0.5)*resol, bounds[3] - (j+0.5)*resol)

        # 空間結合し、マッチした行数を拾う
        pt = gpd.GeoDataFrame({'geometry': [Point(coord)]}, crs=4326)
        raster_data[j][i] = len(gpd.sjoin(target_regions, pt, how='inner', predicate='intersects'))
    print(f"[{datetime.now().strftime('%H:%M:%S.%f')}] finish to calculate row #{j}")



# 結果をラスタファイルに出力する
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

結果は、以下のとおりです。途中で経緯度座標系に対しバッファをあててるところで Warning が出ていますが、今回は意図的にやってるので無視します。

[13:06:42.266684] starting...
[13:06:49.325854] finish to load data
test_1.py:21: UserWarning: Geometry is in a geographic CRS. Results from 'buffer' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  regions = gpd.GeoDataFrame(geometry=points.buffer(distance=0.05), crs=4326)
[13:06:51.889366] finish to generate buffer data
[13:07:47.566205] finish to calculate row #0
[13:08:43.549834] finish to calculate row #1
[13:09:40.043851] finish to calculate row #2
    :
    :

うちのマシンで、1行（1,000点）の処理で1分弱かかっています。単純計算だと16時間とかそのくらいの時間が必要になり、さすがにきついです。

方針2　分割する

大きな問題は、小さな問題に分割せよ

というわけで、大量のデータを扱う際のセオリーに立って、分割統治法を考えてみます。

九州の分析地点において、北海道のポリゴンと交差する可能性はゼロです。つまり地理的に分割し、別々に計算してもまったく問題ありません。不要なポリゴンを削除し $N$ を小さくすること計算量を下げる作戦です。最初にポリゴンを必要な範囲だけ抽出する処理は必要ですが、そのあとの各点での1,000,000回の分析において効いてきます。

先ほどのスクリプトを

-# 対象ポリゴンを全域のポリゴンに
-target_regions = regions
+# 先に範囲を絞る
+target_regions = regions.clip(bounds)

にしてみます。結果は、

[13:31:54.816383] starting...
[13:32:01.867153] finish to load data
test_2.py:21: UserWarning: Geometry is in a geographic CRS. Results from 'buffer' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  regions = gpd.GeoDataFrame(geometry=points.buffer(distance=0.05), crs=4326)
[13:32:04.402002] finish to generate buffer data
[13:32:18.087468] finish to calculate row #0
[13:32:31.779864] finish to calculate row #1
[13:32:45.658077] finish to calculate row #2
    :
    :

1行（1,000点）あたり14秒くらいなので、4倍速くなりました。もともとのデータが日本全国であり、フィルタリングした範囲は1平方度の限られた範囲であったため、もっと速くなると思ったのですが、思ったより速くはなりませんでした。

もちろん、もっと小範囲に区切ると速くなるはずですが、区切りすぎると分割する処理が増えてしまいます。

方針3　空間インデックスを利用し、対象を絞る

――まずは自己紹介をしなくちゃいけないね
　　　　　　私の名前はね、インデックスって言うんだよ？――

空間的な分割統治法は明らかに不要な部分を除外する方法でした。一方 GIS データには空間インデックスという、事前に地物がどのあたりにあるか大まかにマッピングする仕組みがあり、空間インデックスを利用すると、そもそも不要な箇所を分割する必要すらないはずです。

ただし空間インデックスは地物の細かな形状までチェックしているわけではなく、エンベロープに対してインデックス化されているため、一次段階の簡易判定とし、実際に交差するかは別途確認する必要があります。

# 先に範囲を絞る
target_regions = regions.clip(bounds)
# 空間インデックスを作成
regions_sindex = target_regions.sindex

# 個数を管理するための numpy の array
raster_data = np.zeros((width, height), dtype=np.uint32)
for j in range(height):
    for i in range(width):
        # セル中心の座標のため 0.5 を加える
        coord = (bounds[0] + (i+0.5)*resol, bounds[3] - (j+0.5)*resol)

        pt = gpd.GeoDataFrame({'geometry': [Point(coord)]}, crs=4326)
        # 空間インデックスを利用して、エンベロープに交差判定がでたポリゴンを抽出する
        candidated_regions = target_regions.iloc[regions_sindex.intersection(coord), :]
        # ざっくりと抽出したポリゴンに対し、詳細な交差のチェックを行い数を求める
        raster_data[j][i] = len(gpd.sjoin(candidated_regions, pt, how='inner', predicate='intersects'))
    print(f"[{datetime.now().strftime('%H:%M:%S.%f')}] finish to calculate row #{j}")

結果は、以下となりました。

[16:39:32.489957] starting...
[16:39:39.520795] finish to load data
test_3.py:21: UserWarning: Geometry is in a geographic CRS. Results from 'buffer' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  regions = gpd.GeoDataFrame(geometry=points.buffer(distance=0.05), crs=4326)
[16:39:42.054760] finish to generate buffer data
[16:39:55.284986] finish to calculate row #0
[16:40:08.613861] finish to calculate row #1
[16:40:21.695760] finish to calculate row #2
    :
    :

あ、あれ。全然変わってない。交差判定をかけるポリゴンはかなり絞られているはずなのになんでだ……。

pandas (geopandas) の取り扱い方とかがよくわかってなくて、ちゃんと絞り切れていないからでしょうか？？？

方針4　空間インデックスで交差判定

こまけぇこたぁいいんだよ!!

空間インデックスを使った交差判定はエンベロープに対するものなので、交差しているという判定が出たからといって、セル中心点に地物本体が存在しているかはわかりません。わかりませんが、ある程度近くに存在していることは確かです。ちなみに、空間インデックスだけで交差判定を行って集計すると、こんな感じになります。

とはいえ、そのまま $(137^\circ,\ 35^\circ)$ - $(138^\circ,\ 36^\circ)$ の範囲全体で空間インデックスをとり、上記のような結果を使用するのはさすがにやりすぎなので、行ごとに細長い短冊を作成することにします。その短冊に対し空間インデックスを作成し、空間インデックスとして交差判定を行ってみます。

下図は細い短冊状に分割した際のイメージ図となります。セル及びセル中心があり、それに重なるように緑のポリゴンがあります。セル中心点付近に細い短冊状の範囲（赤枠）を設定し、ポリゴンを切り取ります。その切り取られたポリゴンに対し、空間インデックスを張るとポリゴン片のエンベロープ（赤塗り潰し）がインデックス化されます。空間インデックス（ポリゴン片のエンベロープ）に対して交差を判定すると、本来のポリゴン領域と重ならない場合でも交差していると判断されるケースがでてきます。

セル中心点において正確に重なっているかは保証されませんが、非常に細い短冊領域とは多少なりとも重複しているはずなので、まあだいたい正確かと思います。

# 個数を管理するための numpy の array
raster_data = np.zeros((width, height), dtype=np.uint32)
for j in range(height):
    # 当該行において、セル中心点 (+0.5) を含む範囲で細長い短冊をつくる
    target_regions = regions.clip((bounds[0], bounds[3] - (j+0.6)*resol, bounds[2], bounds[3] - (j+0.4)*resol))
    # 空間インデックスを計算する
    regions_sindex = target_regions.sindex

    for i in range(width):
        # セル中心の座標のため 0.5 を加える
        coord = (bounds[0] + (i+0.5)*resol, bounds[3] - (j+0.5)*resol)

        # 空間インデックスを用いて、大まかに重なり数を求める
        raster_data[j][i] = len(regions_sindex.intersection(coord))
    print(f"[{datetime.now().strftime('%H:%M:%S.%f')}] finish to calculate row #{j}")

結果は

[17:05:19.066477] starting...
[17:05:26.080451] finish to load data
test_4.py:21: UserWarning: Geometry is in a geographic CRS. Results from 'buffer' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  regions = gpd.GeoDataFrame(geometry=points.buffer(distance=0.05), crs=4326)
[17:05:28.624656] finish to generate buffer data
[17:05:28.763601] finish to calculate row #0
[17:05:28.798753] finish to calculate row #1
[17:05:28.836869] finish to calculate row #2
    :
    :
[17:06:10.819567] finish to calculate row #997
[17:06:10.866994] finish to calculate row #998
[17:06:10.913629] finish to calculate row #999

対象領域全体の交差判定が42秒ほどで完了しました。

QGIS に読み込んで、適当に着色するとこんな感じです。ポリゴンの境界付近が正確に求められているかはわかりませんが、おおまかにはちゃんと動いているように見えます。