PostGIS recursive intersection between polygons

Question

I'm trying to perform a recursive intersection between all the Polygons in a spatial Table, and obtain the resulting (multi)pololygons and the information about every intersection for each of them.

An image (not really in scale) to explain it: Example

Let's say there are A, B, C squares in a table. I'd like to have A, B, C, A+B, A+C, B+C, A+B+C polygons in output, and I need to know that A+B is the intersection of A and B and so on.

So far I have a query which performs the intersections, but it does not "cut off" the intersected part of the original polygons. For example:

Polygon A should be      A - (A+B) - (A+C) - (A+B+C)
Polygon A+C should be    A+C - (A+B+C)

An image of the result I get now for the A and A+C polygons:

Current WRONG result

Here is a test script, using the squares in the images as data. Looking at the area column, it is clear some recursive ST_Difference is missing, I just can't figure out how. Any idea is welcomed.

-- Create a test table
CREATE TABLE test (
    name text PRIMARY KEY,
    geom geometry(POLYGON)
);

-- Insert test data
INSERT INTO test (name, geom) VALUES 
    ('A', ST_GeomFromText('POLYGON((1 2, 1 6, 5 6, 5 2, 1 2))')),
    ('B', ST_GeomFromText('POLYGON((0 0, 0 4, 4 4, 4 0, 0 0))')),
    ('C', ST_GeomFromText('POLYGON((2 0, 2 4, 6 4, 6 0, 2 0))'));


-- Query    
WITH RECURSIVE 
source (rownum, geom, ret) AS (
    SELECT row_number() OVER (ORDER BY name ASC), ST_Multi(geom), ARRAY[name] FROM test 
),
r (rownum, geom, ret, incroci) AS (
    SELECT rownum, geom, ret, 0 FROM source 
    UNION ALL
    SELECT s.rownum, ST_CollectionExtract(ST_Intersection(s.geom, r.geom), 3), (r.ret || s.ret), (r.incroci + 1) 
        FROM source AS s INNER JOIN r ON s.rownum > r.rownum AND ST_Intersects(s.geom, r.geom) AND ST_Area(ST_Intersection(s.geom, r.geom)) > 0.5
),
result (geom, ret) AS (
    SELECT ST_Union(geom) AS geom, ret FROM r GROUP BY ret
)
SELECT geom, ST_Area(geom) AS area, ret FROM result ORDER BY ret

The window function isn't strictly necessary in this particular example of course, but this code is a simplified version of my real case, which does a few more things on the side.

I'm using PostgreSQL 9.2 and PostGIS 2.0

Jakub Kania Jakub Kania · Accepted Answer · 2013-03-25T23:32:01

ST_DIFFRENCE doesn't have to be recursive, you already have all the polygons so from every geom you have to substract the union of other geoms which contain that ret but ain't equal to it. This works so you should do it kinda like that:

    WITH RECURSIVE 
source (rownum, geom, ret) AS (
    SELECT row_number() OVER (ORDER BY name ASC), ST_Multi(geom), ARRAY[name] FROM test 
),
r (rownum, geom, ret, incroci) AS (
    SELECT rownum, geom, ret, 0 FROM source 
    UNION ALL
    SELECT s.rownum, ST_CollectionExtract(ST_Intersection(s.geom, r.geom), 3), (r.ret || s.ret), (r.incroci + 1) 
        FROM source AS s INNER JOIN r ON s.rownum > r.rownum AND ST_Intersects(s.geom, r.geom) AND ST_Area(ST_Intersection(s.geom, r.geom)) > 0.5
),
result (geom, ret) AS (
    SELECT ST_Difference(ST_Union(r.geom),q.geom) AS geom, r.ret FROM r JOIN (SELECT r.ret,ST_UNION(COALESCE(r2.geom,ST_GeomFromText('POLYGON EMPTY'))) as geom FROM r LEFT JOIN r AS r2 ON r.ret<@r2.ret AND r.ret!=r2.ret GROUP BY r.ret) AS q on r.ret=q.ret GROUP BY r.ret,q.geom
)
SELECT geom, ST_Area(geom) AS area, ret FROM result ORDER BY ret

PostGIS recursive intersection between polygons

1 Answers