I need to get the Difference of an input layer with ~30k polygons and an overlay layer with ~500k polygons (of which ~200k intersect the input layer).
I am trying to learn PostGIS to do this, as it will (hopefully) be much faster than e.g. QGIS. Newbie at PostGIS so probably making stupid mistakes (sorry).
I have found it easy (and fast) to use ST_Difference across the datasets as set out here. E.g.
SELECT COALESCE(ST_Difference(input.geom,overlay.geom),input.geom) FROM input LEFT JOIN overlay ON ST_Intersects(input.geom,overlay.geom)
However, this returns every pairwise Difference. So rather than returning only those elements of the input layer which are not intersected by any element of the overlay layer, it instead gives multiple copies of the input polygons, stacked on top of each other, each with a different bite out of them corresponding to a single intersecting polygon from the overlay layer.
EDIT: Don't bother to read the rest of my original question, which is all about the bad way I tried to solve this. Skip to the excellent answer from @JGH and @geozelot.
I have tried to solve this using a double-loop approach building on the the second answer here. Essentially:
For each input polygon (first loop): Find the subset of the overlay polygons which intersect that input polygon (using ST_Intersects) (Second, nested loop) Run ST_Difference recursively with each intersecting overlay polygon Replace the input polygon with the result of the recursive ST_Difference
My SQL code attempting to implement this is at the bottom of this post and I know it's broken (I am a complete SQL newbie and have made a stupid mistake I don't understand - EXPLAIN output also given below).
Is this a viable approach (if I can fix the code), or is there a simpler/better way to do this? (It seems like a common enough GIS requirement...?)
My current (broken) code:
create function ST_Layer_Difference(input_geom geometry,overlay_geom geometry) returns geometry as $$ declare i integer; j integer; input_array geometry; output_array geometry; local_overlays_array geometry; temp_geom geometry; output_geom geometry; begin input_array := array(select input_geom); FOR i IN 1..array_length(input_array,1) LOOP temp_geom := input_array[i]; local_overlays_array := array(select ST_Intersects(input_array[i],overlay_geom)); FOR j IN 1..array_length(local_overlays_array,1) LOOP temp_geom := ST_Difference(temp_geom,local_overlays_array[j]); END LOOP; output_array[i] := temp_geom; END LOOP; output_geom := unnest(output_array); return output_geom; end; $$ LANGUAGE plpgsql; create table public."output" as select ST_Layer_Difference(input.geom,overlay.geom) from input, overlay;
And the output of explain (showing I've gone wrong on line 1!):
ERROR: syntax error at or near "function" LINE 1: ...FERS false, SUMMARY false, SETTINGS false) create function S... ^ SQL state: 42601 Character: 132