Spatial clustering with PostGIS

Question

I'm looking for spatial clustering algorithm for using it within PostGIS-enabled database for point features. I'm going to write plpgsql function that takes distance between points within the same cluster as input. At the output function returns array of clusters. The most obvious solution is to build buffer zones specified distance around the feature and search for features into this buffer. If such features exist then continue to build a buffer around them, etc. If such features not exist that means cluster building is completed. Maybe there are some clever solutions?

Answer 1

Here is the kmeans-postgresql solution.

Installation: You need to have PostgreSQL 8.4 or greater on a POSIX host system (I wouldn't know where to start for MS Windows). If you have this installed from packages, make sure you also have the development packages (e.g., postgresql-devel for CentOS). Download and extract:

wget http://api.pgxn.org/dist/kmeans/1.0.0/kmeans-1.0.0.zip
unzip kmeans-1.0.0.zip
cd kmeans-1.0.0/

Before building, you need to set the USE_PGXS environment variable (my previous post instructed to delete this part of the Makefile, which wasn't the best of options). One of these two commands should work for your Unix shell:

# bash
export USE_PGXS=1
# csh
setenv USE_PGXS 1

Now build and install the extension:

make
make install
psql -f /usr/share/pgsql/contrib/kmeans.sql -U postgres -D postgis

(Note: I also tried this with Ubuntu 10.10, but no luck, as the path in pg_config --pgxs does not exist! This is probably an Ubuntu packaging bug)

Usage/Example: You should have a table of points somewhere (I drew a bunch of pseudo random points in QGIS). Here is an example with what I did:

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;

the 5 I provided in the second argument of the kmeans window function is the K integer to produce five clusters. You can change this to whatever integer you want.

Below is the 31 pseudo random points I drew and the five centroids with the label showing the count in each cluster. This was created using the above SQL query.

---

You can also attempt to illustrate where these clusters are with ST_MinimumBoundingCircle:

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;

Answer 2

I've written function that calculates clusters of features based on distance between them and build convex hull over this features:

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;

Example of using this function:

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

'poi' - name of layer, 'wkb_geometry' - name of geometry column, 'ogc_fid' - primary key of table, 14000 - cluster distance.

The result of using this function:

Answer 3

So far, the most promising I found is this extension for K-means clustering as a window function: http://pgxn.org/dist/kmeans/

However I haven't been able to install it successfully yet.

---

Otherwise, for basic grid clustering, you could use 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
;