実行環境
Windows10
Python 3.9.5
shapely 1.8.0
geopandas 0.10.2
はじめに
ボロノイ分割とは「平面上の複数の地点の最近隣領域を求める」処理のことです。（ざっくり）
詳しくはこちら等をご覧ください。
・Wikipedia ボロノイ図
・PASCO様の用語集 ボロノイ分割とは？
個人的にはティーセン分割という名前の方がなじみがあるのですが、幾何学の分野ではボロノイ分割というようです。
ここでは、ライブラリ shapely を使ったボロノイ分割の方法について説明します。
入出力ともシェープファイルを想定しています。
ボロノイ分割は scipy 1 を使った実装の方が一般的かと思いますが、GISデータの処理を前提とする場合は shapely を使うとより簡単です。
shapely のバージョン
shapely のバージョンは 現時点最新安定版の 1.8.0 を使います。
1.7.1 以前のバージョンには目的の関数 voronoi_diagram() が含まれていないため、気をつけてください。
テストケース
「どの避難施設が最も近いか」が分かる地図を作ります。
本来はしっかりネットワーク解析等して作るべきですが、今回はボロノイ分割をしたいので、単純に直線距離で最近隣を判定するということにします。
避難施設のポイントデータは 国土数値情報の避難施設データ を使います。
対象都道府県は、最もデータが軽い滋賀県としました。
コード
避難施設のポイントデータをもとに、ボロノイ分割されたポリゴンデータを出力するコードです。
import geopandas as gpd
from shapely.geometry import MultiPoint
from shapely.ops import voronoi_diagram
# ポイントデータのファイルパス
fp_in = r'P20-12_41_GML\P20-12_41.shp'
# 出力ファイルパス
fp_out = r'voronoi_P20-12_41.shp'
# ポイントデータの読込と投影変換
gdf_point = gpd.read_file(fp_in).to_crs(2444)
# 座標を MultiPoint に格納
points = MultiPoint(
[(point.x, point.y) for point in gdf_point['geometry']]
)
# ボロノイ分割 *後述
regions = voronoi_diagram(points)
# ポリゴンを GeoDataFrame に格納
gdf_voronoi = gpd.GeoDataFrame(
regions, columns=['geometry'], crs=gdf_point.crs
)
# ポイントの属性をポリゴンに付与（空間結合）
gdf_voronoi = gpd.sjoin(gdf_voronoi, gdf_point)
# シェープファイルに出力
gdf_voronoi.to_file(fp_out, encoding='utf-8')
voronoi_diagram() の構文
指定のジオメトリオブジェクトからボロノイ図形を作成する関数です。
ユーザマニュアルはこちらです。
引数は 4 つあります。
shapely.ops.voronoi_diagram(
geom,
envelope=None,
tolerance=0.0,
edges=False
)
voronoi_diagram の引数
引数についての簡単な説明です。
すみませんが、やや曖昧な部分があります。
詳細はユーザマニュアル等をご確認ください。
geom
ボロノイ領域のもととなる地点を示す shapely のジオメトリオブジェクトです。
今回のテストケースでは避難施設を MultiPoint として使用しましたが、他のジオメトリタイプでも有効です。
例えば LineString を指定した場合、全ての頂点が使用されます。
envelope
ボロノイ分割する領域の外側境界を示す shapely のジオメトリオブジェクトです。
指定しない場合は、自動で決定されます。
私が確認した限りでは、自動で決定される外側境界のさらに外側に指定した場合のみ有効でした。
tolerance
「計算の堅ろう性を改善するために使用されるスナップ許容値」として指定できるようです。
指定することでエラーが発生することもあるそうなので、使う場合はご注意ください。
edges
出力のジオメトリタイプをラインにしたいときに指定します。
デフォルトは False で、ボロノイ図形がポリゴンで出力されます。
True を指定するとバラバラの線分ラインが出力されます。
なお、先述のコードは edges=True には対応していません。
結果
QGIS で避難施設のポイントデータとボロノイ分割されたポリゴンデータを表示しています。
背景には地理院タイルの標準地図を使用しています。
全体図
佐賀駅周辺
各避難施設の領域が、過不足なく敷き詰められていることが分かります。
また、各避難施設から等距離となる境界で、意図通りに領域が分割されているらしいことも確認できます。
注意点
地理座標系の座標単位は度なので、距離を直接求めることはできません。
国土数値情報のデータは地理座標系（世界測地系 JGD2000）のため、距離を計算できる投影座標系に変換する必要があります。
ここでは、投影座標系（平面直角座標系第2系、EPSGコード2444）に投影変換しています。
さいごに
shapely を使ったボロノイ分割の方法について説明しました。
voronoi_diagram() にマルチポイントを渡すだけなので簡単ですね。