Le portail francophone de la géomatique

image95 · Mon 07 July 2025 11:10

Bonjour,
Suite à mon précédent message :
https://georezo.net/forum/viewtopic.php?id=135748 ,

je me permets d'en rédiger un nouveau afin de préciser/clarifier mon besoin. Merci aux personnes qui avaient pris le temps de me répondre. Grace à elles, j'ai pu avoir des premiers éléments de réponse.
Je travaille avec qgis et postgis. En entrée, j'ai 25 couches polygonales sur un grand territoire(intercommunalité). Dont une donnée très volumineuse (1 million d'objets). Les 24 autres sont beaucoup plus légères (au maximum une centaines d'objets).
Pour information, je précise que certaines de ces données polygonales sont en mode "entités multi parties" et d'autre en mode "morceaux/entités mono parties". J'imagine que celà peut avoir in fine une petite incidence sur la méthode/résultat.
Ces 25 fichiers .shp polygonaux ont des structures de données très variables, non homogènes/non harmonisées. Chaque couche a un champ "data_id" qui permet de définir pour chaque entité son appartenance à la couche. Par exemple, toutes valeurs du champ "data_id" de la première couche sont de valeur '1'. Pour la 2ème couche, les valeurs du champ sont '2', etc.

Mon objectif serait de pouvoir m'inspirer du moins en partie du geotraitement existant Qgis nommé Union multiple :
https://docs.qgis.org/3.40/fr/docs/user … n-multiple
Je colle une capture écran de la doc Qgis.

Mon but est d'avoir in fine un fichier:

- qui serait le résultat des croisements des 25 données d'entrée. Pour reprendre les termes de la doc qgis, que le traitement vérifie les chevauchements entre les entités au sein des 25 couches et crée des entités distinctes pour les parties qui se chevauchent et qui ne se chevauchent pas.
Ce geotraitement "union multiple" me convient pour les parties où il n'y pas de chevauchement (a, NULL ; b, NULL; c, NULL).
Pour les parties où il y a un chevauchement, le géotraitement union qgis créées autant d’entités de chevauchement identiques qu’il y a d’entités qui participent à ce chevauchement. En soi celà ne me dérange pas. Mais vu qu'in fine J'aimerais que dans le fichier résultat, un champ puisse permettre, pour chaque entité, de récupérer la liste des couches d'entrée qui participent/contribuent à cette entité résultat (Dans un but de récupérer l'origine/la source des données). Je me disais qu'il serait peut être préférable qu'il y ait qu'une seule entité créée par zone de chevauchement ?

- J'aimerais que dans le fichier résultat, un champ puisse permettre, pour chaque entité, de récupérer la liste des couches d'entrée qui participent/contribuent à cette entité résultat. Dans un but de récupérer l'origine/la source des données.

- Idéalement, un champ qui permette de récupérer le nombre (COUNT) de couches qui participent à cette entité (a minima 1 couche, au max 25 couches)

-Concernant cette fois la partie attributs/champ non géométrique, j'aimerais pouvoir spécifier la sélection/définir la liste champs que je souhaite vouloir conserver in fine. Je ne souhaite pas conserver la totalité, une partie des champs pour chacune des 25 couches d'entrée.

J'espère être suffisamment clair. J'imagine que cela n'est pas si facile.
J'imagine qu'il est recommandé de faire ce traitement sous postgis plutot que qgis ? Je peux si besoin importer mes 25 fichiers shp dans un bdd postgis. J'imagine aussi qu'il est important d'avoir à l'esprit que le mode "entités multi parties" / "morceaux/entités mono parties"
des couches d'entrée peut influer sur le résultat. Si je pars sur une requête postgis, je me disais qu'il aurait pu être intéressant de forcer le fait que toutes les entités soient en mode mono parties (avec la fonction postgis 'st_dump' ?).

Si une personne pense pouvoir m'aider? Un grand merci.

Dernière modification par image95 (Mon 07 July 2025 11:12)

Nicolas Ribot · Mon 07 July 2025 18:53

Bonsoir

• Importer tous les shp dans postgis
• faire l'union (st_union) de toutes les boundaries de vos polygones (st_boundary) de toutes les tables (utiliser le mot cle SQL UNION ALL pour cela)
• utiliser st_polygonize sur ce gros objet "union des boundaries" pour reconstruire les polygones découpés là ou les lignes s'intersectent
• Prendre les pointOnSurface de cette table de résultat et regarder dans quels polygones initiaux ils tombent: vous pourrez alors associer les attributs des polygones initiaux avec les polygones finaux.

Pour la gestion des champs que vous voulez, il faut mieux préparer cela en amont a partir des 25 tables et définir une structure commune: les champs que vous voulez garder, avec null quand pas de data

Nicolas

Nicolas Ribot · Mon 07 July 2025 19:12

Par exemple, avec 3 couches a, b, c contenant des cercles un peu comme dans l'exemple qgis (déso, je n'arrive pas à joindre des images aux messages ):

Code:

-- chargement data
-- ogr2ogr -f PostgreSQL PG:"dbname=nicolas" layer1.shp a
-- ogr2ogr -f PostgreSQL PG:"dbname=nicolas" layer2.shp b
-- ogr2ogr -f PostgreSQL PG:"dbname=nicolas" layer2.shp c
-- ...

-- union des geom et attributs qui vont bien dans une seule table
drop table if exists layer_union;
create table layer_union as
    select a.data_id, a.geom from a UNION ALL
    select b.data_id, b.geom from b UNION ALL
    select c.data_id, c.geom from c;

-- table de résultat: découpage des intersections entre lignes sous forme de polygones
drop table if exists res;
create table res as
    with tmp as (
    select st_union(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, 4326) as geom
from tmp;

-- id unique
alter table res add column res_id int generated always as identity primary key;

-- index sur les pointOnSurface: indispensable si grosse volumétrie
create index on res using gist (st_pointOnSurface(geom));
analyze res;

-- association des nouveaux polygones avec les originaux mis dans la table layer_union
select res.res_id, array_agg(l.data_id) as 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_id|id_ori |num_sources|
-- +------+-------+-----------+
-- |1     |{1,2}  |2          |
-- |2     |{2}    |1          |
-- |3     |{1,2,3}|3          |
-- |4     |{2,3}  |2          |
-- |5     |{1,3}  |2          |
-- |6     |{3}    |1          |
-- |7     |{1}    |1          |
-- |8     |{1,2}  |2          |
-- |9     |{2}    |1          |
-- |10    |{1,2}  |2          |
-- |11    |{1}    |1          |
-- |12    |{1,3}  |2          |
-- |13    |{3}    |1          |
-- |14    |{2,3}  |2          |
-- |15    |{3}    |1          |
-- |16    |{2}    |1          |
-- |17    |{2,3}  |2          |
-- +------+-------+-----------+

[img]https://i.postimg.cc/k682pP2j/Screenshot-2025-07-07-at-19-12-40.png[/img]

(maj: l'index spatial était défini sur la mauvaise table)

Nicolas

Dernière modification par Nicolas Ribot (Tue 08 July 2025 09:58)

Nicolas Ribot · Tue 08 July 2025 10:11

Pour préciser un peu le processus:

Le fait que certaines couches aient des polygones et d'autres des multipolygones ne gêne pas dans ce traitement (apres chargement dans postgis, suivant l'outil, il y a des fortes chances que toutes les couches soient des multipolygones de toute façon).

D'accord avec vous qu'il faut mieux n'avoir qu'un seul objet a chaque endroit, et pas des copies du meme objet. Ce traitement avec st_polygonize fait cela: il forme des polygones pour chaque contour fermé qu'il trouve dans la grosse multilinestring qu'on lui donne.
Avec le traitement final pour associer un point de chaque polygone final à un polygone initial, on retrouve toutes les infos dont on besoin: la ou les polygones initiaux qui ont conduit à ce polygone final.

Avec 1M de polygones, il se peut que vous ayez à découper le process spatialement pour qu'il passe, ou utiliser des clusters de géométries si les couches comportent des groupes séparés de polygones jointifs

Nicolas

Dernière modification par Nicolas Ribot (Tue 08 July 2025 10:25)

Nicolas Ribot · Tue 08 July 2025 10:50

J'ai testé avec les parcelles de Haute Garonne (1 343 528 polygones) vs 2 autres couches bidon que j'ai créées dans QGIS.
Au final layer_union (union dans une seule table des 3 tables d'input) contient 1 343 534 polygones.

La création de la table res:
[2025-07-08 10:41:45] 1,391,580 rows affected in 1 m 24 s 150 ms

La création de la table d'association donnant pour chaque polygone de résultat tous les polygones de layer_union concernés:
[2025-07-08 10:44:04] 1,391,580 rows affected in 10 s 771 ms

Plutôt performant donc. La perf va dépendre de la complexité des couches d'entrées (notamment, si vous avez de gros polygones étendus, il faut mieux les découper en amont avec st_subdivide (cf votre thread précédent)

(QGIS fait ca bien et vite également, j'ai testé sur les memes data)

Nicolas

Dernière modification par Nicolas Ribot (Tue 08 July 2025 11:43)

image95 · Tue 08 July 2025 11:28

Bonjour,
Je tenais à vous remercier grandement pour vos messages. Je suis en train d'essayer de reproduire votre proposition de méthode pour le moment juste pour 4 données sur les 25.
Je vous partage dans l'immédiat des messages qui parlent de ce cas d'usage. Des posts très proches de votre proposition de méthode. Vous les connaissez peut etre déjà :

https://lists.osgeo.org/pipermail/postg … 46812.html

https://blog.cleverelephant.ca/2019/07/ … rlays.html

https://dr-jts.github.io/postgis-patter … f-polygons

https://github.com/bcgov/harvest-restri … verlay.sql

https://lists.osgeo.org/pipermail/qgis- … 55609.html

Nicolas Ribot · Tue 08 July 2025 11:35

image95 a écrit:

Bonjour,
Je tenais à vous remercier grandement pour vos messages. Je suis en train d'essayer de reproduire votre proposition de méthode pour le moment juste pour 4 données sur les 25.
Je vous partage dans l'immédiat des messages qui parlent de ce cas d'usage. Des posts très proches de votre proposition de méthode. Vous les connaissez peut etre déjà :

https://lists.osgeo.org/pipermail/postg … 46812.html

https://blog.cleverelephant.ca/2019/07/ … rlays.html

https://dr-jts.github.io/postgis-patter … f-polygons

https://github.com/bcgov/harvest-restri … verlay.sql

https://lists.osgeo.org/pipermail/qgis- … 55609.html

Oui le blog de Paul Ramsey est une mine d'or.

image95 · Tue 08 July 2025 13:07

Encore un grand merci. J'arrive à reproduire votre proposition pour 4 données. Et comprends la logique globale. N'étant pas habitué à ce genre de procédure, il me faudra encore un peu de temps pour comprendre précisément la méthode. Je me permets une petite question. Si je dis pas de bêtise, concernant la dernière étape : la clause select permettant de récupérer le tableau d'association des nouveaux polygones avec les originaux mis dans la table layer_union / et de récupérer le comptage du nombre de couches source participant à la superposition.

Pour que ce ne soit pas uniquement un select mais des éléments ajoutés/présents dans la table "géométrique" existante "layer_union" ou bien même dans une nouvelle table géométrique (qu'on pourrait nommer "final"), il faudra ajouter une étape de jointure ?
Vous disiez par ailleurs avoir réussi une approche via qgis également ? Si ce n'est pas trop vous demander, je serai intéressé pour connaitre votre méthode avec qgis, même en quelques mots.
Encore un grand merci.

Dernière modification par image95 (Tue 08 July 2025 13:11)

Nicolas Ribot · Tue 08 July 2025 13:39

Oui la dernière étape dans mon exemple n'était qu'un select, mais dans votre cas, vous voulez surement stocker cela dans une table.

La fin de la procédure complète pourrait etre un truc comme ceci (je laisse les timings des requêtes qui ont tourné chez moi pour vous donner une idée des temps sur une couche 1.3M de polygones):

Code:

-- association des nouveaux polygones avec les originaux mis dans la table layer_union
-- On peut faire left join pour garder tout les polygones créés par st_polygonyze, 
-- ou un JOIN simple pour filtrer les polygones créés sans lien avec les polygones input (par ex des trous dans l'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;
-- [2025-07-08 10:44:04] 1,380,393 rows affected in 10 s 771 ms

create index on res_assoc(res_id);
analyse res_assoc;

-- ménage: on vire de res les polygones qui n'ont pas matché dans res_assoc:
-- ce sont des polygones représentant des trous dans les couches d'input
delete from res
where not exists (
    select null
    from res_assoc ra
    where ra.res_id = res.res_id
    );
-- [2025-07-08 11:06:57] 11,187 rows affected in 566 ms

-- table finale avec les nouveaux polygones et les infos des polygones source, sous forme de jointure: 
-- bien plus rapide de créer une table que de faire un update dans la table res (apres y avoir ajouté les colonnes voulues).
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;

Dans QGIS, j'ai juste lancé l'outil Union Multiple que vous mentionniez dans votre message.

Nicolas

Dernière modification par Nicolas Ribot (Tue 08 July 2025 13:40)

image95 · Tue 08 July 2025 15:28

Un grand merci. J'arrive à reproduire votre proposition.
Concernant Qgis, pour votre information, hier j'ai fait des tests avec union multiple. J'ai eu des messages d'erreurs pour un bon nombre de données.

Code:

Could not add feature with geometry type MultiLineString to layer of type MultiPolygon.

Quelques posts concernant ce message d'erreur :
https://github.com/qgis/QGIS/issues/44801

https://gis.stackexchange.com/questions … er-in-qgis

https://issues.qgis.org/issues/14846

Mais oui, je pense que c'est possible de faire le traitement géométrique avec Qgis. Mais à mon humble avis (je n'ai pas du tout votre expertise), pour reproduire votre procédure complète Postgis dans Qgis via l'interface simple user (récupération des couches et du nombre de couches contribuant/participant aux diverses superpositions notamment, j'ignore si c'est faisable.

Encore merci en tout cas !

image95 · Wed 09 July 2025 10:55

Bonjour,
Je suis en train de faire un test sur les 25 couches. J'ai bien réussi à faire la première étape des unions des geom/attributs des données sources au sein de la table layer_union. Pour l'étape de découpage des intersections entre lignes sous polygones.

Code:

drop table if exists res;
create table res as
    with tmp as (
    select st_union(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;

J'obtiens un message d'erreur :

Code:

ERROR:  lwcollection_construct: mixed dimension geometries: 0/3 

ERREUR:  lwcollection_construct: mixed dimension geometries: 0/3
SQL state: XX000

Il y a peut être des entités sans géométrie ? (0D) et/ou des entités à 3 dimensions ? Il me faut peut etre forcer la dimension 2D avec par exemple un ST_Force2D ?

EDIT : je confirme que ça passe avec un st_force2d :

Code:

-- table de résultat: découpage des intersections entre lignes sous forme de polygones
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;

Dernière modification par image95 (Wed 09 July 2025 11:43)

image95 · Thu 10 July 2025 12:36

Bonjour,
Le script fonctionne sans message d'erreur. Toutefois, j'ai l'impression qu'il y a encore de l'optimisation et/ou surtout de la fiabilisation/stabilisation à apporter. Car j'ai relancé plusieurs fois le même traitement. Et je n'arrive jamais exactement in fine au même nombre d'entités au sein de la table res_final. Lorsque je fais une symbologie categorisée sur le champ 'num_sources' et que j'affiche pour chaque classe le nombre d'entités rattachées, on peut voir en effet que ce nombre est variable d'un calcul à l'autre.

Dernière modification par image95 (Thu 10 July 2025 12:38)

Nicolas Ribot · Fri 11 July 2025 09:26

Bonjour

Hmm c'est curieux. Si vous pouvez isoler un petit dataset pour lequel vous pouvez reproduire ce comportement, ca vaut surement le coup de faire un ticket à Postgis.

Sans vos données, c'est difficile de dire ce qu'il se passe.
Ce sont des données publiques ?

Nicolas

image95 · Fri 11 July 2025 11:29

Bonjour,
Oui il s'agit essentiellement de data publiques. Oui je peux isoler et partager un petit dataset. D'accord je peux faire un ticket postgis.
Je partage une capture écran de mes tests d'hier. Capture écran permettant de comparer le nombre d'objets pour le test v1 et test v2 qui sont des tests identiques. On peut voir qu'il y a une différence qui s'observe dans la table res_final mais egalement pour la table res.
J'ai refait plusieurs tests hier. Avec toujours différents nombres d'objets pour les tables res et tables res_final. Avec souvent des différentiels plus importants que celui que je montre en capture écran.
Nombre d'entités pour chaque table :

Test 1 :
layer_union : 1026194
res : 1462661
res_assoc : 1462661
res_final :1462661

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

Dernière modification par image95 (Fri 11 July 2025 11:51)

image95 · Wed 23 July 2025 15:22

Delete

Dernière modification par image95 (Thu 24 July 2025 09:56)

Communauté	Ressources	L'emploi	Le marché	L'association	Fils RSS
Forum Blogs Wiki Planet Sites amis	Agenda Biblio Fiches AMO-CNIG Paris PCGIS Téléchargements	Offres d'emploi Banque de CV Formations Démarche Métiers	Geo-Entreprises Geo-Communiqué Appels d'offres	L'association Qui sommes-nous ? Nous aider Bienfaiteurs Partenaires Historique Statuts	Tout le site Job Géomatique Services GeoBlogs ...	Faire un don Mentions légales

Le portail francophone de la géomatique

Annonce

#1 Mon 07 July 2025 11:10

Union Multiple (25 couches polygonales)

#2 Mon 07 July 2025 18:53

Re: Union Multiple (25 couches polygonales)

#3 Mon 07 July 2025 19:12

Re: Union Multiple (25 couches polygonales)

Code:

#4 Tue 08 July 2025 10:11

Re: Union Multiple (25 couches polygonales)

#5 Tue 08 July 2025 10:50

Re: Union Multiple (25 couches polygonales)

#6 Tue 08 July 2025 11:28

Re: Union Multiple (25 couches polygonales)

#7 Tue 08 July 2025 11:35

Re: Union Multiple (25 couches polygonales)

#8 Tue 08 July 2025 13:07

Re: Union Multiple (25 couches polygonales)

#9 Tue 08 July 2025 13:39

Re: Union Multiple (25 couches polygonales)

Code:

#10 Tue 08 July 2025 15:28

Re: Union Multiple (25 couches polygonales)

Code:

#11 Wed 09 July 2025 10:55

Re: Union Multiple (25 couches polygonales)

Code:

Code:

Code:

#12 Thu 10 July 2025 12:36

Re: Union Multiple (25 couches polygonales)

#13 Fri 11 July 2025 09:26

Re: Union Multiple (25 couches polygonales)

#14 Fri 11 July 2025 11:29

Re: Union Multiple (25 couches polygonales)

#15 Wed 23 July 2025 15:22

Re: Union Multiple (25 couches polygonales)

Pied de page des forums