基于PostGIS的高级应用（5）-- Polygon Spli
一 案例背景
PostGIS提供了丰富的function用于GIS数据的存储，元数据描述，空间分析，测量，空间图形处理等等，这些函数基本上都很简单，遇到合适的场景时，很容易能知道应该选用哪种function去解决。但有时候的图形处理问题并不是很简单就能实现的，PostGIS核心成员就遇到了社区提出的一个问题：
PostGIS是否有方法能将一个Polygon面切割成若干份小的Polygon面，且每一份的面积差不多大？
输入图形--切割前
输出图形--切割后
其第一反应是：
不可以吧，如此复杂的问题不是sql能解决的。
打脸的是，PostGIS开发者Darafei Praliaskouski解决了这个问题，并分享了解决步骤。本文作者，也就是我，仅仅负责稍微整理下，搬运了下大神们的解决方案，非个人原创。要看原文的朋友，可访问原文：《PostGIS Polygon Splitting》
二 切割步骤
Darafei Praliaskouski提供的切割步骤如下：
-
使用ST_GeneratePoints方法将一个polygon转换成与面积成比例的一系列的点 (点越多，效果越好，大约1000个点为宜)。
ST_GeneratePoints
-
假设计划将Polygon切成k等份，则使用ST_ClusterKMeans方法将这些转换后的点聚合成k簇。
ST_ClusterKMeans
-
使用ST_Centroid方法求出每一簇的的均值中心。
ST_Centroid
-
将求出的均值中心作为ST_VoronoiPolygons方法的输入参数，可以计算出每个点映射出的Polygon面。
ST_VoronoiPolygons
-
使用ST_Intersection将这些映射的面和初始化的Polygon做切割处理，得到结果。
ST_Intersection
灵活使用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;
ST_ClusterKMeans
- 获取每一簇的均值中心
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;
voronoi
- 使用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还是很强大的！！！另外一定要自己练习。。。
