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
C) (can be multipolygon/polygon), transparent blue polygons (three rows in table in this case).
1is fully covered by the admin boundary
2is not fully covered by admin boundaries (
C), so I would like to get this reported. If possible, the hatched area of country
2could 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.
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.