Vjylku

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.

I can't test right now, but I'd say try with

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

;

Run this on simple Polygons in a new spatially indexed table. I used an additional && join, see what the planner makes of that; it might help or not. ST_Difference and ST_Within/ST_Contains should do the same checks, so I included none of the latter as additional spatial relation check.

I think that unioning a possibly low number of smaller difference polygons should be a lot faster than a per se union of all intersecting admin polygons. It shoud also get rid of empty GeometryCollections that can result from ST_Difference.

Note that this query will return areas that are not covered by any admin region, so if a province of a country overlaps with another country, that part will count as covered. You'd need to have attribute reference to avoid this.