PostGIS Recursively Finding Intersections of Polygons to Determine Clusters

Question

My dataset currently contains about 15.000 buildings which partially touch each other. The objective is that those which do touch should belong to a cluster. This "basically" works but the main issue I am facing is the speed which is taking hours and if the dataset of polygons is too large it feels like the recursion is entering an endless loop.

-- generating the clusters table which should feature the row number and an array of building ids which make up the cluster
CREATE TABLE clusters (
   rownum INTEGER NOT NULL,
   ret INTEGER[]
);

-- recursion to find the polygons which are very close to each other by buffering them and finding checking if they intersect
WITH RECURSIVE 
source (rownum, geom, ret) AS (
   SELECT row_number() OVER (ORDER BY id ASC), ST_Multi(geom), ARRAY[id::numeric::integer] FROM buildings), 
result (rownum, geom, ret, incroci) AS (
    SELECT rownum, geom, ret, 1 FROM source 
    UNION ALL
    SELECT s.rownum, s.geom, (r.ret || s.ret), (r.incroci + 1) 
    FROM source AS s 
    INNER JOIN result as r
    ON s.rownum > r.rownum 
    AND ST_Intersects(ST_Buffer(s.geom,0.2), ST_Buffer(r.geom,0.2))
)
-- ultimately insert them to the clusters table if they contain more than 1 building
INSERT INTO clusters(rownum, ret)
SELECT row_number() OVER (ORDER BY rownum ASC), ret FROM result
WHERE array_length(ret, 1) > 1;

-- we create an index on the clusters table array column which will be needed later
CREATE INDEX clusters_idx_ret_with_intarray ON clusters USING GIN(ret gin__int_ops);
CREATE INDEX clusters_idx_rownum ON clusters (rownum);

So basically this script is recursively finding polygons which intersect with each other and afterwards inserts the results into a new table. Currently I am using a subset of the 15.000 buildings of about 1.429 buildings which takes approximately 3 minutes (however if I use a subset of 5.000 buildings it runs for hours and never stops).

The result will look something like this:

gis=# SELECT * FROM clusters ORDER BY ret;
 rownum |              ret
--------+--------------------------------
      8 | {3,2262}
     37 | {3,2262,3456}
    209 | {3,2262,3456,5827}
    451 | {3,2262,3456,5827,7653}
    456 | {3,2262,3456,7653}
    529 | {3,2262,8144}
    570 | {3,2262,8144,8437}
    682 | {3,2262,8144,8913}
   1089 | {3,2262,8144,8913,10436}
   1086 | {3,2262,8144,10436}
    566 | {3,2262,8437}
     38 | {3,3456}
    207 | {3,3456,5827}
    457 | {3,3456,5827,7653}
    455 | {3,3456,7653}
    208 | {3,5827}
    453 | {3,5827,7653}
    955 | {9,9988}
   1372 | {9,9988,11689}
    647 | {12,8740}
    811 | {12,8740,9422}
      9 | {31,2264}
      1 | {33,71}
   1136 | {33,71,10663}
   1137 | {33,10663}
    954 | {62,9988}
   1373 | {62,9988,11689}
   1135 | {71,10663}
     28 | {91,3345}
    297 | {95,6545}
    541 | {95,6545,8246}
    617 | {95,6545,8246,8597}
    615 | {95,6545,8597}
    375 | {101,7135}
      7 | {118,2105}
    101 | {118,2105,4445}
    130 | {118,2105,4806}
    316 | {118,2105,6683}
    ...

The obvious issue now is that the result will contain many subsets of arrays considering other rows (I guess this is got to do with the rownum comparator?).

To overcome this I added a second script which will remove subsets, e.g.

DELETE FROM clusters
WHERE ret NOT IN (
    SELECT ret FROM clusters AS a
    WHERE NOT EXISTS (
        SELECT *
        FROM clusters AS b
        WHERE a.ret <@ b.ret 
        AND a.ret <> b.ret
    )      
); 

Ultimately this will leave me with a table without subsets but it will yet have rows with common array members such as

570 | {3,2262,8144,8437}
682 | {3,2262,8144,8913}
...
..

And to "merge" these I figured to introduce a second recursion script similarly to https://stackoverflow.com/questions/46485699/postgresql-group-by-array-elements-in-common but this is absolutely slow. It will take 30 minutes for this small dataset and if I use a dataset with approx 5.000 buildings it will never stop running after 5 hours or more. The script in particular I am using for the array merges is

WITH RECURSIVE
   cte(rownum, ret) AS (  
      SELECT rownum, ret FROM clusters
      UNION
      SELECT cs.rownum, array_merge(c.ret, cs.ret)
      FROM cte AS c
      JOIN clusters AS cs 
      ON c.rownum <> cs.rownum
      AND c.ret && cs.ret
   )
SELECT count(DISTINCT ret) FROM cte; 
-- keeping it simple but here I join with buildings again and group them by cluster and union them...

Where array_merge is a helper function:

CREATE OR REPLACE FUNCTION public.array_merge(arr1 anyarray, arr2 anyarray)
   RETURNS anyarray
   LANGUAGE SQL IMMUTABLE
AS $function$
   SELECT array_agg(DISTINCT elem ORDER BY elem)
   FROM (
      SELECT unnest(arr1) elem 
      UNION
      SELECT unnest(arr2)
   ) s
$function$;

I have the feeling that I am missing something in the first recursion script and possibly the second is not even needed? Can I anyhow simplify this pipeline?

Some further notes:

• the buildings have a geometry index added

• to the generated clusters table I add a gin__int_ops index on the building_ids array column to speed up the && operator

Answer 1

You really want to use the mighty ST_ClusterDBSCAN:

SELECT *,
       ST_ClusterDBSCAN(geom, 0, 1) OVER() AS clst_id
FROM   <your_table>
;

where clst_id is an integer value representing the cluster a row (geometry) belongs to.

With a eps distance of 0, the function effectively clusters by intersection only.