I've been trying to create a union* of two shapefiles using a Unix-based approach. So far without success.

My two shapefiles can be found here.

The process using Union (analysis) takes approximately 10 minutes in ArcMap on my laptop and creates the desired result. I am looking for a Unix- or cloud-based alternative.

So far I've explored several options:

  1. Geopandas on AWS EC2 instance (large instance). Fails to execute. Stuck after 20 min.
  2. QGIS Doesn't handle invalid geometries that well. Also much slower.
  3. PostGIS on AWS RDS approach:

The query for two similar shapefiles:

CREATE TABLE test.union_v12 AS
-- input data
with polys1 AS (
    y2018m10d01_rh_icep_basins_postgis_v01_v02 t1

polys2 AS (
    y2018m11d14_rh_whymap_to_rds_v01_v01 t1
-- left and right unions
union1 AS (
  SELECT ST_UNION(geom) FROM polys1
union2 AS (
  SELECT ST_UNION(geom) FROM polys2
-- various combinations of intersections
pairs AS (
  SELECT icepbasinid, aqid, ST_INTERSECTION(a.geom, b.geom) geom FROM polys1 a, polys2 b WHERE ST_INTERSECTS(a.geom, b.geom) 
  SELECT icepbasinid, NULL, ST_DIFFERENCE(geom, (SELECT * FROM union2)) geog FROM polys1
  SELECT NULL, aqid, ST_DIFFERENCE(geom, (SELECT * FROM union1)) geog FROM polys2
SELECT * FROM pairs WHERE NOT ST_IsEmpty(geom)

the query works for smaller shapefiles and it's now running for 50 minutes. This approach is not as fast as ArcMap. Also, the query was not so straight forward.

What would be a fast and reliable way to create a "union" (in the ArcMap sense of the word) of two relatively large shapefiles on a Unix platform? A code based approach (CLI / Python) is preferred.

*What ArcMap defines as a "union" is also known as an overlay, geospatial outer join or spatial OR join.

  • 1
    try using st_collect instead of st_union -- that will help performance. other than that this is how I would formulate an ArcGIS Union isn PostGIS – ziggy Nov 15 '18 at 14:12
  • thank you, ST_Collect improves performance but not significantly and not enough to make this approach work. – RutgerH Nov 15 '18 at 14:36
  • 1
    okay well then I would look into shapelys cascaded_union, I have used it once and I believe it is similar to the arcgis union – ziggy Nov 15 '18 at 14:39
  • ok great. I've not used cascade_union before but by reading the docs it seems that it dissolves the shapes. The preferred outcome is similar to how ArcMap works. With separate parts i.e. A\B AUB B\A – RutgerH Nov 15 '18 at 14:46
  • 1
    you could also improve this method. you want to avoid using st_intersection and st_difference. for st_intersection you should do this instead CASE WHEN st_within(c.shape, pp.shape) THEN c.shape ELSE ST_INTERSECTION(c.shape, pp.shape) END AS shape – ziggy Nov 15 '18 at 17:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.