[postgis-users] Complex Overlays

[Eloi Ribeiro]

Hi all,

Did anyone successfully executed the following procedure http://blog.cleverelephant.ca/2019/07/postgis-overlays.html with complex geometries (polygons with holes)?

Below the SQL code I am using (read failing). I still get duplicated (overlapped) polygons and some others are just missing in the final layer.
Any help would be very much appreciated.

Sample data that I'm using: https://eloiribeiro.eu/s/bkoWYyNZJMT8bpW

CREATE SEQUENCE boundaries_gid_seq;

CREATE TABLE boundaries AS
SELECT nextval('boundaries_gid_seq') AS gid,
gid AS gid_old,
(ST_DumpRings(geom)).path AS path,
((ST_DumpRings(geom)).geom) AS geom
FROM complex_circles; -- layer with complex polygons, with holes

ALTER TABLE boundaries ADD PRIMARY KEY(gid);
CREATE INDEX boundaries_geom_idx ON boundaries USING gist(geom);

CREATE SEQUENCE polys_seq;

CREATE TABLE polys AS
SELECT nextval('polys_seq') AS gid,
((ST_Dump(ST_Polygonize(geom))).geom) AS geom
FROM boundaries;

ALTER TABLE polys ADD PRIMARY KEY(gid);
CREATE INDEX polys_geom_idx ON polys USING gist(geom);

ALTER TABLE polys ADD COLUMN count INTEGER DEFAULT 0;
ALTER TABLE polys ADD COLUMN codes text;

UPDATE polys
SET count = p.count,
codes = p.array_agg
FROM (
SELECT p.gid, count(*), array_agg(code)
FROM polys p
JOIN complex_circles c
ON ST_Contains(c.geom, ST_PointOnSurface(p.geom))
GROUP BY p.gid
) AS p
WHERE p.gid = polys.gid;

DELETE FROM polys WHERE count = 0;

Software version:
SELECT version();
"PostgreSQL 12.16 ..."
SELECT postgis_full_version();
"POSTGIS=""3.1.9 95a641b"" [EXTENSION] PGSQL=""120"" GEOS=""3.10.3-CAPI-1.16.1"" PROJ=""7.2.1"" LIBXML=""2.9.1"" LIBJSON=""0.11"""

Eloi
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20230830/d2eb12bc/attachment.htm>