Le portail francophone de la géomatique

Bonjour,

Vaste problème que la correction topologique de couches surfaciques.
Postgis n'est pas vraiment fait pour ca, et n'a pas d'outils magiques pour corriger ces couches.
La topologie est une solution partielle en ce sens que la construction d'une topo plante souvent qd les geom en entrées ne sont pas valides ou meme pas clean topologiquement.
Il y a des méthodes pour créer la topologie dans une fonction plpgsql, puis attraper les erreurs et de tenter de corriger la geom qui pose probleme avant de la réinjecter dans la topo.
Le pb de la topo est que c'est super lent de la créer.

On peut aussi faire par fonctions spatiales, en regardant a quels polygones sont collés les petits polygones à nettoyer, et en les unissant à ces polygones.
Par ex pour les petits trous dans les polygones, il est facile de les unir au polygone qui les contient.

Perso, pour faire ca sur des couches importantes de parcelles, j'ai laché postgis (sniff ) et suis passé par GRASS 7: il a un module pour écrire des topologies postgis directement.

Le process:
charger la couche dans grass
nettoyer la topo (clean, autres)
charger le résultat dans une topologie postgis
reconstruire les objets geo a partir de la topo.

Comme je suis une tache avec GRASS, j'ai jamais réussi l'opération d'associer certains des petits polygones aux voisins, ce que sait faire grass.
Au final, ces petits pg giclent lors de l'export en topo postgis, ce qui allait bien pour mon traitement. Par contre, ca bouche les petits trous et les jointures sont bonnes => couche clean au niveau topologique

Nico

Bonjour,

Pour la topologie je partage l'avis de Nicolat. Outre le fait que c'est "lent" à créer (il y a des astuces pour optimiser), au niveau d'une éventuelle correction la cela devient de la science fiction. J'ai tenté de supprimer automatiquement des arcs communs à des couples petites et grosses surfaces, un ça va mais plusieurs c'est la cata. En fait il faut faire un vaccum sur la base à chaque arc supprimé si on veut avoir qq chose de propre.
J'ai laissé tombé.
Par contre étant moi aussi un GRASS NULL je serai bien content d'étudier la piste .... des docs, tutos ?

Bonjour.
J'ai eu un problème similaire sur une couche d'occupation du sol. Quand j'ai découpé cette couche avec les limites communales pour affecter un code insee, je me suis retrouvé avec des micros polygones. J'ai donc fait une fonction qui fusionnait ces micro polygone avec le polygone "père" mitoyen (qui a le même identifiant). C'est une boucle qui nettoie les polygones <1m².
Certes le résultat en sortie ne colle pas pil poil aux limites communales mais ça évacue un certains nombre de micro entités.
Ca donne ça : (la couche de départ se nomme ocs)

CREATE OR REPLACE FUNCTION public.cagv_clean_ocs()
  RETURNS character varying AS
$BODY$DECLARE
--detail_index text;
inf_zero record;
inf_max record;
id_min integer;
ogc_fid_max integer;
nb_enregistrements integer;
lageom_union geometry;

i integer;
j integer;

BEGIN

i=0;
nb_enregistrements = 0;

FOR inf_zero IN (select ogc_fid, wkb_geometry, id from ocs where surf_m2 < 1)
LOOP
    j=1;
    i = i+1;
    id_min = inf_zero.id;
    nb_enregistrements =  max(rank) from (select ogc_fid, surf_m2, rank()  over (partition by id order by surf_m2 desc) from ocs where id = inf_zero.id and surf_m2 >= 1 )foo;
    
    FOR inf_max IN (select ogc_fid, wkb_geometry, surf_m2, rank()  over (partition by id order by surf_m2 desc) from ocs where id = inf_zero.id and surf_m2 >= 1 )
    LOOP
        lageom_union = null;        
        IF j <= nb_enregistrements and ST_TOUCHES(inf_zero.wkb_geometry, inf_max.wkb_geometry) THEN
            lageom_union = ST_UNION (inf_zero.wkb_geometry, inf_max.wkb_geometry) from ocs where ogc_fid = inf_max.ogc_fid;
            IF geometrytype(lageom_union) = 'POLYGON' THEN        
                update ocs set wkb_geometry = lageom_union, surf_m2 =st_area(lageom_union) where ogc_fid = inf_max.ogc_fid;
                delete from ocs where ogc_fid = inf_zero.ogc_fid;
                EXIT;
            ELSE
                j = j+1;
            END IF;
        ELSE
            j = j+1;
        END IF;    
    END LOOP;

    

--raise notice 'le nombre: %',i;
END LOOP;
return     'ok'; 
END;$BODY$
  LANGUAGE plpgsql VOLATILE SECURITY DEFINER
  COST 100;
ALTER FUNCTION public.cagv_clean_ocs()
  OWNER TO postgres;
COMMENT ON FUNCTION public.cagv_clean_ocs() IS 'Fonction de nettoyage des micro polygones issus du st_intersections pour découper selon les limites admin et affecter le code INSEE';

C'est certainement améliorable mais ça peut être une base de départ pour vous
Cordialement
JP

Bonjour,

A mon tour de filer du code

Pour cleaner des parcelles cadastrales par ex et générer une topologie clean:

Un shell qui execute du SQL postgis et du grass:

...
$PSQL $DBCON -c "SELECT topology.DropTopology('topo_grass')" -d $DBNAME
# lancement de grass en mode job batch, avec notre liste de commandes:
export GRASS_BATCH_JOB="${TMP_DIR}/clean_by_grass.bash"
echo "grass job: ${GRASS_BATCH_JOB}"
$GRASS $LOCATION

unset GRASS_BATCH_JOB
...

Un shell qui contient les commandes GRASS: clean_by_grass.bash:

#!/bin/bash
# nettoyage de la table parcelle et passage en topologie
########################################################################################################################
# Var utilisées
G_PG_CON="PG:dbname=${DBNAME} host=${DBHOST} port=${DBPORT} user=${DBUSER}"

echo "$TMP_TS: chargement de la table parcelle dans GRASS: supression de la couche grass parcelle ..."
g.remove -f type=vector name=parcelle@paris

v.in.ogr --overwrite snap=0.001 input="${G_PG_CON}" layer=parcenr.parcelle output=parcelle geometry=geom

TMP_TS=$(date +"%d-%m-%Y %H:%M:%S")
echo "$TMP_TS: chargement de la table parcelle dans GRASS: DONE"

# topo
TMP_TS=$(date +"%d-%m-%Y %H:%M:%S")
echo "$TMP_TS: construction de la topo postgis depuis grass___"

v.out.postgis --overwrite -l -2 \
    input=parcelle output="${G_PG_CON}" \
    output_layer=parcenr.parcelle_clean \
    options="SRID=2154,TOPOSCHEMA_NAME=topo_grass,TOPO_TOLERANCE=0.001,SPATIAL_INDEX=NO"

TMP_TS=$(date +"%d-%m-%Y %H:%M:%S")
echo "$TMP_TS: construction de la topo postgis depuis grass: DONE"

TMP_TS=$(date +"%d-%m-%Y %H:%M:%S")
echo "$TMP_TS: chargement parcelle clean avec geometrie depuis grass vers pg: DONE"

J'ai noté des pb de relations dans la topologie postgis écrite par GRASS, j'ai pas creusé, mais un update permet de régler cela:

-- grass ne produit pas de relation: on met a jour a partir de parcelle_clean pour
-- pouvoir exploiter la topo
insert into topo_grass.relation
  select (p.topo).id as topogeo_id, (p.topo).layer_id, f.face_id as element_id, (p.topo).type
    from parcenr.parcelle_clean p join topo_grass.face f on (p.topo).id = f.face_id;

Ensuite, et c'est là que c'est "tout pourri", je ne garde que le plus grand polygone parmi les polygones de meme id retournés par grass: les intersections entre parcelles (qui ne devraient pas exister) ont généré des petits polygones, mais attachés à une deux deux parcelles.
GRASS sait faire ca, et c'est même son but, c'est juste que j'ai compris/réussi à lancer les bonnes commandes (merger par catégorie par ex):

create table parcenr.parcelle as
  with tmp as (
      SELECT
        p.fid,
        p.id,
        p.idpar,
        p.code_insee,
        p.zonage,
        p.topo,
        topology.st_getFaceGeometry(topology.getTopologyName((p.topo).topology_id), (p.topo).id) AS geom
      FROM parcenr.parcelle_clean p
  ), tmp1 as (
  select t.*, row_number()
      OVER (PARTITION BY t.id
        ORDER BY st_area(t.geom) DESC) as rn
    from tmp t
  ) select t.id, t.idpar, t.code_insee, t.zonage, t.topo, t.geom
    from tmp1 t
  where rn = 1;

L'avantage est d'avoir au final une topologie clean des parcelles (exploitation des arcs pour trouver l'univers, les routes, les faces de batiments, etc...)
ET une table parcelle en objet qui est clean aussi, le tout chargé BEAUCOUP plus vite que dans Postgis: qq sec quand postgis met 20 min ou plus.

Nico

Rebonsoir,

J'ai hésité à créer un nouveau sujet. Mais il m'a semblé que le sujet était lié aux difficultés rencontrées ici.

Mon script a fait son petit bonhomme de chemin, mais n'est pas terminé.
Je vous posterai les étapes que j'ai choisies comme retour et si j'ai le temps je nous anonymiserai un jeu de données pour qu'on s'amuse à comparer les différentes techniques (hein Nicolas! ).

Je ne suis pas sur de partir sur des choses très catholiques, mais c'est comme ça qu'on apprend!
Je suis donc preneur de retours sur ma méthode (j'ai choisi de ne pas utiliser d'autres fonctions que celles de PostGIS).

Et comme un schéma vaut mieux qu'un long discours, vous trouverez en pièce jointe une capture d'écran.
Voici quelques étapes qui m'en font baver.

En jaune clair, la couche A à triturer.
En bleu clair, la couche B à ne pas intersecter et à laquelle se "recoller".
En rouge le résultat de l'étape 1
En bleu foncé (lignes): résultat de l'étape 2.

Étape 1:

--extraction des limites des polygones ST_boundary
--segmentation en nombreux petits segments sT_Segmentize
--extraction des points ST_DumpPoints
--snapping des points à masque0 ST_Snap
SELECT DISTINCT
    ST_Snap(
        (ST_DumpPoints(ST_Segmentize(ST_Boundary(t.geom),2))).geom
        , f.geom, 2) geom
    , t.id_pk
FROM saisie.saisie_e1 t, masque0 f

Premier point que je n'arrive pas à régler. J'aimerais aimanter tous les points de A aux polygones B (y compris aux arrêtes des polygones et pas seulement aux sommets). Auriez-vous une piste de résolution?

Étape 2:

SELECT DISTINCT ST_ForceRHR(ST_MakeLine(geom)), id_pk
FROM etape1
GROUP BY 2

Je crée des lignes (ST_MakeLine) en vue de créer des polygones (ST_Polygonize), groupés selon les id_pk.
J'ai testé avec ou sans la fonction ST_ForceRHR(), je pensais qu'en forçant la "digitalisation main droite" j'arriverai à supprimer l'effet zig-zag que vous pouvez admirer sur ma PJ.
Mais cette fonction ne résout pas le pb, et ne change rien au résultat!
Deuxième point que je cherche à résoudre: faire comprendre à PostGIS qu'il doit créer des lignes en reliant les points les plus près entre eux, lorsqu'ils partagent le même id_pk. Auriez-vous là aussi des pistes de résolution?

Merci d'avance et bonne soirée.

Dernière modification par Mathieu Denat (Wed 14 February 2018 18:37)