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?
There are at least two good clustering methods for PostGIS: k-means (via
kmeans-postgresql extension) or clustering geometries within a threshold distance (PostGIS 2.2)
1) k-means with
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.1.0/kmeans-1.1.0.zip unzip kmeans-1.1.0.zip cd kmeans-1.1.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;
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;
2) Clustering within a threshold distance with
This aggregate function is included with PostGIS 2.2, and returns an array of GeometryCollections where all the components are within a distance of each other.
Here is an example use, where a distance of 100.0 is the threshold that results in 5 different clusters:
SELECT row_number() over () AS id, ST_NumGeometries(gc), gc AS geom_collection, ST_Centroid(gc) AS centroid, ST_MinimumBoundingCircle(gc) AS circle, sqrt(ST_Area(ST_MinimumBoundingCircle(gc)) / pi()) AS radius FROM ( SELECT unnest(ST_ClusterWithin(geom, 100)) gc FROM rand_point ) f;
The largest middle cluster has a enclosing circle radius of 65.3 units or about 130, which is larger than the threshold. This is because the individual distances between the member geometries is less than the threshold, so it ties it together as one larger cluster.
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:
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 ;
Complementing @MikeT answer...
For MS Windows:
What you will do:
- Tweak the source code to export the kmeans function to a DLL.
- Compile the source code with
cl.execompiler to generate a DLL with
- Put the generated DLL into PostgreSQL\lib folder.
- Then you can "create" (link) the UDF into PostgreSQL through SQL command.
- Download & install/extract requirements.
kmeans.cin any editor:
#includelines define DLLEXPORT macro with:
#if defined(_WIN32) #define DLLEXPORT __declspec(dllexport) #else #define DLLEXPORT #endif
DLLEXPORTbefore each of these lines:
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);
Open Visual C++ Command Line.
In the command line:
- Go to the extracted
- Set your POSTGRESPATH, mine for example is:
SET POSTGRESPATH=C:\Program Files\PostgreSQL\9.5
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"
- Go to the extracted
Now run the SQL command in your database to "CREATE" the function.
CREATE FUNCTION kmeans(float, int) RETURNS int AS '$libdir/kmeans' LANGUAGE c VOLATILE STRICT WINDOW; CREATE FUNCTION kmeans(float, int, float) RETURNS int AS '$libdir/kmeans', 'kmeans_with_init' LANGUAGE C IMMUTABLE STRICT WINDOW;
Here is a way to display in QGIS the result of the PostGIS query given in 2) in this anwser
As QGIS doesn't handle neither geometrycollections nor different datatypes in the same geometry column, I've created two layers, one for clusters and one for clustered points.
First for clusters, you only need polygons, other results are lonely points :
SELECT id,countfeature,circle FROM (SELECT row_number() over () AS id, ST_NumGeometries(gc) as countfeature, ST_MinimumBoundingCircle(gc) AS circle FROM ( SELECT unnest(ST_ClusterWithin(the_geom, 100)) gc FROM rand_point ) f) a WHERE ST_GeometryType(circle) = 'ST_Polygon'
Then for clustered points, you need to transform geometrycollections in multipoint :
SELECT row_number() over () AS id, ST_NumGeometries(gc) as countfeature, ST_CollectionExtract(gc,1) AS multipoint FROM ( SELECT unnest(ST_ClusterWithin(the_geom, 100)) gc FROM rand_point ) f
Some points are at same coordinates so the label could be confusing.
You can use Kmeans solution more easily with ST_ClusterKMeans method that's available in postgis from 2.3 Example:
SELECT kmean, count(*), ST_SetSRID(ST_Extent(geom), 4326) as bbox FROM ( SELECT ST_ClusterKMeans(geom, 20) OVER() AS kmean, ST_Centroid(geom) as geom FROM sls_product ) tsub GROUP BY kmean;
The bounding box of features is used as cluster geometry in the example above. The first image shows the original geometries and the second one is the result of select above.
Bottom up clustering solution from Get a single cluster from cloud of points with maximum diameter in postgis which involves no dynamic queries.
CREATE TYPE pt AS ( gid character varying(32), the_geom geometry(Point))
and a type with cluster id
CREATE TYPE clustered_pt AS ( gid character varying(32), the_geom geometry(Point) cluster_id int)
Next the algorithm function
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 WHERE cluster_id = joined_clusters; --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
WITH subq AS( SELECT ARRAY_AGG((gid, the_geom)::pt) AS points FROM data GROUP BY collection_id) SELECT (clusters).* FROM (SELECT buc(points, radius) AS clusters FROM subq ) y;