2

EDITED QUESTION: PROPER SYNTAX FOR ST_DIFFERENCE BETWEEN TABLES?

I have two huge polygon datasets from shapefiles, call them 'baselevel' coverage and 'deluxe' coverage. They each contain over 2.4 million polygons, with many polygons describing complex coastal features.

The baselevel dataset has complete coverage over the area of interest but less detailed attributes, while the deluxe dataset has gaps and less extensive spatial coverage but more detailed attributes.

I want to create a third dataset which uses all the deluxe polygons, and uses the baselevel coverage to fill the gaps and extend the coverage. I've been approaching this as two steps: produce a 'fill' dataset which is just the base-level coverage minus all the deluxe polygons, then union the 'baselevel-fill' dataset with the 'deluxe' dataset.

I can do this fairly efficiently in grass via v.overlay (create the fill layer with operator=NOT, then union with deluxe using operator=AND), after dicing the shapefiles into much smaller tiles.

Final Update

I'm pasting my final syntax below, though I'd still welcome any optimization suggestions ;-) based on Paul Ramseys answer, as well as some other stackexchange answers and postgis tutorials

GRASS is pretty excellent, and could handle my polygons where a lot of other python, ogr, and qgis approaches. However the fastest I could complete a similar operation on my (very small) test dataset was 154 seconds. Code below completed in under 28 seconds :-)

CREATE TABLE deluxe_plus_baselvl AS
WITH deluxe_cutter AS (
 SELECT ST_Union(d.geom) AS geom, b.gid
 FROM baselevel b JOIN deluxe d ON ST_Intersects(b.geom, d.geom)
 GROUP BY b.gid
 ),
baselevel_cut AS (
 SELECT ST_Difference(b.geom, d.geom) AS geom, b.gid
 FROM baselevel b JOIN deluxe_cutter d ON b.gid = d.gid
 )
SELECT 'baselevel' AS type, b.geom, b.gid
FROM baselevel_cut b
UNION ALL
SELECT 'deluxe' AS type, d.geom, d.gid
FROM deluxe d
UNION ALL
SELECT 'baselevel' AS type, b.geom, b.gid
FROM baselevel b LEFT JOIN deluxe d ON
ST_Intersects(b.geom,d.geom)
where d.gid is null;

Thank you for your help!

1

In order to get the effect you want (the baselevel with all the areas in deluxe removed) you're going to need to first ST_Union deluxe into a single polygon. Then ST_Difference that polygon from each of the polygons in the baselevel that it intersects. It'll be slow, it'll be painful, how slow and painful depends rather on how many vertices there are in the deluxe data set.

Another possibility would be a more piecewise approach, but will still be painful, though perhaps less so in aggregate.

  • For each baselevel polygon get all the deluxe that intersect it
  • Union these deluxe polygons together, now you have a "thing to subtract" for each of your baselevel polygons
  • Do the subtraction, one baselevel polygon at a time
  • Finally just SQL union the original deluxe set with the now differenced baselevels.

Something like this:

WITH deluxe_diffs AS (
   SELECT ST_Union(d.geom) AS geom, b.id
   FROM baselevel b JOIN deluxe d ON ST_Intersects(b.geom, d.geom)
   GROUP BY b.id
)
baselevel_cut AS (
   SELECT ST_Difference(b.geom, d.geom) AS geom, b.id
   FROM baselevel b JOIN deluxe_diffs d ON b.id = d.id
)
SELECT 'baselevel' AS type, b.geom, b.id 
FROM baselevel_cut b
UNION ALL
SELECT 'deluxe' AS type, d.geom, d.id
FROM deluxe d;

Don't quote me on the SQL syntax, it's off top of my head.

  • I see now that this approach would not include baselevel polygons that had no intersections with deluxe ones. I leave including those extra polygons as an exercise to the reader. – Paul Ramsey Sep 13 '13 at 23:20
  • Hmmm. I'll set up some test regions and try out all three approaches. Thank you! – WileyB Sep 14 '13 at 20:56
  • When you create baselevel_cut, are you joining two unrelated datasets on their unique key? baselevel_cut AS...ON b.id = d.id – WileyB Oct 10 '13 at 21:47
  • Each entry in deluxe_diffs is all the deluxe polygons that intersect a particular baselevel polygon, and each is tagged with a baselevel id, so the join can happen in the next subquery on that baselevel id. – Paul Ramsey Oct 11 '13 at 18:47

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.