使用PostGIS进行空间聚类？

#1 楼

PostGIS至少有两种良好的聚类方法：k均值（通过kmeans-postgresql扩展）或阈值距离内的聚类几何形状（PostGIS 2.2）


1）具有kmeans-postgresql


安装：您需要从源代码进行编译和安装，这在* NIX上比Windows容易（我不知道从哪里开始）。如果从软件包中安装了PostgreSQL，请确保您也具有开发软件包（例如，对于CentOS为postgresql-devel）。

下载，提取，构建和安装： />

wget http://api.pgxn.org/dist/kmeans/1.1.0/kmeans-1.1.0.zip
unzip kmeans-1.1.0.zip
cd kmeans-1.1.0/
make USE_PGXS=1
sudo make install

启用数据库中的扩展名（使用psql，pgAdmin等）：

CREATE EXTENSION kmeans;

用法/示例：您应该在某处有一个点表（我在QGIS中绘制了一堆伪随机点）。以下是我所做的示例：集群。您可以将其更改为所需的任何整数。

下面是我绘制的31个伪随机点和五个质心，其中的标签显示每个群集的计数。这是使用上述SQL查询创建的。




您还可以尝试使用ST_MinimumBoundingCircle来说明这些群集的位置： >

SELECT kmeans, count(*), ST_Centroid(ST_Collect(geom)) AS geom
FROM (
  SELECT kmeans(ARRAY[ST_X(geom), ST_Y(geom)], 5) OVER (), geom
  FROM rand_point
) AS ksub
GROUP BY kmeans
ORDER BY kmeans;

2）使用5在阈值距离内聚类


此聚合函数包含在PostGIS 2.2中，并返回一个GeometryCollections数组，其中所有组件彼此之间的距离都在此范围内。

这里是一个示例用法，其中100.0的距离是导致5个不同的群集：

SELECT kmeans, ST_MinimumBoundingCircle(ST_Collect(geom)) AS circle
FROM (
  SELECT kmeans(ARRAY[ST_X(geom), ST_Y(geom)], 5) OVER (), geom
  FROM rand_point
) AS ksub
GROUP BY kmeans
ORDER BY kmeans;

最大的中间簇的包围圆半径为65.3单位或大约130，大于阈值。这是因为成员几何之间的单个距离小于阈值，因此将其绑定为一个较大的簇。

#2 楼

我编写了一个函数，该函数根据特征之间的距离计算特征簇，并在此特征上构建凸包：

CREATE OR REPLACE FUNCTION get_domains_n(lname varchar, geom varchar, gid varchar, radius numeric)
    RETURNS SETOF record AS
$$
DECLARE
    lid_new    integer;
    dmn_number integer := 1;
    outr       record;
    innr       record;
    r          record;
BEGIN

    DROP TABLE IF EXISTS tmp;
    EXECUTE 'CREATE TEMPORARY TABLE tmp AS SELECT '||gid||', '||geom||' FROM '||lname;
    ALTER TABLE tmp ADD COLUMN dmn integer;
    ALTER TABLE tmp ADD COLUMN chk boolean DEFAULT FALSE;
    EXECUTE 'UPDATE tmp SET dmn = '||dmn_number||', chk = FALSE WHERE '||gid||' = (SELECT MIN('||gid||') FROM tmp)';

    LOOP
        LOOP
            FOR outr IN EXECUTE 'SELECT '||gid||' AS gid, '||geom||' AS geom FROM tmp WHERE dmn = '||dmn_number||' AND NOT chk' LOOP
                FOR innr IN EXECUTE 'SELECT '||gid||' AS gid, '||geom||' AS geom FROM tmp WHERE dmn IS NULL' LOOP
                    IF ST_DWithin(ST_Transform(ST_SetSRID(outr.geom, 4326), 3785), ST_Transform(ST_SetSRID(innr.geom, 4326), 3785), radius) THEN
                    --IF ST_DWithin(outr.geom, innr.geom, radius) THEN
                        EXECUTE 'UPDATE tmp SET dmn = '||dmn_number||', chk = FALSE WHERE '||gid||' = '||innr.gid;
                    END IF;
                END LOOP;
                EXECUTE 'UPDATE tmp SET chk = TRUE WHERE '||gid||' = '||outr.gid;
            END LOOP;
            SELECT INTO r dmn FROM tmp WHERE dmn = dmn_number AND NOT chk LIMIT 1;
            EXIT WHEN NOT FOUND;
       END LOOP;
       SELECT INTO r dmn FROM tmp WHERE dmn IS NULL LIMIT 1;
       IF FOUND THEN
           dmn_number := dmn_number + 1;
           EXECUTE 'UPDATE tmp SET dmn = '||dmn_number||', chk = FALSE WHERE '||gid||' = (SELECT MIN('||gid||') FROM tmp WHERE dmn IS NULL LIMIT 1)';
       ELSE
           EXIT;
       END IF;
    END LOOP;

    RETURN QUERY EXECUTE 'SELECT ST_ConvexHull(ST_Collect('||geom||')) FROM tmp GROUP by dmn';

    RETURN;
END
$$
LANGUAGE plpgsql;

使用此函数的示例：

SELECT * FROM get_domains_n('poi', 'wkb_geometry', 'ogc_fid', 14000) AS g(gm geometry)

'poi'-图层名称，'wkb_geometry'-几何列名称，'ogc_fid'-表的主键，14000-簇距离。

使用此功能的结果：

#3 楼

到目前为止，我发现最有前途的是作为窗口功能的K-means集群扩展：http://pgxn.org/dist/kmeans/



否则，对于基本的网格群集，您可以使用SnapToGrid。

SELECT
    array_agg(id) AS ids,
    COUNT( position ) AS count,
    ST_AsText( ST_Centroid(ST_Collect( position )) ) AS center,
FROM mytable
GROUP BY
    ST_SnapToGrid( ST_SetSRID(position, 4326), 22.25, 11.125)
ORDER BY
    count DESC
;

#4 楼

您可以通过2.3中的postgis中可用的ST_ClusterKMeans方法更轻松地使用Kmeans解决方案。示例：
第一个图像显示的是原始几何形状，第二个图像是上述选择的结果。

#5 楼

补充@MikeT答案...

对于MS Windows：

要求：


任何像这样的Visual C ++ Express版本

kmeans-postgresql模块。

您将做什么：


调整源代码以将kmeans函数导出到DLL。
使用cl.exe编译器编译源代码以使用kmeans函数生成DLL。
将生成的DLL放入PostgreSQL \ lib文件夹中。
然后您可以“创建”（链接）UDF到PostgreSQL通过SQL命令。

步骤：


下载并安装/提取要求。在kmeans.c行之后，使用以下命令定义DLLEXPORT宏：

#if defined(_WIN32)
    #define DLLEXPORT __declspec(dllexport)
#else
   #define DLLEXPORT
#endif

将#include放在每一个之前行：

PG_FUNCTION_INFO_V1(kmeans_with_init);
PG_FUNCTION_INFO_V1(kmeans);

extern Datum kmeans_with_init(PG_FUNCTION_ARGS);
extern Datum kmeans(PG_FUNCTION_ARGS);

打开Visual C ++命令行。

在命令行中：


转到提取的DLLEXPORT。
设置我的POSTGRESPATH例如：kmeans-postgresql


运行

cl.exe /I"%POSTGRESPATH%\include" /I"%POSTGRESPATH%\include\server" /I"%POSTGRESPATH%\include\server\port\win32" /I"%POSTGRESPATH%\include\server\port\win32_msvc" /I"C:\Program Files (x86)\Microsoft SDKs\Windows\v7.1A\Include" /LD kmeans.c "%POSTGRESPATH%\lib\postgres.lib"

将SET POSTGRESPATH=C:\Program Files\PostgreSQL.5复制到kmeans.dll

现在在数据库中运行SQL命令以“创建”功能。

#6 楼

这是一种在QGIS中显示此anwser中2）中给出的PostGIS查询结果的方法。两层，一层用于聚类，一层用于聚类点。

首先用于聚类，您只需要多边形，其他结果是孤独的点：
/>然后，对于聚类点，您需要在多点中转换几何体集合： br />

#7 楼

自底向上的聚类解决方案，从postgis中从具有最大直径的点云中获取单个聚类，不涉及任何动态查询。 >

CREATE TYPE pt AS (
    gid character varying(32),
    the_geom geometry(Point))

下一个算法函数

CREATE TYPE clustered_pt AS (
    gid character varying(32),
    the_geom geometry(Point)
    cluster_id int)

用法：

CREATE OR REPLACE FUNCTION buc(points pt[], radius integer)
RETURNS SETOF clustered_pt AS
$BODY$

DECLARE
    srid int;
    joined_clusters int[];

BEGIN

--If there's only 1 point, don't bother with the loop.
IF array_length(points,1)<2 THEN
    RETURN QUERY SELECT gid, the_geom, 1 FROM unnest(points);
    RETURN;
END IF;

CREATE TEMPORARY TABLE IF NOT EXISTS points2 (LIKE pt) ON COMMIT DROP;

BEGIN
    ALTER TABLE points2 ADD COLUMN cluster_id serial;
EXCEPTION
    WHEN duplicate_column THEN --do nothing. Exception comes up when using this function multiple times
END;

TRUNCATE points2;
    --inserting points in
INSERT INTO points2(gid, the_geom)
    (SELECT (unnest(points)).* ); 

--Store the srid to reconvert points after, assumes all points have the same SRID
srid := ST_SRID(the_geom) FROM points2 LIMIT 1;

UPDATE points2 --transforming points to a UTM coordinate system so distances will be calculated in meters.
SET the_geom =  ST_TRANSFORM(the_geom,26986);

--Adding spatial index
CREATE INDEX points_index
ON points2
USING gist
(the_geom);

ANALYZE points2;

LOOP
    --If the smallest maximum distance between two clusters is greater than 2x the desired cluster radius, then there are no more clusters to be formed
    IF (SELECT ST_MaxDistance(ST_Collect(a.the_geom),ST_Collect(b.the_geom))  FROM points2 a, points2 b
        WHERE a.cluster_id <> b.cluster_id
        GROUP BY a.cluster_id, b.cluster_id 
        ORDER BY ST_MaxDistance(ST_Collect(a.the_geom),ST_Collect(b.the_geom)) LIMIT 1)
        > 2 * radius
    THEN
        EXIT;
    END IF;

    joined_clusters := ARRAY[a.cluster_id,b.cluster_id]
        FROM points2 a, points2 b
        WHERE a.cluster_id <> b.cluster_id
        GROUP BY a.cluster_id, b.cluster_id
        ORDER BY ST_MaxDistance(ST_Collect(a.the_geom),ST_Collect(b.the_geom)) 
        LIMIT 1;

    UPDATE points2
    SET cluster_id = joined_clusters[1]
    WHERE cluster_id = joined_clusters[2];

    --If there's only 1 cluster left, exit loop
    IF (SELECT COUNT(DISTINCT cluster_id) FROM points2) < 2 THEN
        EXIT;

    END IF;

END LOOP;

RETURN QUERY SELECT gid, ST_TRANSFORM(the_geom, srid)::geometry(point), cluster_id FROM points2;
END;
$BODY$
LANGUAGE plpgsql