Getting the 'composite' of two polygon layers

Question

I have two table of polygons, and I want to 'flatten' them into one table, and have the result tell me which areas are in one table of polygons, which in the other table of polygons, and which are in both.

I have had a go at doing this in PostGIS, and below present a self-contained example, but is there a better way?

DROP TABLE IF EXISTS a_table;
CREATE TEMP TABLE a_table AS (
    SELECT  'a1' as id, st_geomfromtext('POLYGON((50 50, 50 150, 150 150, 150 50, 50 50))') geom
    UNION ALL
    SELECT  'a2' as id, st_geomfromtext('POLYGON((150 150, 150 200, 200 200, 200 150, 150 150))') geom
    UNION ALL
    SELECT  'a3' as id, st_geomfromtext('POLYGON((10 10, 10 20, 20 20, 20 10, 10 10))') geom
    );
DROP TABLE IF EXISTS b_table;
CREATE TEMP TABLE b_table AS (
    SELECT  'b1' as id, st_geomfromtext('POLYGON((40 40, 40 140, 140 140, 140 40, 40 40))') geom
    UNION ALL
    SELECT  'b2' as id, st_geomfromtext('POLYGON((30 30, 30 40, 40 40, 40 30, 30 30))') geom
    UNION ALL
    SELECT  'b3' as id, st_geomfromtext('POLYGON((155 155, 155 175, 175 175, 175 155, 155 155))') geom
    );
DROP TABLE IF EXISTS intersects;
CREATE TEMP TABLE intersects AS (
    SELECT      a_table.id aid,
                b_table.id bid,
                CASE WHEN   st_within(a_table.geom, b_table.geom) THEN a_table.geom
                ELSE        st_intersection(a_table.geom, b_table.geom)
                END AS      geom
    FROM        a_table
    INNER JOIN  b_table
    ON          (st_intersects(a_table.geom, b_table.geom)
    AND NOT     st_touches(a_table.geom, b_table.geom))
    );
DROP TABLE IF EXISTS non_intersect_diff;
CREATE TEMP TABLE non_intersect_diff AS (
    SELECT      a_table.id aid,
                '' as bid,
                a_table.geom
    FROM        a_table
    WHERE       a_table.id NOT IN (SELECT aid FROM intersects GROUP BY aid)
    );
DROP TABLE IF EXISTS intersect_diff;
CREATE TEMP TABLE intersect_diff AS (
    WITH collect AS (
        SELECT      a_table.id aid,
                    st_union(b_table.geom) geom
        FROM        a_table
        INNER JOIN  b_table
        ON          (st_intersects(a_table.geom, b_table.geom)
        AND NOT     st_touches(a_table.geom, b_table.geom))
        GROUP BY    a_table.id
        )
        SELECT      a_table.id as aid,
                    '' as bid,
                    st_difference(a_table.geom, collect.geom) geom
        FROM        a_table
        INNER JOIN  collect
        ON          a_table.id = collect.aid);
DROP TABLE IF EXISTS non_intersect_diff_reverse;
CREATE TEMP TABLE non_intersect_diff_reverse AS (
    SELECT      '' as aid,
                b_table.id as bid,
                b_table.geom
    FROM        b_table
    WHERE       b_table.id NOT IN (SELECT bid FROM intersects GROUP BY bid)
);
DROP TABLE IF EXISTS intersect_diff_reverse;
CREATE TEMP TABLE intersect_diff_reverse AS (
    WITH collect AS (
        SELECT      b_table.id bid,
                    st_union(a_table.geom) geom
        FROM        b_table
        INNER JOIN  a_table
        ON          (st_intersects(b_table.geom, a_table.geom)
        AND NOT     st_touches(b_table.geom, a_table.geom))
        GROUP BY    b_table.id
        )
        SELECT      '' as aid,
                    b_table.id bid,

                    st_difference(b_table.geom, collect.geom) geom
        FROM        b_table
        INNER JOIN  collect
        ON          b_table.id = collect.bid
);
CREATE TEMP TABLE result AS (
    SELECT * FROM intersects
    UNION ALL
    SELECT * FROM non_intersect_diff
    UNION ALL
    SELECT * FROM intersect_diff
    UNION ALL
    SELECT * FROM non_intersect_diff_reverse
    UNION ALL
    SELECT * FROM intersect_diff_reverse
);
SELECT * FROM result

Answer 1

One way to do this is to create a polygonal coverage from the linework of the inputs, and then spatially join this back to recover the parent attributes. The advantage of using this technique is that it should be more robust in handling polygons which have very small overlaps or coincident linework. It also avoids the need to drop polygon intersection which result in lines.

• Extract linework using ST_Boundary
• Node/dissolve lines using ST_Union
• Polygonize resultants using ST_Polygonize
• Generate an interior point for each resultant using ST_PointOnSurface
• Attach parent attribution by joining on interior points using ST_Contains

WITH poly_a(id, geom) AS (VALUES
    ( 'a1', 'POLYGON ((10 40, 30 40, 30 10, 10 10, 10 40))'::geometry ),
    ( 'a2', 'POLYGON ((70 10, 30 10, 30 90, 70 90, 70 10), (40 40, 60 40, 60 20, 40 20, 40 40), (40 80, 60 80, 60 60, 40 60, 40 80))'::geometry ),
    ( 'a3', 'POLYGON ((40 40, 60 40, 60 20, 40 20, 40 40))'::geometry )
)
,poly_b(id, geom) AS (VALUES
    ( 'b1', 'POLYGON ((90 70, 90 50, 50 50, 50 70, 90 70))'::geometry ),
    ( 'b2', 'POLYGON ((90 30, 50 30, 50 50, 90 50, 90 30))'::geometry ),
    ( 'b2', 'POLYGON ((90 10, 70 10, 70 30, 90 30, 90 10))'::geometry )
)
,lines AS ( 
  SELECT ST_Boundary(geom) AS geom FROM poly_a
  UNION ALL
  SELECT ST_Boundary(geom) AS geom FROM poly_b
)
,noded_lines AS ( SELECT ST_Union(geom) AS geom FROM lines ) 
,resultants AS (  
  SELECT geom, ST_PointOnSurface(geom) AS pip 
    FROM St_Dump(
           ( SELECT ST_Polygonize(geom) AS geom FROM noded_lines ))   
)
SELECT a.id AS ida, b.id AS idb, r.geom
  FROM resultants r
  LEFT JOIN poly_a a ON ST_Contains(a.geom, r.pip) 
  LEFT JOIN poly_b b ON ST_Contains(b.geom, r.pip)
  WHERE a.id IS NOT NULL OR b.id IS NOT NULL;