3

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

enter image description here

Then, the objective is to get:

enter image description here

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:

enter image description here

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: enter image description here

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

2

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.

  • Almost. One comment: The preferred result is that none of the geometries overlap i.e. if you look at the image with the seven colors, the green geometry is a multipolygon with two rectangles. If you query the result layer at a random location, you only find 1 layer i.e. not stacking / overlap. – RutgerH Nov 13 '18 at 11:44
  • can you explain why you use ST_dump and ST_collect ? – RutgerH Nov 13 '18 at 11:56
  • You can check there will be no overlaps, ST_collect is same as st_union only difference is that it does not merge boundaries but it returns geometry collection so to get polygons from geometry collection st_dump is used – Asad Abbas Nov 13 '18 at 12:14
1

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.

0

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);
  • this is how I would do this, although I would do it in one query instead of creating all these tables – ziggy Nov 13 '18 at 18:52
  • @ziggy this has come up a few times in recent weeks, so a solution you're describing would be helpful! – DPSSpatial Nov 13 '18 at 21:55
0

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

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.