Extracting non-overlapping polygons in PostGIS

Question

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.

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.

Answer 1

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;