PostGIS recursive intersection between polygons

Question

I am trying to perform a recursive intersection between all the polygons in a spatial table and get the resulting (many) polygons and information about each intersection for each of them.

Image (not to scale) to explain this:

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

So far I have a query that performs 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) 

Image of the result that I get now for polygons A and A+C :

Here is a test script using squares in images as data. Looking at the area column, it’s clear that some recursive ST_Difference are missing, I just can’t figure out how to do this. Any idea is welcome.

 -- 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 

In this particular example, of course, the window function is not necessary, but this code is a simplified version of my real case, which does a few more things on the side.

I am using PostgreSQL 9.2 and PostGIS 2.0

Answer 1

ST_DIFFRENCE does not have to be recursive, you already have all the polygons, so from each geometry you must subtract the union of other geomes that contain this ret, but are not equal to it. This works, so you should do it whatever you want:

  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