utilisation de shapely : is_ring

Version imprimable

24/11/2011, 13h59
awalter1

utilisation de shapely : is_ring
Bonjour,
Je commence à utiliser le module shapely. Je ne comprends pas la notion de ring malgré le manuel. En particulier l'utilisation de la propriété is_ring me renvoie toujours False :
Code:

1 2 3 4 polygon=Polygon(((0,0),(0,4),(4,4),(4,0),(0,0))) print polygon.is_ring False
Merci
24/11/2011, 14h03
awalter1

De plus, y a t'il une importance ou pas de définir pour un polygone une suite de coordonnées de points avec coord. dernier point = coord. premier point. Dans les exemples, c'est le cas, mais si on ne le fait pas je ne vois pas l'impact !
Merci

Bonjour

is_ring est plus un attribut (@property je crois) à utliser sur les instances de LineString (arcs ouverts ou fermés) et LinearRing (arcs fermés). Encore que sur LinearRing, ça n'a guère d'intérêt, étant donné que is_ring est vrai pour un arc qui a les mêmes première et dernière coordonnées, ce qui est une caractéristique des instances de LinearRing.

Je n'avais même pas essayé de l'utiliser sur un Polygon. Et si cela levait une exception, ça ne me dérangerait pas non plus.

Par contre, tu peux le tester sur l'attribut exterior (périmètre extérieur) et chacun, s'il exsitent, des interiors (trous) : ça te retournera True. Mais ça a autant d'intérêt que de l'utiliser sur un LinearRing.

Pour les Polygon, lors de l'appel du constructeur, tu n'es pas obligé d'avoir les mêmes premier et dernier couples de coordonnées (c'est "rajouté", si besoin, automatiquement). C'est écrit quelque part dans la doc

doc constructeur de Polygon

Citation:

The Polygon constructor takes two positional parameters. The first is an ordered sequence of (x, y[, z]) point tuples and is treated exactly as in the LinearRing case.

et doc constructeur de LinearRing :

Citation:

The sequence may be explicitly closed by passing identical values in the first and last indices. Otherwise, the sequence will be implicitly closed by copying the first tuple to the last index.
Code:

1 2 3 4 5 6 7 8 9 10 11 12 >>> from shapely.geometry import Polygon >>> p1=Polygon(((0,0),(0,1),(1,1),(1,0))) >>> p2=Polygon(((0,0),(0,1),(1,1),(1,0),(0,0))) >>> p1.wkt 'POLYGON ((0.0000000000000000 0.0000000000000000, 0.0000000000000000 1.0000000000000000, 1.0000000000000000 1.0000000000000000, 1.0000000000000000 0.0000000000000000, 0.0000000000000000 0.0000000000000000))' >>> p2.wkt 'POLYGON ((0.0000000000000000 0.0000000000000000, 0.0000000000000000 1.0000000000000000, 1.0000000000000000 1.0000000000000000, 1.0000000000000000 0.0000000000000000, 0.0000000000000000 0.0000000000000000))' >>> >>> p1.equals(p2) True >>> p1.exterior.is_ring True

24/11/2011, 15h36
awalter1

Merci de ces infos.
A présent je veux faire quelques chose de relativement simple ; tester l'intersection de 2 polygones, pourtant je tourne en rond.
J'ai testé les methodes intersects(), intersection(), crosses(), mais je n'ai pas les résultats espérés. Je pense que cela vient du fait que certains polygones ont des frontières communes, ce qui est normal pour moi mais doit être vu comme des polygones disjoints : ce n'est pas la logique de shapely.
Comment puis je savoir si 2 polygones se recouvrent (sans tenir compte des segments communs) ?
Dans le cas ou je fais P1.intersection(P2) pour des polygones qui ont des segments communs, je peux tester si le résultat est une polyline et j'en déduirai qu'il n'y a pas overlap, mais ce n'est pas très efficace.
Merci
24/11/2011, 16h01
plxpy

Citation:

mais ce n'est pas très efficace.

pourquoi ? Tester la classe de l'objet retourné par la méthode intersection est une bonne idée. Ce que tu veux, c'est une intersection qui soit de la classe Polygon ou MultiPolygon (ou qui ait une partie qui soit dans ce cas là: attention aux résultats/intersections de type GeometryCollection !)

Et si c'est au niveau des performances que tu trouves que ce n'est pas très efficace, c'est peut-être que tu as de nombreux polygones et que tu testes tous les couples possibles, non ?

En fait, si tu veux déterminer les polygones qui se "recouvrent" (qui ont une intersection surfacique et non pas linéaire), le mieux est encore d'utiliser les prédicats intersects et touches.

Pour avoir une intersection surfacique, les polygones doivent s'intersecter mais ne pas se toucher :
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 >>> from shapely.geometry import Polygon >>> p1=Polygon(((0,0),(0,1),(1,1),(1,0))) >>> p2=Polygon(((1,0),(1,1),(2,1),(2,0))) >>> >>> p1.intersects(p2) True >>> p1.touches(p2) True >>> >>> p3=Polygon(((0.5,0),(0.5,1),(1.5,1),(1.5,0))) >>> p1.intersects(p3) True >>> p1.touches(p3) False
La "logique" Shapely est héritée de la bibliothèque C/C++ GEOS, portage de JTS (Java Topology Suite, voir les spécifications techniques), elle-même calquée sur les standards OGC. Elle peut parfois sembler déroutante au premier abord mais, en fait, c'est la façon de traiter les cas "finauds" (tordus ?) qu'on rencontre plus souvent qu'on ne le pense (et je sais de quoi je parle !)

25/11/2011, 09h42
awalter1
oui, en creusant un peu on trouve ce qui convient.
au final je fais
Code:

1 2 inter = poly1.intersection(poly2) if type(inter) == typ_polygon:

Tu traites surement beaucoup de cas mais ça ne marche pas avec ces deux polygones par exemple :
Code:

1 2 3 4 5 6 7 8 9 >>> from shapely.geometry import Polygon >>> >>> polygone_en_U = Polygon(((0,0),(0,5),(1,5),(1,1),(2,1),(2,5),(3,5),(3,0))) >>> polygone_horiz = Polygon(((0,3),(0,4),(4,4),(4,3))) >>> >>> inter = polygone_en_U.intersection(polygone_horiz) >>> >>> type(inter) <class 'shapely.geometry.multipolygon.MultiPolygon'>
Mais tout dépend des polygones que tu traites. Peut-être que ce cas de figure ne se présente (mais, perso, j'ai appris à me méfier, dans ce domaine, des cas de figure qui n'arrivent "jamais" ...)

25/11/2011, 09h59
awalter1
oui, en creusant un peu on trouve ce qui convient.
au final je fais
Code:

1 2 3 4 typ_polygon=type(poly1) inter = poly1.intersection(poly2) if type(inter) == typ_polygon: print "c'est bien une intersection de surface"
Maintenant, comment détecter un polygone dont au moins 2 segments se croisent (à moins de tester tous les couples de segments soi même). A priori, ce n'est pas détecter lors de la création du polygone.

Merci
25/11/2011, 10h01
awalter1

j'ai trouvé : on fait poly1.is_valid

Citation:

j'ai trouvé : on fait poly1.is_valid

effectivement, is_valid t'indique si une géométrie est ... valide ou pas.

Par contre, ce que tu fais ne marche pas :
Code:

1 2 3 4 5 6 7 8 9 10 11 12 >>> pu = polygone_en_U >>> ph = polygone_horiz >>> >>> pu.intersects(ph) and not pu.touches(ph) True >>> >>> typ_polygon=type(pu) >>> inter = pu.intersection(ph) >>> if type(inter) == typ_polygon: ... print "c'est bien une intersection de surface" ... >>>
et mes deux polygones sont tout ce qu'il y a de plus valides :
Code:

1 2 3 4 >>> pu.is_valid True >>> ph.is_valid True

25/11/2011, 13h25
awalter1
J'ai une intersection qui me renvoie une linestring et un polygon :
Code:

1 2 print ">>> ",poly1.intersection(poly2) >>> GEOMETRYCOLLECTION (LINESTRING (213551.0 19113.0 ....), POLYGON ((213878.0 19114.0 ....)))
Comment puis je extraire ces deux entités ?
Merci
25/11/2011, 13h43
plxpy

A partir du moment où tu as une "multi-géométrie" (MultiPoint, MultiLineString, MultiPolygon et, à fortiori, une GeometryCollection), les géométries sont accessibles via l'attribut geoms
28/11/2011, 09h21
awalter1
Bonjour,
Je vois que si je fais
Code:

1 2 p=Polygon([[0,0][0,5],[5,5][5,0]]) print p.geoms
j'ai une erreur 'Polygon' object has no attribute geoms. Je comprends alors que geoms ne s'applique que sur des collections de pts, lignes ou polygones (geoms pourrait renvoyer null). Mais alors, comment tester que l'on a pas une figure simple : 1 point, 1 ligne, 1 polygone ? A moins qu'il faille gérer cela par un try/except ?
Merci

PS: J'ai un doute de matin : vendredi soir mon jeu de test faisait apparaitre une intersection qui me donnait une GEOMETRYCOLLECTION (d'ou le mail précédent). Ce matin je relance mon application avec le même jeu de test, j'obtiens un POLYGON ! Shapely est il bien déterministe et fiable ?
28/11/2011, 10h28
plxpy

Citation:

Shapely est il bien déterministe et fiable ?

Sans aucune hésitation : oui !

Je l'utilise depuis plus de 3 ans dans un contexte professionnel sur des données géographiques "musclées" : des objets géographiques de tous types (ponctuel, linéaire, surfacique, collection) et parfois de formes particulièrement "torturées". J'ai développé de nombreuses applications avec et n'ai pas rencontré de problème particulier.

Shapely est basé sur la bibliothèque GEOS, utilisée dans de nombreux projets et outils OpenSource professonnel (PostGis, OpenJump (la version java) et j'en oublie). Si elle était instable ou peu fiable, il y a belle lurette qu'elle serait passée à la trappe.

Après, aucune application ne peut se vanter d'être "bug free" ! Mais si tu démarres dans l'utilisation de Shapely, il est plus raisonnable de chercher l'erreur ailleurs pour commencer.

Citation:

j'ai une erreur 'Polygon' object has no attribute geoms

Normal ! Si tu relis l'un de mes posts, tu verras que j'en ai parlé pour des multi-géométries : MultiPoint, MultiLineString, MultiPolygon et GeometryCollection. Une instance de Polygon (ou de Point ou de LineString ou de LinearRing) n'a pas d'attribut geoms

Citation:

... une intersection qui me donnait une GEOMETRYCOLLECTION (d'ou le mail précédent). Ce matin je relance mon application avec le même jeu de test, j'obtiens un POLYGON !

Désolé mais je n'y crois pas. Tu as forcément loupé quelque chose. Poste des données (celles qui posent problème) si tu veux que j'y jette un coup d'oeil

1 pièce(s) jointe(s)

En rapport avec ton post sur comment tester le type des objets.

polygons est un dictionnaire dont les clés sont les couleurs des polygones (voir image jointe). J'ai chargé au préalable les polygones et je teste les intersections entre le polygone vert (variable vert) et les autres. A la fin, je passe en revue les géométries de l'intersection "GeometryCollection" entre le vert et le marron

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
>>> from shapely.geometry import Point, LineString, Polygon, GeometryCollection
>>>
>>> for couleur in polygons:
...     inter = vert.intersection(polygons[couleur])
...     if inter.is_empty:
...             print "pas d intersection avec le polygone", couleur
...     elif isinstance(inter,Point):
...             print "intersection ponctuelle avec le polygone", couleur
...     elif isinstance(inter,LineString):
...             print "intersection lineaire avec le polygone", couleur
...     elif isinstance(inter,Polygon):
...             print "intersection surfacique avec le polygone", couleur
...     elif isinstance(inter,GeometryCollection):
...             print "intersection compliquee avec le polygone", couleur
...
intersection ponctuelle avec le polygone bleu
intersection surfacique avec le polygone orange
pas d intersection avec le polygone violet
intersection compliquee avec le polygone marron
intersection lineaire avec le polygone rouge
>>>
>>>>>> vert_marron = vert.intersection(polygons['marron'])
>>>
>>> for geom in vert_marron.geoms:
...     if isinstance(geom,Point):
...             print "un point"
...     elif isinstance(geom, LineString):
...             print "un arc"
...     elif isinstance(geom,Polygon):
...             print "un polygone"
...     else:
...             print "autre chose"
...
un point
un arc
un arc
un polygone
>>>

Attention : le code n'est pas complet. Avec d'autres polygones, on pourrait obtenir aussi des MultiPoint, MultiLineString ou MultiPolygon. Dans mon code, je n'ai pas géré ces possibilités. C'est un exemple.

28/11/2011, 12h24
plxpy

Ce serait ballot de ne pas profiter de l'occasion : un article sur Shapely sur le site portail SIG qui fait bien le tour de la question et qui est un très bon complément à la documentation officielle (la nouvelle et l'ancienne).

Pour démarrer, l'ancienne est peut_être la plus appropriée. Moins touffue, elle va à l'essentiel.

Merci pour votre aide, c'est clair.
Et vous aviez raison, j'ai du cafouiller dans mes tests. Je préfère savoir que Shapely est fiable, il est en plus performant, quant à son utilisation ça vient petit à petit, il y a en effet des situations que je n'imaginais pas (trop concentré su ma seul problématique).
Mon intérêt à Shapely vient du besoin de déterminer les "holes and overlaps" sur quelques centaines de volumes (polygones en x,y,z avec z constant). Cela était fait auparavant en C en exploitant les Regions X de la Xlib, mais je crois que j'atteins les limites. Avec Shapely, c'est plus simple et plus puissant.
La partie détection des "overlaps" est ok, je vais m'interresser à la détection des trous.
Merci encore

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
# l_polygon contains a list of polygons
#  l_polygon[i][0] = polygon name
#  l_polygon[i][1] = upper level
#  l_polygon[i][2] = lower level
#  l_polygon[i][3] = list of points coordinates
#
for l in range(1,len(l_level)):
       selected_level = l_level[l]
 
       for i in range(len(l_polygon)):
               if selected_level <= l_polygon[i][1] and selected_level > l_polygon[i][2]:
                       for j in range(i+1,len(l_polygon)):
                               if selected_level <= l_polygon[j][1] and selected_level > l_polygon[j][2]:
                                       inter = l_polygon[i][3].intersection(l_polygon[j][3])
                                       # if intersection
                                       if inter.is_empty == False:
                                               l_inter = []
                                               # check if collection or polygon, exclude point and linestring
                                               if isinstance(inter,Polygon):
                                                       l_inter.append(inter)
                                               elif isinstance(inter,GeometryCollection):
                                                       for geom in inter.geoms:
                                                               if isinstance(geom,Polygon) == True:
                                                                       l_inter.append(geom)
                                               for poly in l_inter:
                                                       print "Intersection between volumes <%s> and <%s> between level
%d and %d" % (l_polygon[i][0],l_polygon[j][0],min_level,selected_level)
                                                       # find the points name corresponding to the overlapping polygon
                                                       l_points_overlap = []
                                                       point_overlap = ""
                                                       for coord in list(geom.exterior.coords)[:-1]:
                                                               for pt in l_points:
                                                                       if (str(coord[0]) == "%s.0" % pt[1]) and (str(coord[1]) ==
"%s.0" % pt[2]):
                                                                               if pt[0] not in l_points_overlap:
                                                                                       point_overlap = pt[0]
                                                               # the point name is not found : get only its coordinates
                                                               if point_overlap == "":
                                                                       l_points_overlap.append((self.ApplicationConvertNumlatlong(coord[0],coord[1])))
                                                               else:
                                                                       l_points_overlap.append(point_overlap)
 
                                                       ch_points_overlap = ""
                                                       for pt in l_points_overlap:
                                                               if type(pt) == type(""):
                                                                       ch_points_overlap = "%s %s" % (ch_points_overlap,pt)
                                                               else:
                                                                       ch_points_overlap = "%s %s/%s" % (ch_points_overlap,pt[0],pt[1])
 
                                                       print "   intersection definition : %s" % (ch_points_overlap)

28/11/2011, 13h07
plxpy

Quelques petites remarques

Attention, Shapely c'est de la 2D. Avec des z constants, on peut peut-être s'y ramener mais faut le faire "à la main". Pour du vrai 3D, il y a par exemple CGAL. C'est codé en C/C++, il y a une API pour Python.

Dans le code (que je n'ai pas encore regardé en détail), le premier truc qui m'a frappé c'est l'utilisation de liste pour grouper les caractéristiques de chaque polygone d'où les [0] ... [3]. Un petit dictionnaire ferait bien l'affaire, avec des clés bien lisibles.

Si tu as beaucoup de polygones à examiner, l'utilisation d'un index spatial peut accélérer les choses. Il y en a un tout trouvé, c'est rtree et c'est aussi Sean Gillies (le "papa" de Shapely) qui l'a écrit.

Sinon, tes "levels" correspondent aux différents z ?
28/11/2011, 17h22
awalter1

oui, du 2D me suffit car le z est constant, et cette dimension corresponds en effet à mes "level".
Un dico rendrait en effet le pg plus lisible. J'en suis dans un premier temps à garantir la faisabilité.
Pour la rapidité, avec un jeu de données opérationnelles, l'exécution est tellement rapide par rapport à la version C antérieur que l'optimisation me semble superflu.
Cdlt

1 pièce(s) jointe(s)

Bonjour,
J'ai mis au point l'implémentation de la détection de trous et de recouvrement entre polygones. Cela marche bien lorsqu'il y a recouvrement d'un polygone sur un autre. Cependant lorsque le recouvrement se fait sur plusieurs polygones, un trou est détecté. Etrange.
Sur le cas fourni ci-dessous, un trou est détecté en :
hole= POLYGON ((214263.6511675126967020 17843.6196834876172943, 214699.0000000000000000 19093.0000000000000000, 214263.6511675126967020 17843.6196834876100183, 214263.6511675126967020 17843.6196834876172943))
Ce qui donne en lat/long un polygone avec 2 points identiques à cause de la précision insuffisante : ('5931N03', '00457E23'), 'AXA2', ('5931N03', '00457E23').
Je ne m'explique pas ce trou.
Si par ailleurs je rapproche AXA2 de P14, le triangle cerclé de rouge est détecté aussi en tant que trou.
Y a t'il une explication qui m'échappe ?
Merci

Mon example comporte 3 polygones:
DR1 = P14 P15 AXA2 AXA1 AXA6B AXA6 AXA6A AXA5 P13 P13A P10 P11 P14
AX1 = AXA1 AXA2 HD2 HD1 HD5 HD6 P16 P61 ZV2 ZV5 ZVB AXA5 AXA6A AXA6 AXA6B AXA1
DR1B = OF24 HD2 HD6 HD5 OF24
(Le point OF24 est normalement en HD1 pour qu'il n'y ai ni trou ni recouvrement).
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 # get the sectors envelope l_polygon = d_polygons.keys() upper,lower,envelop = d_polygons[ l_polygon[0] ] for l in range(len(l_level)): for i in range(1,len(l_polygon)): # add the current polygon upper,lower,points = d_polygons[ l_polygon[i] ] envelop = envelop.union(points) # get the rectangular volume that includes all sectors volumes envelop = envelop.envelope # check if all volumes equals the envelope d_holes = {} for l in range(0,len(l_level)): first = True for i in range(len(l_polygon)): i_v_upper,i_v_lower,i_v_points = d_polygons[ l_polygon[i] ] if l_level[l] <= i_v_upper and l_level[l] > i_v_lower: # add the current polygon if first == True: volumes = i_v_points first = False else: volumes = volumes.union(i_v_points) if first == False: if envelop.equals(volumes) == False: holes = envelop.difference(volumes) l_holes = [] # check if collection or polygon, exclude point and linestring if isinstance(holes,Polygon): l_holes.append(holes) elif isinstance(holes,MultiPolygon): for geom in holes.geoms: if isinstance(geom,Polygon) == True: l_holes.append(geom) for h in l_holes: l_pts = self.ApplicationGetVolumePoints(h,d_points) if d_holes.has_key(l_level[l]): l = d_holes[ l_level[l] ] l.append(l_pts) d_holes[ l_level[l] ] = l else: d_holes[ l_level[l] ] = [l_pts]

01/12/2011, 11h29
plxpy

Peux-tu envoyer un zip de 3 fichiers (1 par polygone) ? Utilise le format wkb comme ça (poly est une instance de Polygon) :

Code:

>>> open('poly.wkb','w').write(poly.wkb.encode('hex'))

Avec le schéma et seulement les points en texte, c'est l'horreur à regarder/débugger

1 pièce(s) jointe(s)

Voici 3 fichiers dans le zip.
En précisant mes tests j'ai mis en avant que l'ordre d'union de polygones ne donne pas le même résultat et c'est cela qui introduit un trou.
Dans mon zip, 1 fichier pour chaque polygone : Other, DR3, DR1B
Dans les traces suivantes, je fais l'union des polygones Other, DR1B, DR3 et j'ai un trou qui apparait.
Je fais ensuite l'union des polygones Other, DR3, DR1B : je n'ai pas de trou:

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
print "union volumes =",allvolumes
 
i_v_upper,i_v_lower,i_v_points = d_polygons[ 'DR1B' ]
volumes = allvolumes.union(i_v_points)		
i_v_upper,i_v_lower,i_v_points = d_polygons[ 'DR3' ]
volumes = volumes.union(i_v_points)
print "union volumes+DR1B+DR3 =",volumes
 
i_v_upper,i_v_lower,i_v_points = d_polygons[ 'DR3' ]
volumes = allvolumes.union(i_v_points)		
i_v_upper,i_v_lower,i_v_points = d_polygons[ 'DR1B' ]
volumes = volumes.union(i_v_points)
print "union volumes+DR3+DR1B =",volumes

et les traces :

Code:

1
2
3
union volumes = POLYGON ((211500.0000000000000000 57600.0000000000000000, 239400.0000000000000000 57600.0000000000000000, 239400.0000000000000000 40800.0000000000000000, 239400.0000000000000000 -7200.0000000000000000, 194400.0000000000000000 -7200.0000000000000000, 194400.0000000000000000 12000.0000000000000000, 194400.0000000000000000 57600.0000000000000000, 211500.0000000000000000 57600.0000000000000000), (214699.0000000000000000 19093.0000000000000000, 213878.0000000000000000 19114.0000000000000000, 213551.0000000000000000 19113.0000000000000000, 213457.0000000000000000 19120.0000000000000000, 212878.0000000000000000 20273.0000000000000000, 212096.0000000000000000 20298.0000000000000000, 212096.0000000000000000 19837.0000000000000000, 211856.0000000000000000 19837.0000000000000000, 211856.0000000000000000 20307.0000000000000000, 211463.0000000000000000 20323.0000000000000000, 210198.0000000000000000 21351.0000000000000000, 210271.0000000000000000 18840.0000000000000000, 211161.0000000000000000 17529.0000000000000000, 212102.0000000000000000 17020.0000000000000000, 213900.0000000000000000 16800.0000000000000000, 214699.0000000000000000 19093.0000000000000000))
union volumes+DR1B+DR3 = POLYGON ((211500.0000000000000000 57600.0000000000000000, 239400.0000000000000000 57600.0000000000000000, 239400.0000000000000000 40800.0000000000000000, 239400.0000000000000000 -7200.0000000000000000, 194400.0000000000000000 -7200.0000000000000000, 194400.0000000000000000 12000.0000000000000000, 194400.0000000000000000 57600.0000000000000000, 211500.0000000000000000 57600.0000000000000000), (214263.6511675126967020 17843.6196834876172943, 214263.6511675126967020 17843.6196834876100183, 214699.0000000000000000 19093.0000000000000000, 214263.6511675126967020 17843.6196834876172943))
union volumes+DR3+DR1B = POLYGON ((211500.0000000000000000 57600.0000000000000000, 239400.0000000000000000 57600.0000000000000000, 239400.0000000000000000 40800.0000000000000000, 239400.0000000000000000 -7200.0000000000000000, 194400.0000000000000000 -7200.0000000000000000, 194400.0000000000000000 12000.0000000000000000, 194400.0000000000000000 57600.0000000000000000, 211500.0000000000000000 57600.0000000000000000))

Dans le 2eme print, on voit que le polygone a un trou.
Cdlt

02/12/2011, 12h49
awalter1

4 pièce(s) jointe(s)

j'ajoute les vues des polygones correspondants: dans l'ordre DR1B, DR3, DR1B+DR3, enfin Other = la zone complete-DR3-DR1B
Cdlt
02/12/2011, 12h53
plxpy

je regarde tout ça (voir tes MP)

5 pièce(s) jointe(s)

Selon l'ordre dans lequel on utilise les polygones envoyés pour calculer l'union, un (micro) trou peut apparaître ou non dans le résultat (voir la localisation sur Fig0.png : la ligne noire). C'est théoriquement/mathématiquement choquant mais ici ça s'explique.

Tout vient du fait qu'il existe des segments qui s'intersectent sans qu'il y ait de points existants dans les données en entrée (voif Fig1.png - les petits carrés jaunes sont les points existants : il en "manque" au niveau du cercle rouge). Autrement dit, le graphe correspondant n'est pas planaire. Ensuite, c'est un problème de pur calcul numérique.

Il y a (visuellement) 3 segments concernés (4 en réalité) : 2 participant au périmètre du polygone bleu, un même segment (double) dans les polygones vert et rouge.

Quand on fait l'union entre le bleu et le rouge, les coordonnées des points d'intersection (points manquants) sont calculées dans un certain contexte géométrique et intégrées dans le résultat de l'union (voir Fig2.png)

Quand on fait l'union entre ce résultat et le dernier polygone vert, ces mêmes coordonnées sont calculées pour le segment du polygone vert dans un tout autre contexte : les deux segments superposés du premier résultat, chacun dans un trou différent. Ensuite, les problèmes inhérents au calcul avec des flottants fait le reste : du côté du polygone vert, on ne retrouve pas exactement les mêmes coordonnées d'où le micro-trou.

C'est un problème incontournable et lié à la méthode de calcul. J'ai d'ailleurs mené le même calcul d'union, à la "main", avec OpenJump qui utilise JTS (donc la version Java, sans code commun avec GEOS/Shapely) : j'obtiens le même micro trou. Shapely n'est pas la cause.

Parmi tous les tests que j'ai effectués, celui-ci me semble parfaitement montrer ce que j'avance. Je suis parti de tes polygones et je leur ai, en bouclant, degré par degré, sur un tour complet, appliqué une rotation changeant ainsi leurs coordonnées tout en gardant une géométrie/configuration relative constante mais se plaçant dans des conditions très variées, numériquement parlant, au niveau des coordonnées.

Résultat : zéro, un ou deux micros trous apparaissent selon la valeur de l'angle. Et, sur les 6 combinaisons possibles (ordre dans lequel on utilise les polygones), on a, selon les cas, 1, 2 ou 4 combinaisons qui provoquent ce(s) trou(s).
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 # -*- coding: utf-8 -*- import sys from itertools import permutations from shapely.wkb import loads from shapely.geometry import Polygon from math import sin, cos, radians from glob import glob def union_has_holes(polygones): errors = list() for o,ordre in enumerate(permutations((0,1,2))): # nouvelle liste avec les polygones dans l'ordre courant polygones_ordonnes = [ polygones[i] for i in ordre ] # calcul de l'union, 2 par 2, dans l'ordre courant union = reduce(lambda a,b:a.union(b),polygones_ordonnes) # si le resultat a un ou des trous, on memorise l'ordre if len(union.interiors): errors.append(ordre) return errors def tourne(coords,i): S, C = sin(radians(i)), cos(radians(i)) newcoords = map(lambda (x,y):(C*x-S*y,S*x+C*y),coords) return newcoords # pour compter les erreurs, configuration par configuration # (ordre des polygones dans le calcul) errors_dict = dict() # les polygones sont charges dans l'ordre DR1B.wkb (0), DR3.wkb (1) puis Other.wkb (2) print sorted(glob('*.wkb')) polygones = [ loads(open(name).read().decode('hex')) for name in sorted(glob('*.wkb')) ] # fichiers format WKT pour stocker les cas avec trous (KO_wkt) et sans trou (OK_wkt) OK_wkt = open('rotation_OK.wkt','w') KO_wkt = open('rotation_KO.wkt','w') # pour compter les differents cas with_holes = without_holes = 0 # coordonnées des différents périmètres (en dur, en connaissant les caractéris- # tiques des polygones et l'ordre du chargement ! EXT_0 = tuple(polygones[0].exterior.coords)[:-1] EXT_1 = tuple(polygones[1].exterior.coords)[:-1] EXT_2 = tuple(polygones[2].exterior.coords) INT_2 = tuple(polygones[2].interiors[0].coords)[:-1] for i in range(360): ext0 = tourne(EXT_0,i) ext1 = tourne(EXT_1,i) ext2 = tourne(EXT_2,i) int2 = tourne(INT_2,i) new_polygones = [Polygon(ext0),Polygon(ext1),Polygon(ext2,(int2,))] # "cast" en tuple pour servir de cle dans le dictionnaire errors = tuple(union_has_holes(new_polygones)) if errors: errors_dict[errors] = errors_dict.get(errors,0) + 1 which_wkt = KO_wkt else: which_wkt = OK_wkt # ecriture des 3 polygones dans le bon fichier for poly in new_polygones: print >> which_wkt, poly.wkt OK_wkt.close() KO_wkt.close() print "les combinaisons d'erreurs :" for k in errors_dict: print "%s : %d" % (k, errors_dict[k]) errors = sum(errors_dict.values()) print "soit %d unions avec trous, %d sans trou" % (errors,360-errors)
Code:

1 2 3 4 5 6 7 8 9 10 plx@sony:~/Bureau$ python rotations.py ['DR1B.wkb', 'DR3.wkb', 'Other.wkb'] les combinaisons d'erreurs : ((0, 1, 2), (0, 2, 1)) : 6 ((0, 2, 1), (2, 0, 1)) : 92 ((0, 1, 2), (0, 2, 1), (1, 0, 2), (2, 0, 1)) : 53 ((0, 1, 2),) : 8 ((0, 1, 2), (1, 0, 2)) : 99 soit 258 unions avec trous, 102 sans trou plx@sony:~/Bureau$
Le script crée également deux fichiers au format WKT qui permettent de localiser les bons et les mauvais cas (la répartition n'est d'ailleurs pas si aléatoire que ça : voir Fig3.png (bon cas, unions sans trou) et Fig4.png (mauvais cas, unions avec trous))

Si tu ne veux plus tomber sur ces problèmes qui arriveront fréquemment, tu as 2 façons de faire.

La plus simple et la moins couteuse est de filtrer, lors de chaque appel à la méthode union, les trous qui ont des aires très petites. Tu crées un nouveau polygone (via le constructeur) en ne mettant plus le(s) périmètre(s) concerné(s) et tu continues.

La seconde (celle qu'on utilise généralement quand on traite des données géographiques) consiste à prétraiter les données. Pour ce que tu veux faire, ça consisterait à faire toute une série d'intersections entre tes arcs (les périmètres intérieurs et extérieurs) et d'ajouter des points dans ces périmètres quand une intersection survient là où il n'existe pas de point dans les données d'origine. C'est un peu plus compliqué mais plus satisfaisant pour l'esprit.

06/12/2011, 09h46
awalter1

Merci beaucoup pour cette analyse complète et détaillée. J'ai compris vos explications et puisque ces effets de bord sont incontournables, je vais chercher une solution de contournement. Parmi les deux propositions :
1) travailler sur les aires très petites : le souci est que l'union de mes zones correspondent à des tailles équivalentes à l'Europe, donc j'imagine que deux segments très proches peuvent avec la leur longueur importante correspondent à des aires non négligeables. Je ne suis pas sûr de pouvoir me fier à l'aire pour en déduire que je suis dans un trou du style de celui que j'ai rencontré, mais faut voir. Je peux aussi me servir du fait que mes trous ne sont pas des polygones parfaits (pour moi) dans le sens ou ils ont des segments confondus, ce qui rends compte d'une configuration aux limites.
2) j'imagine que je vais trouver avec shapely la possibilité de calculer l'intersection entre deux segments, mais cela ne revient t'il pas à faire à la main ce que l'union ne parvient pas à faire correctement ? Mais je vais essayer.
3) n'y a t'il pas une autre solution qui serait de prendre la zone totale et de soustraire mes polygones, le résultat serait les trous (les vrais), à moins que j'introduise les mêmes erreurs.
Je vais explorer ces différentes possibilités et vous tiendrez au courant de mes résultats
Merci encore