Using PostGIS to find the overall difference between two (large) polygon datasets

Question

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

Answer 1

Expanding on this answer, we can build a query that will process each geometry of table 1, then, for each of them (via the lateral join), union all intersecting geometries from the 2nd layer and finally keep either the entire original shape or compute the difference

SELECT ST_Multi(COALESCE(
         ST_Difference(a.geom, blade.geom),
         a.geom
       )) AS geom
FROM   table1 AS a
CROSS JOIN LATERAL (
  SELECT ST_Union(b.geom) AS geom
  FROM   table2 AS b
  WHERE  ST_Intersects(a.geom, b.geom) 
) AS blade
;