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:
- Geopandas on AWS EC2 instance (large instance). Fails to execute. Stuck after 20 min.
- QGIS Doesn't handle invalid geometries that well. Also much slower.
- PostGIS on AWS RDS approach:
The query for two similar shapefiles:
CREATE TABLE test.union_v12 AS -- input data with polys1 AS ( SELECT t1.icepbasinid, t1.geom FROM y2018m10d01_rh_icep_basins_postgis_v01_v02 t1 ), polys2 AS ( SELECT t1.aqid, t1.geom FROM 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) UNION ALL SELECT icepbasinid, NULL, ST_DIFFERENCE(geom, (SELECT * FROM union2)) geog FROM polys1 UNION ALL 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.