#1 Thu 17 October 2013 17:40
- ppluvinet
- Participant assidu
- Lieu: LYON
- Date d'inscription: 6 Aug 2007
- Messages: 609
POSTGIS2 : plus proches voisins un peu spécial
Bonjour,
J'ai une couche de polygones de type 'occupation du sol' où tous les polygones sont jointifs mais ne se chevauchent pas.
Je souhaiterais savoir comment lister les plus proches voisins avec la condition que si deux polygones ont juste 1 vertex commun ils ne soient pas des proches voisins.
J'ai testé la requête suivante. Celle-ci fonctionne mais je me demandais s'il n'y a pas quelque chose de + efficace :
Code:
select a.gid, b.gid from matable a, matable b where a.gid != b.gid and st_intersects(a.geom,b.geom) and st_geometrytype(st_intersection(a.geom,b.geom)) != 'POINT'
Deuxième question ou variante :
Je veux connaitre pour chaque polygone, le voisin ayant le maximum de contact ? Avez-vous une méthode assez efficace sans passer par le calcul d'intersection puis la selection du maximum?
Pascal PLUVINET
Hors ligne
#2 Thu 17 October 2013 19:47
- Nicolas Ribot
- Membre
- Lieu: Toulouse
- Date d'inscription: 9 Sep 2005
- Messages: 1538
Re: POSTGIS2 : plus proches voisins un peu spécial
Bonsoir,
Pour la premiere question, vous pouvez utiliser la fonction st_relate avec une matrice précisant le type de relation que vous voulez entre vos objets. C'est un peu abscons de prime abord, mais c'est tres puissant: les autres fonctions utilisent ceci en interne (st_intersects, st_touches, st_overlaps, etc):
http://postgis.net/docs/manual-2.1/ST_Relate.html
Dans votre cas, il faut exprimer le fait que la relation entre deux exterieurs de polygones est de type linéraire (1-dimension).
Il faut extraire les contours des polygones puis tester ces linestrings.
La deuxième partie de la requete teste si la dimension de l'intersection des intérieurs de deux lignes est de type lineaire (dimension 1)
-- plus proche voisin en se basant sur le contour:
-- extraction du contour des polygones avec prise en compte des trous éventuels.
-- Si les polygones sont simples, st_exteriorRing suffit pour extraire le contour.
with exterior as (
select gid, nom_comm,
st_exteriorRing((st_dumpRings((st_dump(geom)).geom)).geom) as geom
from communes
where code_dept = '31'
)
select a.gid, a.nom_comm, b.gid, b.nom_comm
from exterior a, exterior b
where a.gid <> b.gid
and a.geom && b.geom
and st_relate(a.geom, b.geom, '1********');
La meilleure perf sera obtenue avec une table intermédiaire contenant les contours et sur laquelle un index spatial sera défini (partie with exterior...)
Pour la deuxième partie, je ne vois pas comment calculer le max sans calculer la longueur de l'intersection. Par contre, en extrayant les contours des polygones, ca devrait etre plus rapide que l'intersection entre les polygones.
Nicolas
Hors ligne
#3 Fri 18 October 2013 07:10
- ChristopheV
- Membre
- Lieu: Ajaccio
- Date d'inscription: 7 Sep 2005
- Messages: 3175
- Site web
Re: POSTGIS2 : plus proches voisins un peu spécial
Bonjour
Si vous intégrez vos polygones dans le schéma topology tout devient plus simple.
en utilisant une commande équivalente à
Code:
SELECT topology.TopoGeo_AddPolygon('topo',st_geomfromwkb(:wkb, SRID),0);
Commande qui intègre le polygone (:wkb) dans le layer "topo" créé auparavant par
Code:
SELECT topology.AddTopoGeometryColumn('topo','SchemaName','my_table','the_topo','POLYGON');
Dans ce cadre votre demande est de connaître la face commune à un arc de votre polygone et les plus grand c'est l'arc qui a la plus grande longueur
Christophe
L'avantage d'être une île c'est d'être une terre topologiquement close
Hors ligne
#4 Fri 18 October 2013 11:06
- Nicolas Ribot
- Membre
- Lieu: Toulouse
- Date d'inscription: 9 Sep 2005
- Messages: 1538
Re: POSTGIS2 : plus proches voisins un peu spécial
Bonjour,
Effectivement, le module topology répond bien à ce problème.
Je n'ai pas testé sur la 2.1, mais sur la 2.0, la création d'une topologie prenait beaucoup de temps (plus de 30h sur les communes de France) et n'etait pas viable pour les gros volumes de données.
Nicolas
Hors ligne
#5 Fri 18 October 2013 12:15
- Nicolas Ribot
- Membre
- Lieu: Toulouse
- Date d'inscription: 9 Sep 2005
- Messages: 1538
Re: POSTGIS2 : plus proches voisins un peu spécial
Pascal, une précision concernant ta requete:
st_geometryType renvoie des type SQL/MM: ST_Point, ST_Linestring, etc.
La fonction geometryType (sans le prefixe st_) renvoie, elle, les types OGC: POINT, LINESTRING, etc.
Nicolas
Hors ligne
#6 Fri 18 October 2013 14:11
- ppluvinet
- Participant assidu
- Lieu: LYON
- Date d'inscription: 6 Aug 2007
- Messages: 609
Re: POSTGIS2 : plus proches voisins un peu spécial
Bonjour,
Merci pour vos messages. Je suis en train de faire quelques essais notamment avec St_relate, c'est très puissant ! Merci Nicolas pour ta dernière remarque mais je suis en train d'abandonner cette solution avec st_geometytype.
Les tests que j'ai fait récemment avec l'extension "topology" m'ont un peu dérouté car peu efficace sur des gros volumes de données.
C'est dommage, car semble très puissant et peu être très utile notamment pour la simplification toplogique de polygone.
Je suis aussi sur la 2.0. Pas encore testé également la 2.1.
Je vous fais un retour ces prochains jours...
Dernière modification par ppluvinet (Fri 18 October 2013 14:27)
Pascal PLUVINET
Hors ligne
#7 Fri 18 October 2013 14:13
- ChristopheV
- Membre
- Lieu: Ajaccio
- Date d'inscription: 7 Sep 2005
- Messages: 3175
- Site web
Re: POSTGIS2 : plus proches voisins un peu spécial
Bonjour,
@Nicolas
J'intègre mon cadastre Edigéo en topologie, effectivement c'est un peu long si on travaille "brut". En fait il faut faire un vaccum full régulièrement ça accélère les choses. Pour une vingtaine de lots soit environ 6000 polygones (c'est une moyenne) il me faut une heure. J'en suis à 452 971 polygones intégrés et il n'y a aucun soucis.
Et comme d'habitude si l'intégration est un peu longue en revanche les temps de réponse aux requêtes sont très, mais très rapides.
A+
Christophe
L'avantage d'être une île c'est d'être une terre topologiquement close
Hors ligne
#8 Fri 18 October 2013 15:04
- Nicolas Ribot
- Membre
- Lieu: Toulouse
- Date d'inscription: 9 Sep 2005
- Messages: 1538
Re: POSTGIS2 : plus proches voisins un peu spécial
Oui ca fait une belle volumétrie !
Du coup, tu conseilles de créer une topologie par petits incréments, plutot que d'un coup ? ca va plus vite ?
(je teste en créant la topo la facon suivante:)
select topology.ST_CreateTopoGeo('topo1',ST_Collect(geom))
from communes;
Nicolas
Hors ligne
#9 Fri 18 October 2013 16:25
- ChristopheV
- Membre
- Lieu: Ajaccio
- Date d'inscription: 7 Sep 2005
- Messages: 3175
- Site web
Re: POSTGIS2 : plus proches voisins un peu spécial
Bonjour,
C'est pas ST_CreateTopoGeo mais CreateTopoGeom
En fait j'utilise la fonction citée plus haut (TopoGeo_AddPolygon) polygone par polygone, cela construit directement, les nœuds, arcs et faces. Et ça renseigne tous les champs. Cela te renvoie l'ID de la face crée et tu fais le lien avec ton objet métier comme ceci
Code:
INSERT INTO schemaname.my_table (the_topo) VALUE (topology.CreateTopoGeom('topo',3,layerid, {{idface,3}}')) RETURNING idobjet;
Il faut juste faire attention au nombre de faces crées (si tu as une intersection entre deux polygones il y a au moins deux faces crées et ton objet métier faisant référence au polygone original comportera deux faces au moins, bon il suffit de faire un topoelement qui va bien du style
Code:
{{idface1,3},{idface2,3},...}
Voilà l'idée c'est de dissocier la création de la topologie et celle de l'objet métier possédant une topologie.
Tous les n polygones tu fais un vaccum full (tu vas voir qu'il y a suppression de ligne sur topo.nœud, arc et face et autre topo.relation) et ça roule impec.
Christophe
L'avantage d'être une île c'est d'être une terre topologiquement close
Hors ligne
#10 Mon 21 October 2013 12:36
- Nicolas Ribot
- Membre
- Lieu: Toulouse
- Date d'inscription: 9 Sep 2005
- Messages: 1538
Re: POSTGIS2 : plus proches voisins un peu spécial
Bonjour Christophe,
Merci pour ces explications !
Nicolas
Hors ligne
#11 Wed 12 March 2014 08:44
- ChristopheV
- Membre
- Lieu: Ajaccio
- Date d'inscription: 7 Sep 2005
- Messages: 3175
- Site web
Re: POSTGIS2 : plus proches voisins un peu spécial
Bonjour,
Je remonte ce sujet pour signaler (entre autre à Nicolas) que l'une des lenteur de l'intégration topologique vient d'un oubli dans la fonction st_addfacesplit, le ticket est géré par Sandro Santilli. Pour ceux que ça intéresse et que l'anglais rebute ci dessous le code de la fonction modifiée.
J'ajoute qu'après de nombreux test il apparaît clairement qu'une façon d’accélérer les choses est de lancer plusieurs instances d'intégration sur le modèle cité plus haut.
Code:
CREATE OR REPLACE FUNCTION topology._st_addfacesplit(atopology character varying, anedge integer, oface integer, mbr_only boolean) RETURNS integer AS $BODY$ DECLARE fan RECORD; newface INTEGER; sql TEXT; isccw BOOLEAN; ishole BOOLEAN; BEGIN IF oface = 0 AND mbr_only THEN RETURN NULL; END IF; SELECT null::int[] as newring_edges, null::geometry as shell INTO fan; SELECT array_agg(edge) FROM topology.getringedges(atopology, anedge) INTO STRICT fan.newring_edges; -- You can't get to the other side of an edge forming a ring IF fan.newring_edges @> ARRAY[-anedge] THEN RETURN 0; END IF; sql := 'WITH ids as ( select row_number() over () as seq, edge from unnest(' || quote_literal(fan.newring_edges::text) || '::int[] ) u(edge) ), edges AS ( select CASE WHEN i.edge < 0 THEN ST_Reverse(e.geom) ELSE e.geom END as g FROM ids i left join ' || quote_ident(atopology) || '.edge_data e ON(e.edge_id = abs(i.edge)) ORDER BY seq) SELECT ST_MakePolygon(ST_MakeLine(g.g)) FROM edges g;'; EXECUTE sql INTO fan.shell; isccw := NOT ST_OrderingEquals(fan.shell, ST_ForceRHR(fan.shell)); IF oface = 0 THEN IF NOT isccw THEN RETURN NULL; END IF; END IF; IF mbr_only AND oface != 0 THEN -- Update old face mbr (nothing to do if we're opening an hole) IF isccw THEN -- { sql := 'UPDATE ' || quote_ident(atopology) || '.face SET mbr = ' || quote_literal(ST_Envelope(fan.shell)::text) || '::geometry WHERE face_id = ' || oface; EXECUTE sql; END IF; -- } RETURN NULL; END IF; IF oface != 0 AND NOT isccw THEN -- { -- Face created an hole in an outer face sql := 'INSERT INTO ' || quote_ident(atopology) || '.face(mbr) SELECT mbr FROM ' || quote_ident(atopology) || '.face WHERE face_id = ' || oface || ' RETURNING face_id'; ELSE sql := 'INSERT INTO ' || quote_ident(atopology) || '.face(mbr) VALUES (' || quote_literal(ST_Envelope(fan.shell)::text) || '::geometry) RETURNING face_id'; END IF; -- } -- Insert new face EXECUTE sql INTO STRICT newface; -- Update forward edges sql := 'UPDATE ' || quote_ident(atopology) || '.edge_data SET left_face = ' || newface || ' WHERE left_face = ' || oface || ' AND edge_id = ANY (' || quote_literal(array( select +(x) from unnest(fan.newring_edges) u(x) )::text) || ')'; EXECUTE sql; -- Update backward edges sql := 'UPDATE ' || quote_ident(atopology) || '.edge_data SET right_face = ' || newface || ' WHERE right_face = ' || oface || ' AND edge_id = ANY (' || quote_literal(array( select -(x) from unnest(fan.newring_edges) u(x) )::text) || ')'; EXECUTE sql; IF oface != 0 AND NOT isccw THEN -- { -- face shrinked, must update all non-contained edges and nodes ishole := true; ELSE ishole := false; END IF; -- } -- Update edges bounding the old face sql := 'UPDATE ' || quote_ident(atopology) || '.edge_data SET left_face = CASE WHEN left_face = ' || oface || ' THEN ' || newface || ' ELSE left_face END, right_face = CASE WHEN right_face = ' || oface || ' THEN ' || newface || ' ELSE right_face END WHERE ( left_face = ' || oface || ' OR right_face = ' || oface || ') AND NOT edge_id = ANY (' || quote_literal( array( select abs(x) from unnest(fan.newring_edges) u(x) )::text ) || ') AND '; IF ishole THEN sql := sql || 'NOT '; END IF; sql := sql || '( ' || quote_literal(fan.shell::text) --********************************** rajout && geom || ' && geom AND _ST_Contains(' || quote_literal(fan.shell::text) -- We only need to check a single point, but must not be an endpoint || '::geometry, ST_Line_Interpolate_Point(geom, 0.2)))'; EXECUTE sql; -- Update isolated nodes in new new face sql := 'UPDATE ' || quote_ident(atopology) || '.node SET containing_face = ' || newface || ' WHERE containing_face = ' || oface || ' AND '; IF ishole THEN sql := sql || 'NOT '; END IF; sql := sql || 'ST_Contains(' || quote_literal(fan.shell::text) || '::geometry, geom)'; EXECUTE sql; RETURN newface; END $BODY$ LANGUAGE plpgsql VOLATILE COST 100; ALTER FUNCTION topology._st_addfacesplit(character varying, integer, integer, boolean) OWNER TO postgres; GRANT EXECUTE ON FUNCTION topology._st_addfacesplit(character varying, integer, integer, boolean) TO postgres;
Christophe
L'avantage d'être une île c'est d'être une terre topologiquement close
Hors ligne
#12 Wed 12 March 2014 10:25
- Nicolas Ribot
- Membre
- Lieu: Toulouse
- Date d'inscription: 9 Sep 2005
- Messages: 1538
Re: POSTGIS2 : plus proches voisins un peu spécial
Bonjour Christophe,
Merci pour cette info.
Nicolas
Hors ligne