Overlay of overlapping polygonal datasets : postgis bug/lack of reliability/stabilization ?

[celati Laurent]

 Good afternoon,
I'm taking the liberty to write a message beacause i feel a bug with
postgis. As part of a  "Overlay of overlapping polygonal datasets"
usecase.
>From 28 postgis polygon vectors.I try to node and polygonize to make a flat
coverage, and then join back to the parent layers to determine
attribution (including counts).

For information, i work with :

> "POSTGIS=""3.5.2 3.5.2"" [EXTENSION] PGSQL=""170""
> GEOS=""3.13.0-CAPI-1.19.0"" SFCGAL=""SFCGAL 1.5.2, CGAL 5.6.1, BOOST
> 1.84.0"" PROJ=""8.2.1


I succeed in executing a postgis script without error message. It is
working.
However, i notice a lack of  reliability/stabilization. Because when I've
rerun the same process several times, i never end up with exactly the same
number of features in my intermediate/final result tables.
I'm taking the liberty to share you the sql script.
The screenshot compares the number of objects for test v1 and test v2,
(which are identical tests). We can see that there is a difference in the
res_final table, but also in the res table. I ran several tests agai. Still
with different numbers of objects for the res and res_final tables. Often
with larger differences than the one shown in the screenshot.

Number of entities for each table:
*Test 1*:
layer_union table: 1026194
res table : 1462661
 res_assoc table : 1462661
 res_final table : 1462661

*Test 2*
 layer_union table : 1026194
res table 1462645
 res_assoc table : 1462645
 res_final table : 1462635

I share below/and attach the script :

--Import all shp in postgis db
-- union ALL des geom et attributs des 28 data sources dans une seule table
drop table if exists layer_union;
create table layer_union as
select inrae_2014.data_id, inrae_2014.geom from inrae_2014 UNION ALL
select aesn_2006.data_id, aesn_2006.geom from aesn_2006 UNION ALL
select aesn_2019.data_id, aesn_2019.geom from aesn_2019 UNION ALL
--(...)etc.

-- res table
drop table if exists res;
create table res as
with tmp as
(select st_union(ST_Force2D(st_boundary(geom))) as geom from layer_union
)
select (st_dump(st_collectionextract(st_polygonize(geom), 3))).path[1] as
id,
       (st_dump(st_collectionextract(st_polygonize(geom),
3))).geom::geometry(polygon, 2154) as geom
from tmp;

-- res table id unique
alter table res add column res_id int generated always as identity primary
key;
-- res table index on pointOnSurfacee
create index on res using gist (st_pointOnSurface(geom));
analyze res;

-- res_assoc table
--JOIN simple for filter polygons without link with input polygons (for
instance : holes for data input)
drop table if exists res_assoc;
create table res_assoc as
select res.res_id, array_agg(l.data_id) as data_id_ori, count(distinct
l.data_id) as num_sources
from res join layer_union l on st_contains(l.geom,
st_pointonsurface(res.geom))
group by res.res_id;
-- res_assoc table : index creation
create index on res_assoc(res_id);
analyse res_assoc;

----cleaning: we remove from the res table the polygons that did not match
in res_assoc:
-- these are polygons representing holes in the input layers
delete from res
where not exists (
select null
from res_assoc ra
where ra.res_id = res.res_id);


-- -- Final table with the new polygons and the source polygon information,
as a join:
-- Much faster to create a table than to update the res table (after adding
the desired columns).
drop table if exists res_final;
create table res_final as
select ra.res_id, data_id_ori, num_sources, geom::geometry(polygon, 2154)
as geom
from res_assoc ra join res r on ra.res_id = r.res_id;

Thanks so much
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20250711/5abc3b12/attachment.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Compare_res_final.jpg
Type: image/jpeg
Size: 16925 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20250711/5abc3b12/attachment.jpg>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: script_GB_comment.sql
Type: application/octet-stream
Size: 2095 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20250711/5abc3b12/attachment.obj>