Zenn
🐘

PostGISでポリゴンを分割

2022/09/06に公開
PostgreSQL
#
postgis
tech

市町村ポリゴンを取り扱うとき、1つのポリゴンではサイズが大きすぎて、複数の適度なサイズに分割したい場合もあると思います。SQL文を駆使すれば、よしなに分割することが可能です。

実行環境

  • Apple M1(macOS Monterery)
  • PostgreSQL 13.6
  • PostGIS 3.2

予めデータベースに入れておいた行政区画ポリゴンを選択します。
今回は北海道北見市を例としてますが、ここの箇所は各自環境によります。

CREATE TABLE area AS
    SELECT * FROM administrative_border_polygon WHERE name = 'kitami'

ST_GeneratePointsでポリゴン内にランダムでポイントを生成します。生成されるポイントはMultiPointのため、ST_Dumpで各ポイントを抽出します。geometryの抽出は.geomで行います。

CREATE TABLE area_pts AS
    SELECT (ST_Dump(ST_GeneratePoints(geom,2000))).geom AS geom
    FROM area


ST_ClusterKMeansでクラスタを生成します。ST_ClusterKMeansはWindow関数のためover句を付ける必要があります。

CREATE TABLE area_pts_clustered AS
    SELECT geom, ST_ClusterKMeans(geom, 20) over () AS cluster
    FROM area_pts

各クラスタのポイントをST_CollectでMultiPoint化したうえで、ST_Centroidで各クラスタの重心を求めます。

CREATE TABLE area_centers AS
    SELECT cluster, ST_Centroid(ST_collect(geom)) AS geom
    FROM area_pts_clustered
    GROUP BY cluster

ST_VoronoiPolygonsで各クラスタ重心を母点したボロノイ分割を行います。

CREATE TABLE area_volonoi AS
    SELECT (ST_Dump(ST_VoronoiPolygons(ST_collect(geom)))).geom AS geom
    FROM area_centers


ボロノイ分割したgeometryをST_Intersectionにて当初のポリゴン形状で抽出します。

CREATE TABLE area_divided AS
    SELECT ST_Intersection(a.geom, b.geom) AS geom
    FROM area a
    CROSS JOIN area_volonoi b

最後に、今までのSQL文を一つにまとめましょう。

CREATE VIEW area_divided AS 
    SELECT ST_Intersection(a.geom, b.geom) AS geom 
    FROM (
        SELECT * 
        FROM administrative_border_polygon 
        WHERE name = 'kitami'
    ) a 
    CROSS JOIN (
        SELECT (ST_Dump(ST_VoronoiPolygons(ST_collect(geom)))).geom AS geom
        FROM (
            SELECT cluster, ST_Centroid(ST_collect(geom)) AS geom
            FROM (
                SELECT geom, ST_ClusterKMeans(geom, 10) over () AS cluster
                FROM (
                    SELECT (ST_Dump(ST_GeneratePoints(geom,2000))).geom AS geom
                    FROM (
		        SELECT * FROM administrative_border_polygon
		        WHERE name = 'kitami')
                    AS area) 
                AS area_pts) 
            AS area_pts_clustered2 
            GROUP BY cluster
        )
        AS area_centers
    ) b

参考サイト
https://blog.cleverelephant.ca/2018/06/polygon-splitting.html

Discussion