一 案例背景

  PostGIS提供了豐富的function用於GIS資料的儲存，後設資料描述，空間分析，測量，空間圖形處理等等，這些函式基本上都很簡單，遇到合適的場景時，很容易能知道應該選用哪種function去解決。但有時候的圖形處理問題並不是很簡單就能實現的，PostGIS核心成員就遇到了社群提出的一個問題：

PostGIS是否有方法能將一個Polygon面切割成若干份小的Polygon面，且每一份的面積差不多大？

其第一反應是：

不可以吧，如此複雜的問題不是sql能解決的。

打臉的是，PostGIS開發者Darafei Praliaskouski解決了這個問題，並分享瞭解決步驟。本文作者，也就是我，僅僅負責稍微整理下，搬運了下大神們的解決方案，非個人原創。要看原文的朋友，可訪問原文：《PostGIS Polygon Splitting》

二 切割步驟

Darafei Praliaskouski提供的切割步驟如下：

• 使用ST_GeneratePoints方法將一個polygon轉換成與面積成比例的一系列的點 (點越多，效果越好，大約1000個點為宜)。
  
• 假設計劃將Polygon切成k等份，則使用ST_ClusterKMeans方法將這些轉換後的點聚合成k簇。
  
• 使用ST_Centroid方法求出每一簇的的均值中心。
  
• 將求出的均值中心作為ST_VoronoiPolygons方法的輸入引數，可以計算出每個點對映出的Polygon面。
  
• 使用ST_Intersection將這些對映的面和初始化的Polygon做切割處理，得到結果。
  

靈活使用PostGIS的方法，將如此複雜的問題，簡單的解決了，堪稱完美。

四 實踐總結

  百聞不如一見，百看不如一試試，本文作者就是我，看完覺得很贊，於是決定抄抄看看，如何將南京切割成大小相等的十個面，感興趣的朋友可以按照我的步驟也可以測試測試。
準備測試資料：

• 將測試資料寫入臨時表：

 create table nanjing as select name,geom from city where name=`南京市`;

• 面轉換為點：

CREATE TABLE nanjing_points AS SELECT (ST_Dump(ST_GeneratePoints(geom, 2000))).geom 
AS geom FROM nanjing;

• 點聚合成簇（看原文方法的朋友請注意他的ST_ClusterKMeans寫錯了）

CREATE TABLE nanjing_pts_clustered AS 
SELECT geom, ST_ClusterKMeans(geom, 10) over () AS cluster FROM nanjing_points;

• 獲取每一簇的均值中心

CREATE TABLE nanjing_centers AS SELECT cluster, ST_Centroid(ST_collect(geom)) AS geom
  FROM nanjing_pts_clustered GROUP BY cluster;

• 使用voronoi演算法生成面

  CREATE TABLE nanjing_voronoi AS
  SELECT (ST_Dump(ST_VoronoiPolygons(ST_collect(geom)))).geom AS geom
  FROM nanjing_centers;

• 使用ST_Intersection方法切割

CREATE TABLE nanjing_divided AS
  SELECT ST_Intersection(a.geom, b.geom) AS geom
  FROM nanjing a
  CROSS JOIN nanjing_voronoi b;

一個個寫sql生成的臨時表寫邏輯方便，但是使用起來比較費勁，我們可以寫個function去處理，我使用了臨時表，怕事務併發衝突，加了uuid字尾。其實這個邏輯不是很複雜的話，多套用幾個with中間表也可以，但是寫多了不是很清晰，我就暫時套用上面表的邏輯改成臨時表做了個事務，測試通過了：

create extension "uuid-ossp";--建立下uuid的擴充套件

create or replace function freegis_polygon_split(
    in split_geom geometry(Polygon),--輸入的面
    in split_num int,--分割的數量
    out geom geometry(Polygon)--輸出切割的面
) returns setof geometry as 
$$
  
declare  
    rec record;
    temp_points text;
    temp_ClusterKMeans text;
    temp_ClusterCentroid text;
    temp_VoronoiPolygons text;
    
begin  
    --防止併發的時候，臨時表名稱衝突
    temp_points:=`temp_points`||uuid_generate_v4();
    temp_ClusterKMeans:=`temp_ClusterKMeans`||uuid_generate_v4();
    temp_ClusterCentroid:=`temp_ClusterCentroid`||uuid_generate_v4();
    temp_VoronoiPolygons:=`temp_VoronoiPolygons`||uuid_generate_v4();

    --生成點
    execute format(`create temp table "%s" on commit drop as 
    SELECT row_number() over() as gid,(ST_Dump(ST_GeneratePoints($1, 2000))).geom`,temp_points) using split_geom;
    --點成簇
    execute format(`create temp table "%s" on commit drop as SELECT t.geom, ST_ClusterKMeans(t.geom, $1) over () AS cluster from "%s" t`,temp_ClusterKMeans,temp_points) using split_num;
    --簇的中心點
    execute format(`create temp table "%s" on commit drop as SELECT t.cluster, ST_Centroid(ST_collect(t.geom)) AS geom FROM "%s" t GROUP BY t.cluster`,temp_ClusterCentroid,temp_ClusterKMeans);
    --voronoi構造面
    execute format(`create temp table "%s" on commit drop as SELECT (ST_Dump(ST_VoronoiPolygons(ST_collect(t.geom)))).geom AS geom FROM "%s" t`,temp_VoronoiPolygons,temp_ClusterCentroid);
    --intersection切割
    for rec in execute format(`SELECT ST_Intersection($1, b.geom) AS geom FROM "%s" b`,temp_VoronoiPolygons) using split_geom loop
        geom:=rec.geom;
        return next;
    end loop;
    return;
end;  

$$
 language plpgsql strict;  

測試下通過了：

結語：PostGIS還是很強大的！！！另外一定要自己練習。。。