项目需求根据某些条件进行聚类统计,然后返回聚类的外轮廓,这里主要用到POSTGIS的两个算法,一个是聚类统计功能,一个是提取外轮廓的功能。
1. 聚类统计
Postgis主要实现并提供了四种聚类方法,前两个为窗口函数,后两个为聚合函数:
ST_ClusterKMeans – 该函数是窗口函数,主要是用K-means(K均值聚类)算法进行聚类,算法原理比较简单,容易实现,主要适用于点样本数据,如果是多边形则用多边形的中心点计算,只有一个参数K簇,就是事先明确了要把这一堆样本数据聚成K个类,然后去计算一个目标函数达到最优解,因此K值对结果影响比较大,不太好选取,而且该算法而且适合凸的数据集,这也算法容易收敛;
ST_ClusterDBSCAN – 该函数也是窗口函数,主要利用DBSCAN算法对输入的地理要空间素进行聚类,该算法基于密度进行聚类,直觉上更加符合认知,主要是通过地理要素分布的紧密程度决定,同一类别的样本之间是紧密相连的,不同样本是分离的。该算法相对复杂,比较常用,聚类效果较好,。
ST_ClusterIntersecting – 该函数是一个聚合函数,方法比较好理解,顾名思义就是把相交关联关系的地理要素作为一个类簇。具体的函数声明和用法,看官网的例子介绍很容易看懂,也可以拿数据进行测试实验,不是很常用,不做重点;
ST_ClusterWithin – 该函数也是一个聚合函数,也很好理解,顾名思义通过一个距离参数,将距离之内的地理要素聚成一个类别,距离之外的就是非同类。详细可参考官网,也不做重点;
综合自己的需求,采用ST_ClusterDBSCAN函数
函数定义,有三个参数,(窗口函数中的每一个地理要素,搜索半径,最小点数),根据以上分析应该不难理解:
integer ST_ClusterDBSCAN(geometry winset geom, float8 eps, integer minpoints);
例:
SELECT ST_ClusterDBSCAN(geom, eps := 0.000179, minpoints := 2) OVER () AS cluster_id,* FROM 表名
这里的eps需要根据表的srid来调整,我这边表的srid2为4490,是经纬度坐标,所以这里的0.000179也是经纬度的距离,大约20米,返回的cluster_id就是聚类后的类id,可以通过它将同一聚类的记录取出来。
2.外轮廓计算
- 计算凸包
计算凸包十分简单,使用ST_ConvexHull即可,但同时可能失去很多特征点,生成的任务区比较粗,只能是一个大致的范围。
select ST_ConvexHull ((select ST_Collect(shape) from 表名))
- 计算凹包
计算凹包十分简单相对复杂一些,使用ST_ConcaveHull,但相对于凸包会多保留一些特征点,生成的任务区保留了一些弯曲变化。
这个函数有三个参数:
- geometry geom ,几何体
- float target_percent,特征点抛弃的比例,0.1-0.99,越小特征保留的越多,但计算效率越低。当等于1时,和凸包的计算效果是一样的。
- boolean allow_holes ,是否允许有洞
select ST_ConcaveHull ((select ST_Collect(shape) from 表名),0.7)
根据需求,我们选择凹包的计算方式
3.实现需求
这里聚类和计算外轮廓分步进行计算
- 首先进行聚类查询
SELECT ST_ClusterDBSCAN(geom, eps := 0.000179, minpoints := 2) OVER () AS cluster_id,* FROM allrecord
这里allrecord我是用了with子句的查询结果,目的是实现按照项目需求条件进行查询,这里不再展示
- 在java程序中完成聚类统计
Map<Integer, List<实体类>> collect = 聚类查询结果list.stream().filter(tcql -> tcql.getCluster_id() != null).collect(Collectors.groupingBy(实体类::getCluster_id));
然后分别将上述map中的id分别取出来,后续计算外轮廓是用
- 计算外轮廓
这里直接将mapper的配置写出来
<select id="selectConcaveHull" resultType="java.util.Map">
with concaveResult as(
select ST_ConcaveHull ((select ST_Collect(geom) from 表名 where id in
<foreach collection="idList" item="id" separator="," open="(" close=")">
#{id}
</foreach>
),0.7) geom
)
select st_astext(geom) geomtext,st_area(st_transform(geom,4528)) geomarea from concaveResult
</select>
由于这里还计算了面积,所以还是采用了with子句的方式,现将凹包的geom查询出来,然后再取text形式和面积