Seeking PostGIS equivalent to GeoPandas union overlay?

Question

What is the PostGIS / PostgreSQL equivalent of the GeoPandas Overlay Union operation?

res_union = geopandas.overlay(df1, df2, how='union')

Say I have two tables, 1 green, 1 red with each with two rows with geometry (polygons):

Then, the objective is to get:

which is a table with 7 rows/ features. See the GeoPandas docs for more info.

I've tried spatial joins and using the ST_Union and ST_Intersection to no avail. I managed to reproduce the equivalent of

res_intersection = geopandas.overlay(df1, df2, how='intersection')

Which is the equivalent of:

using the following SQL command:

SELECT
    table1.letter,
    table2.number,
    ST_Intersection(table1.geom,table2.geom)
FROM 
    table1,
    table2
WHERE
    ST_Intersects(table1.geom,table2.geom)

One solution is to create a table for the intersections and one for the symmetrical differences and than outer join. However I don't know how to create the Symmetrical Difference:

SELECT
    test.hybas_nld.pfaf_id,
    test.gadm_nld.gid_1,
    ST_SymDifference(test.hybas_nld.geom,test.gadm_nld.geom)
FROM 
    test.hybas_nld,
    test.gadm_nld
WHERE
    ST_????(test.hybas_nld.geom,test.gadm_nld.geom)

with the result of:

Possible duplicate with the exception that I would like to keep the polygons where only geometry of table 1 or 2 exists (symmetrical difference).

What would be the simplest and fastest way to get to the final result (Union)?

Answer 1

You can try something like this

select(st_dump(st_collect(st_symdifference(t1.geom,t2.geom),st_intersection(t1.geom,t2.geom)))).geom 
from t1 inner join t2 on st_intersects(t1.geom,t2.geom)

Note: t1 and t2 are names of tables containing polygon geometry and name of geometry column in my case is geom in both tables, it may be different in your case.

Answer 2

Download the PostGIS Addons from this link: https://github.com/pedrogit/postgisaddons

Install by running the postgis_addons.sql file.

Test by running the postgis_addons_test.sql file.

Here is a self contained example of a problem similar to your one:

WITH geomtable AS (
  SELECT 1 id, ST_GeomFromText('POLYGON((0 0, 0 2, 2 2, 2 0, 0 0), (0.2 0.5, 0.2 1.5, 0.8 1.5, 0.8 0.5, 0.2 0.5))') geom
  UNION ALL
  SELECT 2 id, ST_GeomFromText('POLYGON((1 0.2, 1 1, 3 1, 3 0.2, 1 0.2))') geom
  UNION ALL
  SELECT 3 id, ST_GeomFromText('POLYGON((1.5 0.8, 1.5 1.2, 2.5 1.2, 2.5 0.8, 1.5 0.8))') geom
  UNION ALL
  SELECT 4 id, ST_GeomFromText('MULTIPOLYGON(((3 0, 3 2, 5 2, 5 0, 3 0)), ((4 3, 4 4, 5 4, 5 3, 4 3)))') geom
)
SELECT a.id, unnest(ST_SplitAgg(a.geom, b.geom, 0.00001)) geom
FROM geomtable a,
     geomtable b
WHERE ST_Equals(a.geom, b.geom) OR
      ST_Contains(a.geom, b.geom) OR
      ST_Contains(b.geom, a.geom) OR
      ST_Overlaps(a.geom, b.geom)
GROUP BY a.id;

In your query, just replace geomtable with your two tables:

WITH geomall AS (
  SELECT * FROM table1
  UNION ALL
  SELECT * FROM table2
)
SELECT unnest(ST_SplitAgg(a.geom, b.geom, 0.00001)) geom
FROM geomall a,
     geomall b
WHERE ST_Equals(a.geom, b.geom) OR
      ST_Contains(a.geom, b.geom) OR
      ST_Contains(b.geom, a.geom) OR
      ST_Overlaps(a.geom, b.geom)
GROUP BY a.geom;

Should work for thousands of polygons and when there are more than two overlaps.

Answer 3

Union is the equivalent of three operations

1) Intersection
result : 3 polygons
2) Difference(1,2)
result: 1 polygon, 1 multipolygon
3) Difference(2,1)
result: 1 polygon, 1 multipolygon

-- input data
with polys1 AS (
  SELECT 1 df1, ST_GeomFromText('Polygon((0 0, 2 0, 2 2, 0 2, 0 0))') g
  UNION ALL
  SELECT 2, ST_GeomFromText('Polygon((2 2, 4 2, 4 4, 2 4, 2 2))')
  UNION ALL
  SELECT 3, ST_GeomFromText('Polygon((-2 2, -4 2, -4 4, -2 4, -2 2))')
),
polys2 AS (
  SELECT 1 df2, ST_GeomFromText('Polygon((1 1, 3 1, 3 3, 1 3, 1 1))') g
  UNION ALL
  SELECT 2, ST_GeomFromText('Polygon((3 3, 5 3, 5 5, 3 5, 3 3))')
),
-- intersections
intersections AS (
  SELECT df1, df2, ST_INTERSECTION(a.g, b.g) i, a.g AS g1, b.g AS g2 
  FROM polys1 a, polys2 b WHERE ST_INTERSECTS(a.g, b.g)
),
-- per-row union of intersections with this row
diff1 AS (
  SELECT df1, ST_UNION(i) i FROM intersections GROUP BY df1
),
diff2 AS (
  SELECT df2, ST_UNION(i) i FROM intersections GROUP BY df2
),
-- various combinations of intersections
pairs AS (
  SELECT df1, df2, i AS g FROM intersections
  UNION ALL
  SELECT 
    p.df1,
    NULL,
    CASE
      WHEN i IS NULL THEN g 
      ELSE ST_DIFFERENCE(g, i)
    END
  FROM polys1 p LEFT JOIN diff1 d ON p.df1 = d.df1
  UNION ALL
  SELECT
    NULL,
    p.df2,
    CASE
      WHEN i IS NULL THEN g
      ELSE ST_DIFFERENCE(g, i)
    END
  FROM polys2 p LEFT JOIN diff2 d ON p.df2 = d.df2  
)
SELECT * FROM pairs WHERE NOT ST_IsEmpty(g);

Answer 4

I saw this question you posted Why is Union in ArcMap much faster than other approaches? and the query you used for the PostGIS "Union" could be approached a little differently. I based this query off the answer you have above. Maybe this method will be faster for you

with polys1 AS (
  SELECT 1 df1, ST_GeomFromText('Polygon((0 0, 2 0, 2 2, 0 2, 0 0))') geom
  UNION ALL
  SELECT 2, ST_GeomFromText('Polygon((2 2, 4 2, 4 4, 2 4, 2 2))')
  UNION ALL
  SELECT 3, ST_GeomFromText('Polygon((-2 2, -4 2, -4 4, -2 4, -2 2))')
),
polys2 AS (
  SELECT 1 df2, ST_GeomFromText('Polygon((1 1, 3 1, 3 3, 1 3, 1 1))') geom
  UNION ALL
  SELECT 2, ST_GeomFromText('Polygon((3 3, 5 3, 5 5, 3 5, 3 3))')
)

SELECT  df1,  df2, ST_INTERSECTION(a.geom, b.geom) shape 
FROM polys1 a join polys b ON ST_INTERSECTS(a.geom, b.geom)
union
SELECT df1,null df2,coalesce(st_multi(st_collectionextract(st_difference(a.g,(SELECT st_union(b.geom)
                                                                    FROM polys2 b
                                                                    WHERE st_dwithin(a.geom, b.geom, .001) )), 3)), a.geom) shape
FROM polys1 a
union
SELECT null df1,df2,coalesce(st_multi(st_collectionextract(st_difference(b.geom,(SELECT st_union(a.geom)
                                                                    FROM polys1 a
                                                                    WHERE st_dwithin(a.geom, b.geom, .001) )), 3)), b.geom) shape
FROM polys2 b