Applying ST_ClusterDBSCAN repeatedly until cluster sizes small enough

Question

I am looking for a way to spatially clusters thousands of datapoints (potentially millions) based on distance, such that each cluster contains less than 5000 points.

This is a similar question to Problems with ST_ClusterDBSCAN cluster sizes . I would like to build upon the provided answer by using WITH RECURSIVE to automatically continue splitting clusters until they are all bellow a size.

This is the query I came up with (not complete):

WITH RECURSIVE clusterize(cid, csize, autopoi_ids, eps) AS (
    SELECT cid, csize, unnest(poi_ids) as poi_id, eps
    FROM (
        SELECT cid, count(*) as csize, array_agg(id) as poi_ids, 0.05 as eps
        FROM (
            SELECT id, ST_ClusterDBSCAN(geometry, eps := 0.05, minpoints := 3) over () AS cid
            FROM stats_autopoistat
        ) clusters
        GROUP BY cid
    ) q

    UNION ALL

    SELECT cid, csize, unnest(poi_ids) as poi_id, eps
    FROM (
        SELECT cid, count(*) as csize, array_agg(id) as poi_ids, ( SELECT eps/2.0 FROM clusterize LIMIT 1 )/2.0 as eps
        FROM (
            SELECT id, (SELECT max(cid) FROM clusterize) + ST_ClusterDBSCAN(geometry, eps := ( SELECT eps/2.0 FROM clusterize LIMIT 1), minpoints := 0) over () AS cid
            FROM clusterize
            WHERE csize > 5000
        ) clusters
        GROUP BY cid
    ) q
)
SELECT *
-- here filter out non-max cids for each poi_id
FROM clusterize limit 1000

However, it seems I am unable to refer to the recursive CTE inside a subquery, as Postgres complains with:

ERROR:  recursive reference to query "clusterize" must not appear within a subquery
LINE 15: ..., array_agg(id) as poi_ids, ( SELECT eps/2.0 FROM clusterize...

I would like to know if this can even be come with WITH RECURSIVE given the limitations I encountered above.

The reason I want to accomplish this within Postgres and not Python is that the number of points to cluster will continue increasing. The table already has about 1 million rows, and I would like to avoid loading all this data into Python.

Answer 1

No time for more improving or testing, but: for a single, more generic recursive term, and possibly better performance, try

WITH RECURSIVE
    params AS (                      -- convenience variables for testing parameters
        SELECT  10 AS max_size,      -- max. cluster size
                1 AS min_points,     -- 'min_points' parameter
                1 AS eps,            -- 'eps' distance parameter
                0.1 AS fraction      -- decreasing fraction of/to 'eps' parameter
    ),
    clst AS (
        SELECT  ARRAY[a._clst_id] AS _clst_ids,
                1 - (1 * (SELECT fraction FROM params)) AS _eps,
                ST_Collect(a.geom) AS geom
        FROM    (
            SELECT  id,
                    ST_SetSRID(ST_MakePointM(ST_X(geom), ST_Y(geom), id), 4326) AS geom,
                    ST_ClusterDBSCAN(geom, (SELECT eps FROM params), (SELECT min_points FROM params)) OVER() AS _clst_id
            FROM    <pts>
        ) AS a
        GROUP BY
                _clst_id
        UNION ALL
        SELECT  CASE WHEN ST_NumGeometries(b.geom) > (SELECT max_size FROM params)
                    THEN a._clst_ids || b._clst_id
                    ELSE NULL
                END AS _clst_ids,
                a._eps - (a._eps * (SELECT fraction FROM params)) AS _eps,
                b.geom AS geom
        FROM    clst AS a
        CROSS JOIN LATERAL (
            SELECT  ST_Collect(c.geom) AS geom,
                    c._clst_id
            FROM    (
                SELECT  dmp.geom,
                        ST_ClusterDBSCAN(dmp.geom, a._eps, (SELECT min_points FROM params)) OVER() AS _clst_id
                FROM    LATERAL ST_DumpPoints(a.geom) AS dmp
            ) c
            GROUP BY
                    c._clst_id
        ) b
        WHERE   ST_NumGeometries(a.geom) > (SELECT min_size FROM params)
    )
SELECT  ST_M(geom)::INT AS id,
        ST_Force2d(geom) AS geom,
FROM    (
    SELECT  ROW_NUMBER() OVER() AS clst_id,
            (ST_DumpPoints(geom)).geom
    FROM    clst
    WHERE   _clst_ids IS NULL
) q
;

This approach ST_Collects points based on their _clst_id and recursively processes those (each row in clst) with ST_NumGeometries > max_size using a LATERAL JOIN. If a cluster has reached max_size, it get's NULL as _clst_ids to mark it as a finished cluster.

I used params.fraction = 0.1 to decrease the eps distance, which is pretty intense; smaller values will yield more precise results, but increase the execution time (probably) exponentially.

Since geometry aggregation makes it a pain in the to keep attributes along the way, and a join on geometric equality with very large tables to retrieve the original attributes is costly, I write the id of each point into the M coordinate of the points and extract them later. This only works with numeric values.

If you are interested in MultiPoint geometries per cluster, just remove those parts and the dump in the final query.

---

It would probably be a better idea to write a function for this; I couldn't say if a DO ... WHILE loop would perform better than the WITH RECURSIVE implementation, but you could work with attributes a lot better (and probably more performant, especially if you are interested in other original attributes than the id).

Answer 2

I have been able to work around the limitation by "pre-computing" the values for eps and inferring reasonable values of the other subqueries which where previously referring to the recursive CTE.

Note that the new solution may create clusters larger than wanted (5000 in the query below) if you run-out of "pre-computed" values. This helps ensure that the query does not infinitely loop if more than 5000 points have the exact same coordinates (in which case the clustering algorithm would always assign them the same cid).

This is the working solution:

WITH RECURSIVE
    row_count AS (
        SELECT count(*) as count FROM stats_autopoistat
    ), magic_constants AS (
        select *
        FROM ROWS FROM (
            generate_series(1,4,1),
            generate_series(0, 3 * (SELECT count FROM row_count), (SELECT count FROM row_count)),
            unnest(ARRAY[0.05, 0.025, 0.0125, 0.00625])
        ) AS t(iter, starting_cid, eps)
    ),
    clusterize(cid, iter, csize, poi_ids, eps) AS (
        SELECT cid::bigint, 1 as iter, count(*) as csize, array_agg(id) as poi_ids, 0.05 as eps
        FROM (
            SELECT id, ST_ClusterDBSCAN(geometry, eps := 0.05, minpoints := 3) over () AS cid
            FROM stats_autopoistat
        ) clusters
        GROUP BY cid

        UNION ALL

        SELECT cid, (min(iter) + 1)::integer as iter, count(*) as csize, array_agg(id) as poi_ids, min(eps) as eps
        FROM (
            SELECT s.id, c.iter, eps.eps, eps.starting_cid + ST_ClusterDBSCAN(geometry, eps := eps.eps, minpoints := 1) over () AS cid
            FROM clusterize c
            JOIN magic_constants eps ON (c.iter = eps.iter)
            LEFT JOIN stats_autopoistat s ON (s.id = ANY(c.poi_ids))
            WHERE csize > 5000
        ) clusters
        GROUP BY cid
    ), poi_cluster AS (
        SELECT DISTINCT ON (poi_id) poi_id, cid, csize
        FROM (
            SELECT cid, iter, csize, unnest(poi_ids) as poi_id
            FROM clusterize
        ) q
        ORDER BY poi_id, iter DESC
    )
SELECT cid, ST_ConcaveHull(ST_Collect(geometry), 0.99) as geometry
FROM poi_cluster c
JOIN stats_autopoistat s ON (c.poi_id = s.id)
WHERE cid IS NOT NULL AND csize >= 3
GROUP BY cid