3

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.

enter image description here

  • 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.

enter image description here

  • how many geometries are there in both tables? and can those admin polygons overlap more than one country? getting non-intersections/differences is always tricky, but if there'd be a 1 to 1 intersection relation between each country and admin area, one might get around expensive unions. also, having only non-multi geometries boosts the index search and vertex traversal significantly. – ThingumaBob Feb 21 at 11:01
  • Thanks for getting on it. Country table has 240 features. Admin boundary table has 164000 features. Admin polygons can overlap more than one country, so unfortunately, there isn't 1:1 between each country and admin area. Admin areas have 70% polygons and 30% multipolygons. Thanks for the tip to explode those into polygons. – Alex Tereshenkov Feb 21 at 11:14
  • You could try to explode the multipolygon, then check with st_contains+ st_intersects which polygon are contained by which admin area and finally use a group by to regroup the information by country. Don't forget to create a spatial index – obchardon Feb 21 at 11:27
  • @obchardon, thanks for posting. I don't think I follow you on check with st_contains+ st_intersects which polygon are contained by which admin area. Could you please elaborate briefly - or even post an answer? – Alex Tereshenkov Feb 21 at 11:36
  • so you actually have a lower admin level than 0, e.g. provinces? shouldn't there be some sort of reference attribute between both tables? i.e. admin polygons A & B each are provinces of country 1? – ThingumaBob Feb 21 at 11:37
2

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'...

1

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;
  • And in my experience there is no need to use only one special software, it is best to combine their application in practice ...: -) ... – Cyril Feb 22 at 8:55
  • Thanks for posting! Could you please briefly update your answer elaborating on what exactly you meant by In the beginning, I would still remove the polygons from the table2 from the processing that fall into the table1? Thanks! – Alex Tereshenkov Feb 22 at 11:32
  • 1
    In your first drawing, I discovered that you have green (table2) polygons and some of them come in completely in blue (table1) and, as a result, they should not be used in further processing, i.e. we are only interested in those polygons that hit the blue polygon border (table1) ... maybe I misunderstood Figure 1 ... – Cyril Feb 22 at 11:39
  • 1
    But after you added Figure 2, now you need to leave the green polygons that touch the blue (find external fragments - you now have the answer to these actions), and add those green polygons that do not touch and do not fall into the blue polygons and merge them if necessary ... – Cyril Feb 22 at 11:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.