2

I need to get the Difference of a single large-area polygon (input layer) and a very numerous set of polygons (overlay layer). The overlay layer has ~500k polygons, of which ~200k intersect the input layer. I am trying to use PostGIS to do this in the hope that it can be much faster than e.g. running a Difference algorithm in QGIS (which is unfeasibly slow).

I've got the basic approach to do this with ST_Difference from here, i.e. just run ST_Difference on a left join on ST_Intersects (using Coalesce to avoid losing polygons from the input layer that don't intersect with the overlay layer - although as it happens in this particular case there are no such polygons to lose as there is only a single input polygon which does intersect).

So my current query is:

SELECT COALESCE(ST_Difference(input.geom,overlay.geom),input.geom)
FROM input LEFT JOIN overlay ON ST_Intersects(input.geom,overlay.geom)

Unsurprisingly this query is taking a very long time, as it's calculating ~200k differences.

I've read about speeding up intersections enormously by skipping ST_Intersection where polygons are entirely contained within one another e.g. here and here. I.e. using ST_CoveredBy / ST_Within / ST_Contains etc. (I haven't completely sussed out the difference between these yet.)

I feel like I should be able to do something similar here, as the vast majority of the overlay polygons are entirely contained within the input polygon. But I can't get my head around the query code to use to do the same with ST_Difference.

Improve this question
11
  • 1
    Years back this had been a scenario that made me think outside the box (and PostGIS/GIS can be a sticky box...); for some client, I had to do this with +7.5M highly complex, but small polygons on an even more complex continental Eurasia Polygon. Turned out it takes less than 5 seconds to add the formatted coordinate strings of them all as inner rings to the WKT of the large Polygon...
    – geozelot
    Commented Mar 17, 2021 at 11:39
  • 2
    The WKT format is a string representation of a geometry, and a Polygon is defined as a set of rings (outer, plus any amount of inner rings) which are made up by coordinates; making holes into a Polygon would thus mean to simply add some rings. Anyways, subdividing is the ultimate improvement, it wasn't clear from the question if your Polygon needs to stay as one.
    – geozelot
    Commented Mar 17, 2021 at 13:54
  • 2
    Do you want to remove ALL the overlay polygons from the input polygon, or compute the difference of the input and EACH overlay polygon?
    – dr_jts
    Commented Mar 17, 2021 at 17:34
  • 1
    @dr_jts I want to remove all the overlay polygons from the input polygon. You've seen (and helpfully commented) on a subsequent question of mine, which has led to getting the solution. I learnt a key step from this question though - the critical importance of subdividing the input so that bounding box checks can speed things up massively. Thanks.
    – MDB
    Commented Mar 17, 2021 at 22:51
  • 1
    Subdividing the input will definitely increase performance, as noted. I'd be a bit worried that numeric roundoff would lead to small anomalies in the result, however (i.e. "gores" where the added intersections are not quite identical across two adjacent subdivided pieces). But this should happen rarely, and may not matter for your use case.
    – dr_jts
    Commented Mar 17, 2021 at 23:58

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.