4

I want to get the difference between a bounding box (single rectangle) and a collection of polygons (a bunch of complex geometries) within that bounding box.

That is, say my bounding box encompasses the USA, and my polygons are the urban centers polygons, I want to get the geometry corresponding to the regions within my bounding box that are NOT a urban center.

I can easily get the urban centers within my bounding box:

WITH bb_uc AS (
    SELECT urban_centers.geometry as geometry
    FROM urban_centers
    WHERE ST_INTERSECTS(urban_centers.geometry,
                        ST_MakeEnvelope(-124.7844079, 24.7433195, -66.9513812, 49.3457868, 4326))
)
SELECT bb_uc.geometry
FROM bb_uc;

This query runs in around 10 seconds in my case, and returns ~100k rows with polygon geometries.

Now, when I try to get the "complement" of this spatial dataset with relation to my bounding box:

WITH bb_uc AS (
    SELECT urban_centers.geometry as geometry
    FROM urban_centers
    WHERE ST_INTERSECTS(urban_centers.geometry,
                        ST_MakeEnvelope(-124.7844079, 24.7433195, -66.9513812, 49.3457868, 4326))
)
SELECT ST_DIFFERENCE(ST_MakeEnvelope(-124.7844079, 24.7433195, -66.9513812, 49.3457868, 4326), ST_UNION(bb_uc.geometry)) as geometry
FROM bb_uc

This query takes forever to run (I stopped manually after 3h). I had not imagined that this operation would be that complex (basically, apply a mask to my bounding box)

As an aside, note that I'm using an "unnecessary" WITH statement because I sometimes add a WHERE clause in there on my results.

What can I do to dramatically increase my performance to get to the "complementary" geometry of my intersection? ST_Difference man page mentions not to use geometryCollection, hence my use of ST_Union, which is the operation taking forever. However, I can't think of a way around this.

postgis_version: 2.5 USE_GEOS=1 USE_PROJ=1 USE_STATS=1

11
  • 1
    those polygons should not overlap, right, since they represent connected urban areas? in that case, ST_Collect will get you the same MultiPolygon without costly dissolving. you might get better results if you create a table with that MultiPolygon, but that depends on a few things, so try it. in the end, those operations on large and complex geometries are computational intensive, so do expect it to run a while. interestingly, if those polygons don't overlap, one could easily construct a valid OGC defined Polygon by simply adding them as inner rings to the coordinate array...
    – geozelot
    Commented Jul 25, 2019 at 19:55
  • one question, though: why would you want that polygon? I could only come up with very few scenarios where you actually can't find a more performant substitute...
    – geozelot
    Commented Jul 25, 2019 at 19:58
  • Well, the question to answer is: "Show me all the locations that are not urban centers", hence why that particular polygon is important. This polygon will then be intersected with some other datasets. If there is a simple work-around I'm missing, and it's very possible, I'm all ears! ST_Collect could work. I'm trying it now but am already disappointed by the performance (still running after 5 minutes). Commented Jul 25, 2019 at 20:05
  • ...and the snappy answer would be: the white area between the black polygons? ,) i.e. on a map you can work with contrasting background color. and get e.g. the area by subtraction.
    – geozelot
    Commented Jul 25, 2019 at 20:12
  • I definitely agree. But now, consider this scenario: "Show me all the locations that are not urban centers, and then show me which of these locations have a population density between 2 and 4 people/km2". In this case, you need to obtain the geometry for the not-urban-centers. Commented Jul 25, 2019 at 20:27

1 Answer 1

2

SO,

If anyone has problems resolving this issue, proceed as follows:

my source geodata is a polygon called poly_extent and administrative polygons without holes called adm_polygons see Figure 1.

enter image description here

Run the script

WITH 
    tbla AS (SELECT (ST_DumpPoints(geom)).geom FROM adm_polygons),
    tblb AS (SELECT ((ST_Dump(ST_VoronoiPolygons(ST_Collect(geom)))).geom) geom FROM tbla),
    tblc AS (SELECT ST_Intersection(a.geom, b.geom) geom FROM tblb a JOIN poly_extent b ON ST_Intersects(a.geom, b.geom)),
    tbld AS (SELECT ST_Difference(a.geom, b.geom) geom FROM tblc a JOIN adm_polygons b ON ST_Intersects(a.geom, b.geom) AND ST_Overlaps(a.geom, b.geom))
    SELECT ST_Union(geom) geom FROM tbld;

See the result in Figure 2.

enter image description here

Good luck to everyone :-),

This script is called - ST_CarvesPolygons

So, the sequel follows...

dr_jts (Martin Davis https://github.com/dr-jts) paid attention to the architecture of the geoinstrument given in the answer above and started giving me "signals" :-)...

So I had to think again and this is what we have got so far :-).

The architecture of the geoinstrument could most likely look like this:

WITH
  tbla AS (SELECT ST_Union(geom) geom FROM poly_table1),
  tblb AS (SELECT ST_Union(geom) geom FROM poly_table2)
          SELECT DISTINCT ST_Difference(a.geom, b.geom) geom FROM tbla a JOIN tblb b ON ST_Intersects(a.geom, b.geom) 
  AND ST_Overlaps(a.geom, b.geom);

The geospatial function might look like this:

CREATE OR REPLACE FUNCTION ST_CarvesPolygons(geom1 GEOMETRY, geom2 GEOMETRY)
RETURNS GEOMETRY AS
$BODY$
WITH
    tbla AS (SELECT ST_Union($1) geom),
    tblb AS (SELECT ST_Union($2) geom)
            SELECT DISTINCT ST_Difference(a.geom, b.geom) geom FROM tbla a JOIN tblb b ON ST_Intersects(a.geom, b.geom) 
    AND ST_Overlaps(a.geom, b.geom);
$BODY$
LANGUAGE SQL

Running the function so far looks like this, and don't be intimidated by the cumbersome syntax yet, the main thing in our case is logic and the expected result:

SELECT ST_CarvesPolygons(ST_Union(a.geom), ST_Union(b.geom)) geom FROM poly_table1 a JOIN poly_table2 b ON true;

If you want to get individual objects in the table run the script:

SELECT (ST_Dump(geom)).geom geom FROM (SELECT
ST_CarvesPolygons(ST_Union(a.geom), ST_Union(b.geom)) geom FROM poly_table1 a JOIN poly_table2 b ON true) foo;

Now you can swap the polygon tables and get the expected results. So far the speed of execution of the new function differs from the corresponding query with the WITH construct, but perhaps in the long run we can deal with this circumstance...

Don't try to replace ST_Union() with ST_Collect() in this geoinstrument, because the speed of execution of the operation may drop dozens of times...

Thanks again Martin Davis for your "beacon" that made me revisit this geo-tool...

One bright head is good, but two (one bright and one cheerful) is even better :-)...

Original geospatial solutions...

Translated with www.DeepL.com/Translator (free version)

4
  • This solution is very clever, but I don't think it handles overlapping input polygons. (The difference step removes each separate overlapping polygon from a Voronoi cell, instead of the union of the overlapping polygons).
    – dr_jts
    Commented Sep 20, 2021 at 19:08
  • @dr_jts, I agree, I was just playing around with the tools and thought this design could handle the required task, and also noticed that using the ST_VoronoiPolygons function in PostGis is very expensive and so left it as is, I even wanted to create a function from this design, but it didn't work :-). Commented Sep 20, 2021 at 19:51
  • Yes, I was thinking ST_VoronoiPolygons is quite expensive for doing this. I think a better approach might be to grid the extent box.
    – dr_jts
    Commented Sep 20, 2021 at 20:29
  • Thus, we need to think a little bit about how to optimize this design... Commented Sep 20, 2021 at 20:58

Your Answer

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

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