3

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

  • 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... – ThingumaBob Jul 25 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... – ThingumaBob Jul 25 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). – William Abma Jul 25 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. – ThingumaBob Jul 25 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. – William Abma Jul 25 at 20:27
1

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 ta as (SELECT (ST_DumpPoints(geom)).geom FROM adm_polygons), tb as (SELECT ((ST_Dump(ST_VoronoiPolygons(ST_Collect(geom)))).geom) bgeom FROM ta), tc as (SELECT ST_Intersection(a.bgeom, b.geom) cgeom FROM tb a JOIN poly_extent b ON ST_Intersects(a.bgeom,b.geom)), td as (SELECT ST_Difference(c.cgeom, d.geom) dgeom FROM tc c JOIN adm_polygons d ON ST_Intersects(c.cgeom, d.geom) AND ST_Overlaps(c.cgeom, d.geom)) SELECT ST_Union(dgeom) geom FROM td; `

See the result in Figure 2.

enter image description here

Good luck to everyone :-),

OS ...

P.S. I will name this tool - ST_CarvesPolygons.

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.