Find what polygons are not fully covered by union of polygons from another layer in PostGIS

Question

I have two tables in a PostGIS database

>>> select postgis_version();
"2.5 USE_GEOS=1 USE_PROJ=1 USE_STATS=1"

One represents the country borders and another one is administrative boundaries. Both are of global coverage covering the whole world. Both tables have spatial indices on the geometry columns.

My goal is to find what countries are not fully covered by the administrative boundaries and, if feasible to achieve within a reasonable time, the location/geometry of part of country geometry where it is not covered by the administrative boundaries.

To give you an example.

• Two countries (marked as 1 (consists of three polygons) and 2 (consists of two polygons), dark green polygons (two rows in table in this case).

• Three admin boundary polygons (marked as A, B, C) (can be multipolygon/polygon), transparent blue polygons (three rows in table in this case).

• Country 1 is fully covered by the admin boundary A and B. Country 2 is not fully covered by admin boundaries (C), so I would like to get this reported. If possible, the hatched area of country 2 could be reported as well.

This sounded like a job for ST_Difference function however my problem is that it is just too slow running on the global dataset.

Running this SQL query on Singapore takes 4 seconds. This is because I have just a few admin boundaries located around the country which I've extracted beforehand and I use only single country geometry. I get the results I expect.

drop table sg_missing_;
WITH SG AS (SELECT iso_cc, geom FROM countries_agg WHERE iso_cc = 'SG'),
    AD AS (SELECT ST_Union(geom) AS geom FROM admins_sample) 
SELECT ST_Difference(c.geom, ad.geom) AS geom
INTO sg_missing_
FROM AD ad, SG c

Running the same query on a larger country, such as Indonesia, took 15 hours, so there is something wrong with this approach. I was trying to play smart and get into the analysis only those admin boundaries that intersect the boundary box (envelope) of the country, but it still was running for 10 hours for Indonesia.

--preparing a single country bbox to use
SELECT ST_Envelope(geom) AS geom
INTO country_bbox
FROM countries_agg
where iso_cc = 'ID';

--preparing single country table and union of admin boundary geometries that intersect the given country envelope
WITH 
IND AS 
(SELECT iso_cc, geom FROM countries_agg WHERE iso_cc = 'ID'),
AD AS 
(SELECT ST_Union(ST_Intersection(IND.geom, AD.geom)) AS geom 
                                 FROM AD AD, country_bbox IND) 

--get geometry of country shape that is not covered by the admin boundaries
SELECT ST_Difference(c.geom, ad.geom) AS geom
INTO sg_missing
FROM AD ad, IND c

It is acceptable if this query would take a few days to run as long as all countries could be analysed in this time. How should such a query be written, if possible in PostGIS? I am open to solutions that involve other tools and software, it does not have to be done solely in PostGIS.

EDIT:

Running

SELECT  <id>,
        ST_Union(geom) AS geom
FROm    (
    SELECT  a.<id>,
            ST_Difference(a.geom, b.geom) AS geom
    FROM    <countries> AS a
    JOIN    <admins> AS b
      ON    a.geom && b.geom
    WHERE   ST_Intersects(a.geom, b.geom)
) q
GROUP BY
        id
;

does not report those country polygons that are located fully outside of admin boundaries, that is, they do not intersect any admin boundaries polygons at all. The bright green polygons should have been included along with the tiny red part of the country polygon that was included due to intersection with the admin boundaries polygon.

Answer 1

Unfortunately, all intersecting admin polygons need to be considered to create the combined ST_Difference with the country polygons:

SELECT  id,
        geom
FROM    (
    SELECT  a.id,
            CASE
                WHEN ST_Intersects(a.geom, ST_Collect(b.geom))
                    THEN ST_Difference(a.geom, ST_Collect(b.geom))
                ELSE a.geom
            END AS geom
    FROM    countries AS a
    LEFT JOIN
            admin AS b
      ON    a.geom && b.geom
    GROUP BY
            a.id
) q
WHERE GeometryType(geom) != 'GEOMETRYCOLLECTION';

Run this on simple Polygons in a new spatially indexed table and use a LEFT JOIN and according CASE to get and handle geometries that don't intersect. I also included a WHERE filter on GeometryType to exclude empty GEOMETRYCOLLECTION from ST_Difference.

---

My previous answer, with the query in your edit, did not consider multiple admin polygons to create the difference; unioning the individual results will actually lead to the original polygon being returned. It's exactly that what I was referring to when I said 'non-intersections are tricky'...

Answer 2

At the moment, in order to quickly solve your problem, I suggest combined approach (PostGIS, QGIS)

1)

CREATE TABLE tbl_union AS
SELECT table1.geom
FROM table1
UNION ALL
SELECT table2.geom
FROM table2;

2) Dissolving polygon:

CREATE TABLE tbl_union_byf AS
SELECT ST_Buffer(
         (ST_Dump(
           ST_Union(
             ST_Buffer(geom, 0.0000000001)))).geom, -0.0000000001) AS geom
FROM
tbl_union;

(John Powell :-)...)

3) Now use the QGIS Vector tool > Geoprocessing > Symmetrical Difference... and cut out tbl_union_byf table2 and this result.

In the beginning, I would still remove the polygons from the table2 from the processing that fall into the table1

That's how the decision of your question should have been, in my opinion ...

CREATE TABLE land_clip AS
SELECT DISTINCT a.id, ST_Difference(a.geom, b.geom) geom
FROM <countries> a JOIN <admins> b
ON ST_Overlaps(a.geom, b.geom)
GROUP BY a.id, a.geom, b.geom;