美文网首页PostgreSQL开源开源GIS+空间数据应用
基于PostGIS的高级应用(5)-- Polygon Spli

基于PostGIS的高级应用(5)-- Polygon Spli

作者: 遥想公瑾当年 | 来源:发表于2018-09-25 16:33 被阅读347次

    一 案例背景

      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还是很强大的!!!另外一定要自己练习。。。

    相关文章

      网友评论

      • li_9569:请问一下,作者是的什么方法把sql查询结果可视化的?😀
        li_9569:@遥想公瑾当年 噢,谢谢啦
        遥想公瑾当年:@li_9569 一般arcgis直连pg,或者pg导出shp。但主键和字段类型严格。

      本文标题:基于PostGIS的高级应用(5)-- Polygon Spli

      本文链接:https://www.haomeiwen.com/subject/duqloftx.html