<div dir="ltr">
<div dir="ltr"><div>Good afternoon,</div><div>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. <div>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).<br><div><br></div></div></div><div>For information, i work with : <br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">"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 </blockquote><div> </div></div><div>I succeed in executing a postgis script without error message. It is working. </div><div>However, i notice a lack of
<span lang="en"><span><span>reliability/stabilization.</span></span></span><span lang="en"> <span><span>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. <br>I'm taking the liberty to share you the sql script. </span></span></span></div><div>
<div dir="ltr"><span lang="en"><span><span>The screenshot compares the number of objects for test v1 and test v2, (which are identical tests).</span></span> <span><span>We can see that there is a difference in the res_final table, but also in the res table.</span></span><span><span>
I ran several tests agai.</span></span> <span><span>Still with different numbers of objects for the res and res_final tables.</span></span> <span><span>Often with larger differences than the one shown in the screenshot.</span></span><span><span> </span></span></span></div><div dir="ltr"><span lang="en"><span><span><br></span></span></span></div><div dir="ltr"><span lang="en"><span><span>Number of entities for each table:</span></span></span></div><div dir="ltr"><span lang="en"><span><span><b>Test 1</b>:</span></span></span></div><div dir="ltr"><span lang="en"><span><span>layer_union table: 1026194 </span></span></span></div><div dir="ltr"><span lang="en"><span><span>res table : 1462661</span></span></span></div><div dir="ltr"><span lang="en"><span><span> res_assoc table : 1462661</span></span></span></div><div dir="ltr"><span lang="en"><span><span> res_final table : 1462661 </span></span></span></div><div dir="ltr"><span lang="en"><span><span><br></span></span></span></div><div dir="ltr"><span lang="en"><span><span><b>Test 2</b></span></span></span></div><div dir="ltr"><span lang="en"><span><span><b> </b>layer_union table : 1026194 </span></span></span></div><div dir="ltr"><span lang="en"><span><span>res table 1462645</span></span></span></div><div dir="ltr"><span lang="en"><span><span> res_assoc table : 1462645</span></span></span></div><div dir="ltr"><span lang="en"><span><span> res_final table : 1462635</span></span></span></div><div dir="ltr"><span lang="en"><span><span><br></span></span></span></div><div><span lang="en"><span><span>I share below/and attach the script : <br></span></span></span><br><blockquote>--Import all shp in postgis db<br>-- union ALL des geom et attributs des 28 data sources dans une seule table<br>drop table if exists layer_union;<br>create table layer_union as<br>select inrae_2014.data_id, inrae_2014.geom from inrae_2014 UNION ALL<br>select aesn_2006.data_id, aesn_2006.geom from aesn_2006 UNION ALL<br>select aesn_2019.data_id, aesn_2019.geom from aesn_2019 UNION ALL<br>--(...)etc.<br><br>-- res table<br>drop table if exists res;<br>create table res as<br>with tmp as <br>(select st_union(ST_Force2D(st_boundary(geom))) as geom from layer_union<br>)<br>select (st_dump(st_collectionextract(st_polygonize(geom), 3))).path[1] as id,<br> (st_dump(st_collectionextract(st_polygonize(geom), 3))).geom::geometry(polygon, 2154) as geom<br>from tmp;<br><br>-- res table id unique<br>alter table res add column res_id int generated always as identity primary key;<br>-- res table index on pointOnSurfacee<br>create index on res using gist (st_pointOnSurface(geom));<br>analyze res;<br><br>-- res_assoc table<br>--JOIN simple for filter polygons without link with input polygons (for instance : holes for data input)<br>drop table if exists res_assoc;<br>create table res_assoc as<br>select res.res_id, array_agg(l.data_id) as data_id_ori, count(distinct l.data_id) as num_sources<br>from res join layer_union l on st_contains(l.geom, st_pointonsurface(res.geom))<br>group by res.res_id;<br>-- res_assoc table : index creation <br>create index on res_assoc(res_id);<br>analyse res_assoc;<br><br>----cleaning: we remove from the res table the polygons that did not match in res_assoc:<br>-- these are polygons representing holes in the input layers<br>delete from res<br>where not exists (<br>select null<br>from res_assoc ra<br>where ra.res_id = res.res_id);</blockquote><br>-- -- Final table with the new polygons and the source polygon information, as a join:<br>-- Much faster to create a table than to update the res table (after adding the desired columns).<br>drop table if exists res_final;<br>create table res_final as<br>select ra.res_id, data_id_ori, num_sources, geom::geometry(polygon, 2154) as geom<br></div><div>from res_assoc ra join res r on ra.res_id = r.res_id;</div><div><br></div><div>Thanks so much</div></div></div><div class="gmail-yj6qo gmail-ajU"><div id="gmail-:1dh" class="gmail-ajR" role="button" tabindex="0" aria-label="Afficher le contenu abrégé" aria-expanded="false"><img class="gmail-ajT" src="https://ssl.gstatic.com/ui/v1/icons/mail/images/cleardot.gif"></div></div>
<br></div>