3

Summary

I have a PostgreSQL table of some 80M polygons, and I am attempting to trace paths of connectivity from known starting points within areas of interest. The data is three-dimensional-ish, organized into 2D layers. Polygons may intersect with other polygons on their own layers, as well as "adjacent" layers, with adjacency defined using a many-to-many table.

This is an existing system, and I'm trying to improve its performance by shifting the heavy lifting of intersection testing over to Postgres. I've managed to knock a three-hour job down to 45 minutes using the approach I'm about to describe, and the results are correct, but I wonder if it could be done even faster if I actually knew what I was doing with postgres.

I present this test case of tracing one known large net from a single polygon, with the tracing query pulled out of a function which manages a work queue, and other boring and probably irrelevant bits of overhead. This example query runs in ~160 seconds.

Obvious Considerations

  • I am using a spatial index on the polygon data.
  • I've seen another answer suggesting ST_ClusterDBSCAN for something like this, but it seems that function works entirely in-memory, and I think my dataset is too large for that to work here (unless I'm mistaken about how it works).
  • Maybe my entire approach is foolish? I'm completely open to suggestions for a better way to come at this problem.

Table Definitions

-- Polygon table
CREATE TABLE polygon(
    ID serial8 primary key,
    layer_id BIGINT NOT NULL REFERENCES layerrow ON DELETE CASCADE,
    geom GEOMETRY(polygon,0) NOT NULL,
    indices BYTEA NOT NULL);
CREATE INDEX ON polygon(layer_id, id);
CREATE INDEX ON polygon USING GIST(layer_id, geom);

-- Layer connection table
CREATE TEMPORARY TABLE tmp_layer_connections(
    id_from BIGINT NOT NULL,
    id_to BIGINT NOT NULL
);
--(Omitted query to load tmp_layer_connections with 30-ish rows)
CREATE INDEX tmp_layer_connections_id_from_idx ON tmp_layer_connections(id_from);

-- Trace results table
CREATE TEMPORARY TABLE tmp_trace_nets(
    polygon_id BIGINT PRIMARY KEY,
    net_id BIGINT NOT NULL
);
CREATE INDEX tmp_trace_nets_net_id_idx ON tmp_trace_nets(net_id);

Test Case Query

This is the core of the trace procedure, annotated a bit. I'm trying to perform indexed spatial joins in a WITH clause, then use a recursive query to trace outwards from a starting point. The stuff surrounding this query all runs quite quickly, and just runs this query on various known starting points, incrementing the net ID number on each cycle. The hardcoded values establish the test case instead of working for the general case.

WITH RECURSIVE trace_map AS (
    -- Build a trace map to support recursive intersection tests on the polygon table
    SELECT p1.id AS poly_from, p2.id AS poly_to FROM polygon AS p1
    INNER JOIN polygon AS p2 -- Join polygon against itself
        ON ST_INTERSECTS(p1.geom, p2.geom) -- Find intersecting polygons
        AND ST_INTERSECTS(p2.geom, (SELECT n.bounds_nm FROM netlist AS n WHERE n.id = 49)) -- Examine only area of interest
        AND p1.id <> p2.id -- Don't return the source polygon
    INNER JOIN tmp_layer_connections AS tc -- Join layer connections table to intersect only touching layers
        ON p1.layer_id = tc.id_from
        AND p2.layer_id = tc.id_to
), trace_path(trace_id) AS (
    SELECT 229346947::BIGINT -- Start tracing from a polygon in a known large net
    UNION
        -- Recursively union in distinct intersecting polygons from allowed layers, using our trace map
        SELECT trace_map.poly_to
        FROM trace_map, trace_path
        WHERE trace_map.poly_from = trace_path.trace_id
)
-- Dump the discovered polygons as belonging to net ID 1
INSERT INTO tmp_trace_nets
SELECT tp.trace_id, 1
FROM trace_path AS tp;

Explain Analyze

Here's the raw EXPLAIN ANALYZE output:

"Insert on tmp_trace_nets  (cost=7123.91..7124.13 rows=11 width=16) (actual time=160072.642..160072.642 rows=0 loops=1)"
"  CTE trace_path"
"    ->  Recursive Union  (cost=0.00..7123.91 rows=11 width=8) (actual time=0.005..159948.349 rows=40136 loops=1)"
"          ->  Result  (cost=0.00..0.01 rows=1 width=8) (actual time=0.001..0.001 rows=1 loops=1)"
"          ->  Hash Join  (cost=660.73..712.37 rows=1 width=8) (actual time=203.168..203.183 rows=115 loops=787)"
"                Hash Cond: ((tc.id_from = p1.layer_id) AND (tc.id_to = p2.layer_id))"
"                InitPlan 1 (returns $1)"
"                  ->  Index Scan using netlist_pkey on netlist n  (cost=0.28..8.29 rows=1 width=120) (actual time=0.069..0.071 rows=1 loops=1)"
"                        Index Cond: (id = 49)"
"                ->  Seq Scan on tmp_layer_connections tc  (cost=0.00..28.50 rows=1850 width=16) (actual time=0.002..0.003 rows=19 loops=787)"
"                ->  Hash  (cost=652.35..652.35 rows=6 width=24) (actual time=203.161..203.161 rows=1347 loops=787)"
"                      Buckets: 1024  Batches: 1  Memory Usage: 9kB"
"                      ->  Nested Loop  (cost=0.85..652.35 rows=6 width=24) (actual time=0.426..202.435 rows=1347 loops=787)"
"                            ->  Nested Loop  (cost=0.44..84.75 rows=10 width=376) (actual time=0.015..1.992 rows=51 loops=787)"
"                                  ->  WorkTable Scan on trace_path  (cost=0.00..0.20 rows=10 width=8) (actual time=0.000..0.022 rows=51 loops=787)"
"                                  ->  Index Scan using polygon_pkey on polygon p1  (cost=0.44..8.46 rows=1 width=376) (actual time=0.037..0.037 rows=1 loops=40136)"
"                                        Index Cond: (id = trace_path.trace_id)"
"                            ->  Index Scan using polygon_layer_id_geom_idx on polygon p2  (cost=0.42..56.75 rows=1 width=376) (actual time=0.476..3.923 rows=26 loops=40136)"
"                                  Index Cond: ((geom && p1.geom) AND (geom && $1))"
"                                  Filter: ((p1.id <> id) AND st_intersects(geom, $1) AND st_intersects(p1.geom, geom))"
"                                  Rows Removed by Filter: 10"
"  ->  CTE Scan on trace_path tp  (cost=0.00..0.22 rows=11 width=16) (actual time=0.008..159957.471 rows=40136 loops=1)"
"Planning Time: 3.006 ms"
"Execution Time: 160074.990 ms"

Here's a link to a more readable version of that EXPLAIN ANALYZE: https://explain.depesz.com/s/Gdlm

Here are screenshots from that link for relevant areas:

Parsed EXPLAIN ANALYZE results

EXPLAIN ANALYZE stats

Commentary

99.3% of the query time is spent on index scans of the polygon table for the ST_INTERSECTS() calls, which makes intuitive sense to me, as I'd expect the majority of the work for this approach to be repeated intersection tests as it walks the path of polygons.

I don't see anything obviously horrible in the EXPLAIN ANALYZE, but as I mentioned, I only barely know what I'm doing, here.

If anyone can suggest a way to make this approach faster, or suggest a better approach to the problem of net tracing, I would be extremely grateful.

Second Approach

On receiving suggestions to give ST_ClusterDBSCAN another try, I've found that I was wrong. I can cluster across large datasets without having to load it all into memory. I've been trying to rewrite my above approach using ST_ClusterDBSCAN instead of my recursive intersection approach.

My new problem is that my data is pseudo-3D (divided into layers, only some allowed to connect) and the clustering can only handle one layer at a time. I thought to partition the polygons by the their layer ID while clustering, then find the interconnection points between the layers, then use those interconnecting polygons to join the clusters across the layers. However, I've hit a brick wall right at the end, and my SQL-fu doesn't seem strong enough to get across the finish line. I feel like I have all the data available for me, but I can't quite manipulate it into the form I need it.

Here's the query so far (the hardcoded values here just establish a test case):

WITH clusters AS (
    SELECT polygon.id AS polygon_id, polygon.layer_id, ST_CLUSTERDBSCAN(geom,0,0) OVER (PARTITION BY layer_id) AS cluster_id
    FROM polygon
    WHERE layer_id = ANY((SELECT DISTINCT id_from FROM tmp_layer_connections))
    AND ST_INTERSECTS(geom, ((SELECT n.bounds_nm FROM netlist AS n WHERE n.id = 172)))
), interconnections AS (
    SELECT p1.id AS from_id, p2.id AS to_id
    FROM polygon AS p1
    INNER JOIN polygon AS p2 ON ST_INTERSECTS(p1.geom, p2.geom)
        AND ST_INTERSECTS(p2.geom, (SELECT n.bounds_nm FROM netlist AS n WHERE n.id = 172))
        AND p1.id <> p2.id
    INNER JOIN tmp_layer_connections AS lc
        ON p1.layer_id = lc.id_from
        AND p2.layer_id = lc.id_to
        AND lc.id_from <> lc.id_to
), segment_links AS (
    SELECT DISTINCT cl_from.layer_id AS from_layer_id, cl_from.cluster_id AS from_cluster_id, cl_to.layer_id AS to_layer_id, cl_to.cluster_id AS to_cluster_id
    FROM interconnections AS ic
    INNER JOIN clusters AS cl_from on cl_from.polygon_id = ic.from_id
    INNER JOIN clusters AS cl_to on cl_to.polygon_id = ic.to_id
)
-- ?????

Here's my thinking on the WITH progression:

  1. "clusters" runs ST_CLUSTERDBSCAN on all polygons in the area of interest, segmenting by layer_id. It results in a beautiful data set where each polygon in the AOI gets clustered with intersecting polygons on the same layer, each cluster identified by cluster_id. I went with minimum cluster size 0 so that lone polygons still get cluster_ids, since they might touch another layer, and are still relevant. Note that the cluster_id count restarts at 0 for each layer, so what I'm calling a "segment" of polygons is a combination of layer_id and cluster_id. Each polygon now belongs to a segment.
  2. "interconnections" creates a many-to-many relationship, with each row indicating two polygons on different layers which touch, according to the AOI and the layer interconnection rules. My hope us to use these connecting polygons to unite the segments and ultimately answer the question "what are the sets of polygon IDs which intersect through the 3D structure?"
  3. This is where it starts to fall apart. "segment_links" is my attempt to massage the data into its final form. It takes each row from the interconnections table and uses the polygon membership in the segments to join in the segments on each side of the connection. I ultimately have a collection of rows indicating segment connectivity. For example, the row {0,9,1,3} indicates "layer-0,cluster-9 connects to layer-1,cluster3." This is still many to many, so there may be many entries for 0,9 to touch clusters in other layers. Reciprocal connections are present, so there's also a {1,3,0,9} row, but it's easy to be rid of those with a modification to "interconnections," if necessary.

However, this is where I'm stuck. I don't know how take these connections and query out the "sets" of connections. To get something like "Net 1: {0,9},{1,3},{0,8},{2,300}" and then similar sets for the other connected segments. If I could do that, then it would be trivial to use "clusters" to turn the collections of segments into collections of connected polygon IDs, which is ultimately what I need.

Is there anything I can do to pull these sets out of segment_links? I tried ending with a recursive query primed with the first row of segment_links to follow the from/to relationships, but that only lets me assemble the set of segments which includes the first row of the segment_links table. I can't figure out how to extract all the individual sets represented.

A more general way to pose this question might be this. Given a table with columns n1 and n2 which looks like this:

n1  n2
1   2
2   1
1   3
3   1
2   4
4   2
10  11
11  10
11  12
12  11

Is there any way to query this data to extract the sets {1,2,3,4},{10,11,12}?

If that's a dead end, is there instead a better way I could organize this stack of WITH clauses?

14
  • 1
    Interesting question, although slightly confused by your end goal here. Is it to group touching polygons into bigger polygons?
    – ziggy
    Apr 5 at 18:00
  • 1
    I am wondering if ST_Subdivide on the polygons would speed this up. Apr 5 at 20:07
  • 2
    the tracing is throwing me off a bit, typically that type of analysis is run with linestring. you can give st_clusterdbscan a try, I've used that for heavy cluster grouping. not 80 million records but a few million for sure
    – ziggy
    Apr 5 at 22:31
  • 1
    I would create a new table with the subdivided polygons and carry everything over. Apr 6 at 4:36
  • 1
    Thinking about this again, I think Dbscan setting eps to 0 ( gis.stackexchange.com/questions/343514/… ) which will give you the connecting polygons is the way to go and you won't require recursion. Your worry regarding memory can be lifted by possibly breaking it down to disjoiint regions or using a beefier machine? Apr 6 at 4:41

1 Answer 1

2

Regarding my example table "test" with columns n1 and n2, I did come up with a query which can get what I was looking for:

SELECT DISTINCT (
    WITH RECURSIVE collector(n) AS (
        SELECT test.n1
    UNION
        SELECT test.n2
        FROM test, collector
        WHERE collector.n = test.n1
    )
    SELECT ARRAY_AGG(collector.n ORDER BY collector.n) AS n_groups FROM collector
) FROM test

For each row of the "test" table, this recursive query walks the n1->n2 relationships reachable from that starting n1, sort of like searching through a graph. The results are thrown into a sorted array, so that rows with duplicate nets will be trimmed down to a single instance each. When run on my "test" table, this query yields:

n_groups
{1,2,3,4}
{10,11,12}

Which is great. However, the fact that it has to assemble and sort an array for every row of the test table makes it lousy for efficiency. I adapted this approach to the test case query I described in the original question, and got this:

WITH clusters AS (
    SELECT polygon.id AS polygon_id, polygon.layer_id, ST_CLUSTERDBSCAN(geom,0,0) OVER (PARTITION BY layer_id) AS cluster_id
    FROM polygon
    WHERE layer_id = ANY((SELECT DISTINCT id_from FROM tmp_layer_connections))
    AND ST_INTERSECTS(geom, ((SELECT n.bounds_nm FROM netlist AS n WHERE n.id = 146)))
), interconnections AS (
    SELECT p1.id AS from_id, p2.id AS to_id
    FROM polygon AS p1
    INNER JOIN polygon AS p2
    ON ST_INTERSECTS(p1.geom, p2.geom)
        AND ST_INTERSECTS(p1.geom, (SELECT n.bounds_nm FROM netlist AS n WHERE n.id = 146))
        AND ST_INTERSECTS(p2.geom, (SELECT n.bounds_nm FROM netlist AS n WHERE n.id = 146))
        AND p1.id <> p2.id
    INNER JOIN tmp_layer_connections AS lc
        ON p1.layer_id = lc.id_from
        AND p2.layer_id = lc.id_to
        AND lc.id_from <> lc.id_to
), segment_links AS (
    SELECT DISTINCT
        cl.layer_id AS from_layer_id,
        cl.cluster_id AS from_cluster_id,
        cl_to.layer_id AS to_layer_id,
        cl_to.cluster_id AS to_cluster_id
    FROM clusters AS cl
    LEFT JOIN interconnections AS ic
        ON cl.polygon_id = ic.from_id
    LEFT JOIN clusters AS cl_to
        ON ic.to_id = cl_to.polygon_id
), nets AS (    
    SELECT DISTINCT(
        WITH RECURSIVE results(layer_id, cluster_id) AS (
            SELECT sl.from_layer_id, sl.from_cluster_id
        UNION
            SELECT sl.to_layer_id, sl.to_cluster_id
            FROM results AS r, segment_links AS sl
            WHERE r.layer_id = sl.from_layer_id
            AND r.cluster_id = sl.from_cluster_id
            AND sl.to_layer_id IS NOT NULL
        )
        SELECT ARRAY_AGG(array[results.layer_id,results.cluster_id] ORDER BY results.layer_id, results.cluster_id) AS by_segment FROM results
    ) FROM segment_links as sl
), nets_exploded AS (
    SELECT
        ROW_NUMBER() OVER (ORDER BY nets.by_segment) AS net_number,
        REDUCE_DIM(nets.by_segment) AS segment
    FROM nets
)
SELECT net_number, cl.polygon_id
FROM nets_exploded AS ne
INNER JOIN clusters As cl ON ne.segment[1] = cl.layer_id AND ne.segment[2] = cl.cluster_id

This query produces the correct answer, but scales extremely poorly. EXPLAIN ANALYZE shows that it spends the majority of its time performing the many redundant sorts. Worse, with even a small set of polygons, it exceeds PostgreSQL's working memory limit, and performs the sort on-disk.

My question has been answered, but I think this approach is dead in the water without a better way to find the sets of linked clusters. Some experimentation suggests that the "clusters" and "interconnections" intermediate tables are created very quickly, so it does seem like ST_ClusterDBSCAN is the correct answer, but my efforts to massage the data into my preferred form run horribly slowly database-side.

I will next work on just returning the raw results of the "clusters" and "interconnections" queries, and manipulating the data with a more suitable language. I have high hopes for that approach, and I think this question is now done. Thank you very much to everyone who offered ideas and commentary.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.