Le portail francophone de la géomatique

Bonsoir,

Je vais tester l'option décrite par Nicolas, à savoir ajuster sur base des paires de points proches. Je ferai un retour des résultats.

A toute fin utile, j'ai extrait un subset de 4 provinces et 13 communes présentant le cas que je décris (divergences de limites avec gap/overlap - la couche des provinces est la référence). Il n'y a pas de relation de clé étrangère entre les deux tables donc elles peuvent être chargées dans l'ordre souhaité. Les séquences SQL entrainent la création des tables dans le schéma public.

Merci pour le retour

Geoffroy

Elle est un peu pourrie ma méthode je pense

Ce qui marcherait bien c'est de travailler avec les arcs des communes et provinces, et reconstruire le set des arcs provinces pour l'exterieur et arc communes internes aux provinces, puis de reconstruire les pg avec st_polygonize.

... a suivre

Nicolas

Bonsoir Nico,

J'ai trouvé l'idée séduisante et avait débuté une implémentation sur cette base. Une première chose que j'ai faite est de ne retenir que les vertex "candidats" de la limite "commune", c'est-à-dire les vertex qui s'intersectent avec une commune voisine qui elle-même appartient à une province différente (en focalisant l'analyse sur la cod_mun='189110001', on dégage 70 points candidats qui doivent s'ajuster)

Ensuite, je n'ai pas eu de difficultés à trouver les vertex correspondants les plus proches sur la limite provinciale.

Évidemment, un des manquements de la méthode est qu'on ne tient pas compte du comportement même de la limite provinciale entre deux points "les plus proches" (cf. PJ). Des vertex sont donc manquants pour voir la limite de commune s'adapter à la limite de province.

En étudiant le cas, je me suis dit qu'il y avait moyen de récupérer tous les vertex du tronçon de limite provinciale situés entre le premier et le dernier point "les plus proche sur la limite provinciale", ces deux points étant connus par le traitement KNN. Cet ensemble de vertex viendrait remplacer le tronçon de limite communale situé entre le premier et dernier point "les plus proches sur la limite communale". Ensuite on reconstruit le polygone.

J'en étais là quand j'ai vu ton second post, que je vais lire attentivement. Je n'ai donc pas implémenté la recherche de tous les vertex provinciaux, mais en tout cas cela pourrait théoriquement tenir la route.

Je place ici l'implémentation SQL permettant de retrouver le vertex provincial le plus proche pour chacun des 70 vertex communaux candidats. Le code n'est pas optimisé du tout et des trucs ne servent à rien mais je le mets tel qu'il est en construction (pour si ça intéresse qq)

WITH 
com_target as (
    select cod_mun, mun_cod_prov, nombre_com, geom from dts_commune where cod_mun='189110001'
),
com_target_pt as (
    select cod_mun, mun_cod_prov, nombre_com, (ST_DumpPoints(geom)).geom as geom_com_pt from com_target 
),
com_neigh as (
    select cod_mun, mun_cod_prov, nombre_com, geom as geom_com_neigh from  dts_commune
    where st_intersects(geom, (select geom from com_target)) and mun_cod_prov<> (select mun_cod_prov from com_target)
),
-- Vertex candidats de la commune = vertex qui s'interesectent avec une commune voisine qui appartiene à une province différente
com_target_pt_cand as (
 SELECT row_number() OVER (ORDER BY com_target_pt.geom_com_pt) AS gid, com_target_pt.geom_com_pt
   FROM com_target_pt, com_neigh
  WHERE st_intersects(com_target_pt.geom_com_pt, com_neigh.geom_com_neigh)
),
prov_target as (
    select * from dts_province where cod_prov = '18911'
),
prov_target_pt as (
    select cod_prov, nom_prov, (ST_DumpPoints(geom)).geom as geom_prov_pt from dts_province where cod_prov = '18911'
),
prov_neigh as (
    select * from dts_province where st_intersects(geom, (select geom from dts_province where cod_prov='18911')) and cod_prov <> '18911'
)
select gid,  geom_com_pt, geom_prov_pt, dist
from com_target_pt_cand 
CROSS JOIN LATERAL 
  (SELECT
      cod_prov, geom_prov_pt,
      ST_Distance(geom_prov_pt, geom_com_pt) as dist
      FROM prov_target_pt
      ORDER BY geom_com_pt <-> geom_prov_pt
     LIMIT 1
   ) AS closest_vert_prov

Je garde l'idée de la sélection des vertex au frais et je vais étudier le travail réalisé à partir des arcs des deux couches. En tout cas merci beaucoup pour ce retour très détaillé, vraiment super. Je passe à l'action et ferai un retour

Geoffroy

Bonsoir,

J'ai reproduit la méthode d'arcs décrites par Nicolas. Très clair et avec la garantie de respecter la géométrie de la couche référence.

L'enjeu est de bien sélectionner les arcs internes. Pour généraliser la méthodologie à l'ensemble du dataset, je partirai sur un critère de sélection du style "arc communal interne = arc touchant le contour communal en maximum deux points. En effet je crains que la méthode basée sur la distance moyenne d'un arc communal à un arc provincial ne sélectionne pas tous les tronçon du dataset global. Je vais tester cela.

Par contre, j'ai été attiré par un autre élément. Lors de la reconstruction des polygones à partir de la multiligne (via ST_Polygonize), je m'aperçois qu'une limite n'est pas prise en compte et donc qu'un polygone plus important est créé, au lieu d'être séparé en deux (cf. PJ). Pourtant, ladite limite est bien dans la polyligne retournée en "tmp3" et topologiquement tout semble correct.

Dans le screenshot proposé par Nicolas, on retrouve bien les deux polygones.

J'ai examiné ce tronçon pour comprendre ce qu'il se passait mais sans vraiment trouver une explication, si ce n'est que la séquence des vertex ferait en sorte que la polyligne n'est pas considérée comme fermée.

@Nicolas, si vous avez l'occasion, pouvez-vous confirmer que chez vous, le polygone que j'identifie en rouge dans le screenshot est bien coupé en deux? Cela me permettrait de voir si il s'agit d'une limitation

Merci

Geoffroy

J'avais justement un PG13 / PostGIS 3.1.1 – Geos 3.8.1 sous la main pour des tests.

Après y avoir passé le WE, on observe clairement des différences de comportement entre les deux versions de PostGIS. Sur la v2.4, le résultat du st_polygonize ne renvoie pas un polygone et sur le v3.1, ce sont deux polygones qui sont absents (le même que sur la v2.4 + un petit autre situé juste à sa droite).

Je partage en PJ le dataset des arcs après snapping sur la v Postgis 3.1. Les polygones "à problème" reposent sur les arcs internes idarc = 18 et idarc=29

Si on détermine la distance de ces arcs à l'arc provincial le plus proche à partir des nœuds début-fin, le résultat indique 0m pour les deux arcs

with tmp as (
select idarc, st_startpoint(geom) as startpt_geom, 
    st_endpoint(geom) as endpt_geom
    from dts_arc_snap
    )
     select p.idarc, t.dist_start_prov,dist_end_prov
    from tmp p
             cross join lateral (
        select a.geom <-> p.startpt_geom as dist_start_prov                
        from dts_prov_arcs a
        order by a.geom <-> p.startpt_geom
        limit 1)t
            cross join lateral (
        select a.geom <-> p.endpt_geom as dist_end_prov                
        from dts_prov_arcs a
        order by a.geom <-> p.endpt_geom
        limit 1)s
        where idarc in (18,29)

Par contre l'exécution de ST_Touches renvoie false sauf pour le point d'arrivée de l'arc 18. Mais c'est également valable pour d'autres arcs qui se ferment correctement

select *, st_touches(a.geom, b.geom) from dts_arc_snap a,  dts_prov_arcs b where idarc in(18,29)

Et donc les deux polygones sont absents

        
with
 tmp as (
    select st_union(geom) as geom
    from dts_arc_snap
) select (st_dump(st_polygonize(geom))).geom
from tmp;

J'ai prolongé un des arcs pour le faire s'intersecter avec les limites provinciales (aux pts start/end). Le polygone est alors correctement construit.

Je vais continuer à investiguer  afin d'améliorer la qualité du snapping avant d'envisager de prolonger les arcs artificiellement.

Geoffroy

Sur cette version:

POSTGIS="3.0.3 6660953" [EXTENSION] PGSQL="120" GEOS="3.8.1-CAPI-1.13.3" SFCGAL="1.3.9" PROJ="7.2.0" GDAL="GDAL 3.2.0, released 2020/10/26" LIBXML="2.9.4" LIBJSON="0.15" LIBPROTOBUF="1.3.3" WAGYU="0.4.3 (Internal)" TOPOLOGY RASTER

J'obtiens le meme résultat que vous (pg manquants)
ca semble lié à la version de GEOS: soit vous mettez à jour, soit vous prolongez artificiellement les extrémités des arcs internes (ex ici: https://lists.osgeo.org/pipermail/postg … 6831.html)

(Mon conseil: mettez à jour: la derniere version postgis/geos est vraiment chouette !)

Nicolas

De rien

Vous pouvez aussi rapatrier les données depuis la prod sur votre serveur local, faire le traitement, et renvoyer le tout sur la prod (pg_dump ... | psql ...)

Nico