2

I am trying to perform what I thought ought to be a simple operation in PostGIS. I would like to use a set of circle polygons to split a set of square polygons and extract the resulting non-overlapping polygons. I have drawn a simplistic version of this below. I would like to extract non-overlapping polygons representing areas A, B, C, D, E, F, and G. Note that the datasets involved are large with tens of thousands of circles.

Problem representation

I have tried using ST_Intersection but that gives me overlapping polygons. I have tried reducing those overlapping polygons to lines and re-creating non-overlapping polygons (using ST_Boundary, ST_Union and ST_Polygonize) this works initially but then gives me Topology exceptions ('ERROR: GEOSUnaryUnion: TopologyException: found non-noded intersection between LINESTRING xx and LINESTRING xx'. I'm guessing due to issues with excess coordinate places) when I start running it for very large numbers of circles. I have tried getting rid of the topology errors using ST_MakeValid, ST_Node, ST_Buffer(geom,0), ST_SnapToGrid etc. as suggested elsewhere on here, all with no success.

So can anyone provide an alterative approach to solving this, seemingly simple problem? I get the feeling I'm missing something obvious here and getting too lost in trying to unpick the topology exception when I think there ought to be a simplier solution to the intersect.

6
  • 2
    Have a look at this blog post - the technique can be adapted: blog.cleverelephant.ca/2019/07/postgis-overlays.html
    – dr_jts
    Aug 2 at 18:07
  • Have you tried st_difference? To use it, you would union all the geom for the circles, then set the geom on the polygons to st_difference(poly.geom, allcircles.geom)
    – jbalk
    Aug 2 at 19:03
  • something like select st_difference(poly.geom, allcircles.geom), poly.id from your_poly_table poly, (select st_union(geom) from your_circles_table) allcircles;
    – jbalk
    Aug 2 at 19:06
  • Sorry - I misunderstood the question. I think st_split may actually be useful. Something like this select st_split(poly.geom, circlelines.geom), poly.id from your_poly_table poly, (select st_union(st_exteriorring(geom)) from your_circle_table) circlelines;
    – jbalk
    Aug 2 at 19:13
  • I forgot to note that you'll need to either dump or extract the resulting geometry collection after the split - st_collectionextract(st_split(poly.geom, allcircles.geom), 3) or (st_dump(st_split(poly.geom, allcircles.geom))).geom
    – jbalk
    Aug 2 at 22:58

1 Answer 1

0

Install the PostGIS Addons and run this query:

WITH rectangletable AS (
  SELECT 1 id, ST_GeomFromText('POLYGON((0 0, 0 3, 2 3, 2 0, 0 0))') geom
  UNION ALL
  SELECT 2 id, ST_GeomFromText('POLYGON((3 0, 3 3, 5 3, 5 0, 3 0))') geom
  UNION ALL
  SELECT 3 id, ST_GeomFromText('POLYGON((6 0, 6 3, 8 3, 8 0, 6 0))') geom
), circletable AS (
  SELECT 1 id, ST_Buffer(ST_GeomFromText('POINT(4 3)'), 2.5) geom
  UNION ALL
  SELECT 2 id, ST_Buffer(ST_GeomFromText('POINT(4 0)'), 2.5) geom
)
, splitted AS (
  SELECT unnest(ST_SplitAgg(a.geom, b.geom)) geom
  FROM rectangletable a,
       circletable b
  WHERE ST_Equals(a.geom, b.geom) OR
        ST_Contains(a.geom, b.geom) OR
        ST_Contains(b.geom, a.geom) OR
        ST_Overlaps(a.geom, b.geom)
  GROUP BY ST_AsEWKB(a.geom)
)
SELECT DISTINCT ON (ST_AsEWKB(geom)) geom
FROM splitted;

If it works as expected, replace the two tables with your own tables. That would become simply:

WITH splitted AS (
  SELECT unnest(ST_SplitAgg(a.geom, b.geom)) geom
  FROM rectangletable a,
       circletable b
  WHERE ST_Equals(a.geom, b.geom) OR
        ST_Contains(a.geom, b.geom) OR
        ST_Contains(b.geom, a.geom) OR
        ST_Overlaps(a.geom, b.geom)
  GROUP BY ST_AsEWKB(a.geom)
)
SELECT DISTINCT ON (ST_AsEWKB(geom)) geom
FROM splitted;
1
  • Thank you Pierre Racine. I will look into whether this is possible on the linux setup I am using.
    – Ddee
    Aug 3 at 14:43

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.