7

I have a set of 1000 polygons that overlap and I want to count the areas of overlap. I am using this blog post's instructions to create multilines from the 1000 polygons then use this line file to generate polygons for the overlapping areas and then count the overlaps. http://boundlessgeo.com/2014/10/postgis-training-creating-overlays/

This works fine for datasets of 100 or fewer polygons but hangs on anything more than that on the first step. It seems like the ST_Union is what's so slow, but ST_Collect seems not to work here because it generates a multiline with far fewer nodes and doesn't create non-overlapping polygons when used in the second step. It's too bad because ST_Collect is so fast - but doesn't give the same result. Anyone have ideas on how I can modify this process to work faster on a dataset of 1000 polygons?

Step 1

CREATE TABLE boundaries_polygons1000 AS
SELECT ST_Union(ST_ExteriorRing(wkb_geometry)) AS geom
FROM polygons1000;

Step 2

CREATE SEQUENCE polyseq_polygons1000;
CREATE TABLE polys_polygons1000 AS
SELECT nextval('polyseq_polygons1000') AS id, (ST_Dump(ST_Polygonize(geom))).geom AS geom
FROM boundaries_polygons1000;

Step 3

ALTER TABLE polys_polygons1000 ADD COLUMN count INTEGER DEFAULT 0;
UPDATE polys_polygons1000 set count = p.count
FROM (
  SELECT count(*) AS count, p.id AS id  
  FROM polys_polygons1000 p 
  JOIN polygons1000 c 
  ON ST_Contains(c.wkb_geometry, ST_PointOnSurface(p.geom)) 
  GROUP BY p.id
) AS p
WHERE p.id = polys_polygons1000.id; 
  • A hunch: cascaded union doesn't work with linestrings (yet?) in PostGIS, so it's using unary union which is much slower (you've probably seen the GEOSUnaryUnion warning if you hit CTRL-C during the query). I'll try to think of an answer if I have time. – Rob Skelly Aug 21 '15 at 23:27
  • Thanks, Rob - that is what I have seen if I cancel during the query GEOSUnaryUnion. Any advice you have would be helpful. – Lee R. Aug 22 '15 at 0:08
6

It seems you can use ST_Node to do this faster than ST_Union. I used this query, which I think gives the same result:

CREATE TABLE boundaries AS
SELECT (ST_Dump(ST_Node(ST_Collect(ST_ExteriorRing(geom))))).geom AS geom
FROM circles;

What I'm doing is collecting the exterior rings and then noding them, which splits them into individual segments, ending at the intersections. With 1000 circles, it completes in well under a minute.

The final query in step 3 is also slow. I sped it up by building indices on the geometries for both tables and using a bounding box filter (&&) in the ON clause.

Judging by the image, it works.

Image of circles and counts.

Original answer

Maybe I misunderstand, but if you only want to count the overlaps, why make any new geometries? You can get a count of intersections:

SELECT count(*) 
FROM geoms a, geoms b 
WHERE a.gid < b.gid AND ST_Intersects(a.geom ,b.geom)

This basically checks the intersection between every geometry and every other geometry, but prevents an intersection from being counted twice (a && b, but not b && a). It's still n^2 in the worst case, but it should be fine. It will use the index.

  • Thanks. I'm trying to find all overlapping parts of polys and count how many times each part is overlapped. I think this from the blog post above better explains what I'm trying to do, "Calculating the overlapping parts of a pair of shapes is easy, using the ST_Intersection() function in PostGIS, but that only works for pairs, and doesn’t capture the areas that have no overlaps at all. How can we handle multiple overlaps and get out a polygon set that covers 100% of the area of the input sets? By taking the polygon geometry apart into lines, and then building new polygons back up." – Lee R. Aug 21 '15 at 21:23
  • Apologies -- that's what I get for not reading carefully! – Rob Skelly Aug 21 '15 at 21:33
  • Edited with an answer to the actual question ;) – Rob Skelly Aug 23 '15 at 8:31
  • Thank you for the answer. This creates multiple features instead of just one multiline like in the boundlessgeo.com example and these won't give the same result in the third step. I realized that step one works in under a minute on 1000s of circles on a new macbook, so this is an issue with my particular starting polygons. I'll spend some time understanding why this is slow on my polygons and/or why this might be slow on polygons that aren't circles. – Lee R. Aug 25 '15 at 14:02
  • If you run the example with 100 circles, and then run it again substituting my query, you'll see that the final result is identical even though the count on the boundaries table is >1 while it's 1 for the original algorithm (because the dump in step 2 does nothing -- it's already dumped). If you'd be willing to share a sample of your data, I can try too. – Rob Skelly Aug 25 '15 at 15:48

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.