shapely : union de polygones

Bonjour,
Dans mon application python 2.1.1, j'utilise la lib shapely 1.2.13.
Je fais l'union de polygones, et j'ai une erreur. Si je fais cette union dans un autre ordre je n'ai pas d'erreur. Voici l'extrait de mon code :

Code:

---

1 2 3 4 5 6 7 8 9

l_polygon = self.d_polygons.keys() print "d_polygons[ l_polygon[0]]=",self.d_polygons[ l_polygon[0] ] envelop = self.d_polygons[ l_polygon[0] ] i = 1 for p in reversed(l_polygon[1:]): points = self.d_polygons[ p ] print "i=%d l_polygon[i]=%s" % (i,l_polygon[i]) i += 1 envelop = envelop.union(points)

---

Le résultat est ici OK
si je fais :

Code:

---

1 2 3 4 5 6 7 8 9

l_polygon = self.d_polygons.keys() print "d_polygons[ l_polygon[0]]=",self.d_polygons[ l_polygon[0] ] envelop = self.d_polygons[ l_polygon[0] ] i = 1 for p in l_polygon[1:]: points = self.d_polygons[ p ] print "i=%d l_polygon[i]=%s" % (i,l_polygon[i]) i += 1 envelop = envelop.union(points)

---

j'ai l'erreur :

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

>linux awalter 48>: python ./HMI/holes_and_overlaps.py -u awalter -p ipas0 -s ACE2011B -d ALLSIMUL -a FCBALDB -m Y IPAS_ENV file used : </ipas2/opeint/IPAS/ACE2011B/IPAS/tmp/IPAS_ENV_NESU.data> d_polygons[ l_polygon[0]]= POLYGON ((186581.0000000000000000 23067.0000000000000000, 186621.0000000000000000 22995.0000000000000000, 186709.0000000000000000 23005.0000000000000000, 186725.0000000000000000 23003.0000000000000000, 186716.0000000000000000 23158.0000000000000000, 186727.0000000000000000 23412.0000000000000000, 186819.0000000000000000 23716.0000000000000000, 186840.0000000000000000 23805.0000000000000000, 186877.0000000000000000 23973.0000000000000000, 186920.0000000000000000 24402.0000000000000000, 186961.0000000000000000 24447.0000000000000000, 186879.0000000000000000 24756.0000000000000000, 186679.0000000000000000 25123.0000000000000000, 186177.0000000000000000 23879.0000000000000000, 186581.0000000000000000 23067.0000000000000000)) i=1 l_polygon[i]=A039 i=2 l_polygon[i]=A032 ... i=123 l_polygon[i]=A072 i=124 l_polygon[i]=A049B Traceback (most recent call last): File "./HMI/holes_and_overlaps.py", line 2270, in <module> hole_overlap = Application__(options.user, options.password, options.simulation, options.database, options.airspace, options.watch_dog, options.mmi) File "./HMI/holes_and_overlaps.py", line 134, in __init__ self.I_main = Main__(self) File "./HMI/holes_and_overlaps.py", line 943, in __init__ envelop = envelop.union(points) File "/opt/IPASTEST/lib/python2.7/site-packages/shapely/geometry/base.py", line 343, in union return geom_factory(self.impl['union'](self, other)) File "/opt/IPASTEST/lib/python2.7/site-packages/shapely/topology.py", line 53, in __call__ "This operation produced a null geometry. Reason: unknown") shapely.geos.TopologicalError: This operation produced a null geometry. Reason: unknown >

---

Je n'ai pas trop d'idée pour m'en sortir si ce n'est pas déterministe.
Merci

J'utilise la version 3.3.0 de geos.

Bonjour awalter1

Citation:

---

Envoyé par awalter1

si ce n'est pas déterministe.

---

Dans shapely/Geos, l'union de polygones est censée être indépendante de l'ordre dans lequel on traite les polygones.

Mais, la topologie et les opérations d'union, d'intersection, etc... est, pour ce que je connais, le domaine où on peut le plus trainer un bug pendant très longtemps et qu'il vous "pête" à la figure des mois après. La complexité des données traitées fait le reste.

Je suis preneur des polygones qui posent problème (format shapefile / wkt / wkb / autre ...) pour y jeter un oeil

Bonjour plxpy, content de vous retrouver sur le forum.
Dans la boucle qui fusionne mes polygones un à un, j'ai localisé le volume qui introduit le pb (à mon avis), j'ai imprimé l'enveloppe de départ et la valeur de ce polygone :

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12

l_polygon = self.d_polygons.keys() envelop = self.d_polygons[ l_polygon[0] ] for p in l_polygon[1:]: # add the current polygon points = self.d_polygons[ p ] print "envelop=",envelop print "polygon = ",p print points try: envelop = envelop.union(points) except TopologicalError, err: print "err=",err

---

Le résultat est le suivant :

Code:

---

1 2 3 4 5

envelop= MULTIPOLYGON (((186177.0000000000000000 23879.0000000000000000, 186679.0000000000000000 25123.0000000000000000, 187108.0000000000000000 25611.0000000000000000, 187443.0000000000000000 25800.0000000000000000, 187799.0000000000000000 25998.0000000000000000, 187850.0000000000000000 26030.0000000000000000, 188098.0000000000000000 26164.0000000000000000, 188060.0000000000000000 25420.0000000000000000, 187750.0000000000000000 24635.0000000000000000, 187707.0000000000000000 24769.0000000000000000, 187698.0000000000000000 24751.0000000000000000, 187627.0000000000000000 24529.0000000000000000, 187285.0000000000000000 24351.0000000000000000, 186961.0000000000000000 24447.0000000000000000, 186920.0000000000000000 24402.0000000000000000, 186877.0000000000000000 23973.0000000000000000, 186840.0000000000000000 23805.0000000000000000, 186819.0000000000000000 23716.0000000000000000, 186727.0000000000000000 23412.0000000000000000, 186716.0000000000000000 23158.0000000000000000, 186725.0000000000000000 23003.0000000000000000, 186709.0000000000000000 23005.0000000000000000, 186621.0000000000000000 22995.0000000000000000, 186581.0000000000000000 23067.0000000000000000, 186140.0000000000000000 21764.0000000000000000, 185860.0000000000000000 21973.0000000000000000, 185861.0000000000000000 22876.0000000000000000, 185906.0000000000000000 23297.0000000000000000, 186177.0000000000000000 23879.0000000000000000)), ((193200.0000000000000000 23400.0000000000000000, 198001.0000000000000000 23400.0000000000000000, 198000.0000000000000000 23115.0000000000000000, 198000.0000000000000000 22860.0000000000000000, 198000.0000000000000000 21746.0000000000000000, 198000.0000000000000000 21122.0000000000000000, 193078.0000000000000000 19813.0000000000000000, 192546.0000000000000000 19678.0000000000000000, 192600.0000000000000000 20040.0000000000000000, 192623.0000000000000000 20162.0000000000000000, 193200.0000000000000000 23400.0000000000000000)), ((190098.0000000000000000 17306.0000000000000000, 189907.0000000000000000 17336.0000000000000000, 189949.0000000000000000 17780.0000000000000000, 189997.0000000000000000 18164.0000000000000000, 190027.0000000000000000 18480.0000000000000000, 190082.0000000000000000 19023.0000000000000000, 190158.0000000000000000 19017.0000000000000000, 192535.0000000000000000 19616.0000000000000000, 192390.0000000000000000 18630.0000000000000000, 192380.0000000000000000 18580.0000000000000000, 192180.0000000000000000 17640.0000000000000000, 192120.0000000000000000 17400.0000000000000000, 192060.0000000000000000 17280.0000000000000000, 191940.0000000000000000 17160.0000000000000000, 191820.0000000000000000 17040.0000000000000000, 191790.0000000000000000 17035.0000000000000000, 191700.0000000000000000 17018.0000000000000000, 191399.0000000000000000 17098.0000000000000000, 190098.0000000000000000 17306.0000000000000000)), ((192900.0000000000000000 23701.0000000000000000, 192818.0000000000000000 23784.0000000000000000, 192682.0000000000000000 24020.0000000000000000, 192616.0000000000000000 24270.0000000000000000, 192277.0000000000000000 23790.0000000000000000, 191545.0000000000000000 22173.0000000000000000, 190084.0000000000000000 19031.0000000000000000, 190100.0000000000000000 19200.0000000000000000, 189120.0000000000000000 20040.0000000000000000, 189120.0000000000000000 22800.0000000000000000, 188734.0000000000000000 22800.0000000000000000, 188650.0000000000000000 22947.0000000000000000, 189006.0000000000000000 23789.0000000000000000, 189121.0000000000000000 24060.0000000000000000, 189143.0000000000000000 25345.0000000000000000, 189566.0000000000000000 25410.0000000000000000, 189785.0000000000000000 25430.0000000000000000, 190070.0000000000000000 25455.0000000000000000, 189950.0000000000000000 22800.0000000000000000, 190500.0000000000000000 23395.0000000000000000, 191919.0000000000000000 25494.0000000000000000, 191940.0000000000000000 25290.0000000000000000, 192011.0000000000000000 25177.0000000000000000, 192236.0000000000000000 25018.0000000000000000, 192508.0000000000000000 24709.0000000000000000, 192585.0000000000000000 24539.0000000000000000, 192599.0000000000000000 24415.0000000000000000, 192914.0000000000000000 24970.0000000000000000, 192900.0000000000000000 23701.0000000000000000))) polygon = A020A POLYGON ((190048.0000000000000000 22906.0000000000000000, 190070.0000000000000000 25455.0000000000000000, 189950.0000000000000000 22800.0000000000000000, 190048.0000000000000000 22906.0000000000000000))   envelop= MULTIPOLYGON (((186177.0000000000000000 23879.0000000000000000, 186679.0000000000000000 25123.0000000000000000, 187108.0000000000000000 25611.0000000000000000, 187443.0000000000000000 25800.0000000000000000, 187799.0000000000000000 25998.0000000000000000, 187850.0000000000000000 26030.0000000000000000, 188098.0000000000000000 26164.0000000000000000, 188060.0000000000000000 25420.0000000000000000, 187750.0000000000000000 24635.0000000000000000, 187707.0000000000000000 24769.0000000000000000, 187698.0000000000000000 24751.0000000000000000, 187627.0000000000000000 24529.0000000000000000, 187285.0000000000000000 24351.0000000000000000, 186961.0000000000000000 24447.0000000000000000, 186920.0000000000000000 24402.0000000000000000, 186877.0000000000000000 23973.0000000000000000, 186840.0000000000000000 23805.0000000000000000, 186819.0000000000000000 23716.0000000000000000, 186727.0000000000000000 23412.0000000000000000, 186716.0000000000000000 23158.0000000000000000, 186725.0000000000000000 23003.0000000000000000, 186709.0000000000000000 23005.0000000000000000, 186621.0000000000000000 22995.0000000000000000, 186581.0000000000000000 23067.0000000000000000, 186140.0000000000000000 21764.0000000000000000, 185860.0000000000000000 21973.0000000000000000, 185861.0000000000000000 22876.0000000000000000, 185906.0000000000000000 23297.0000000000000000, 186177.0000000000000000 23879.0000000000000000)), ((193200.0000000000000000 23400.0000000000000000, 198001.0000000000000000 23400.0000000000000000, 198000.0000000000000000 23115.0000000000000000, 198000.0000000000000000 22860.0000000000000000, 198000.0000000000000000 21746.0000000000000000, 198000.0000000000000000 21122.0000000000000000, 193078.0000000000000000 19813.0000000000000000, 192546.0000000000000000 19678.0000000000000000, 192600.0000000000000000 20040.0000000000000000, 192623.0000000000000000 20162.0000000000000000, 193200.0000000000000000 23400.0000000000000000)), ((190098.0000000000000000 17306.0000000000000000, 189907.0000000000000000 17336.0000000000000000, 189949.0000000000000000 17780.0000000000000000, 189997.0000000000000000 18164.0000000000000000, 190027.0000000000000000 18480.0000000000000000, 190082.0000000000000000 19023.0000000000000000, 190158.0000000000000000 19017.0000000000000000, 192535.0000000000000000 19616.0000000000000000, 192390.0000000000000000 18630.0000000000000000, 192380.0000000000000000 18580.0000000000000000, 192180.0000000000000000 17640.0000000000000000, 192120.0000000000000000 17400.0000000000000000, 192060.0000000000000000 17280.0000000000000000, 191940.0000000000000000 17160.0000000000000000, 191820.0000000000000000 17040.0000000000000000, 191790.0000000000000000 17035.0000000000000000, 191700.0000000000000000 17018.0000000000000000, 191399.0000000000000000 17098.0000000000000000, 190098.0000000000000000 17306.0000000000000000)), ((192900.0000000000000000 23701.0000000000000000, 192818.0000000000000000 23784.0000000000000000, 192682.0000000000000000 24020.0000000000000000, 192616.0000000000000000 24270.0000000000000000, 192277.0000000000000000 23790.0000000000000000, 191545.0000000000000000 22173.0000000000000000, 190084.0000000000000000 19031.0000000000000000, 190100.0000000000000000 19200.0000000000000000, 189120.0000000000000000 20040.0000000000000000, 189120.0000000000000000 22800.0000000000000000, 188734.0000000000000000 22800.0000000000000000, 188650.0000000000000000 22947.0000000000000000, 189006.0000000000000000 23789.0000000000000000, 189121.0000000000000000 24060.0000000000000000, 189143.0000000000000000 25345.0000000000000000, 189566.0000000000000000 25410.0000000000000000, 189785.0000000000000000 25430.0000000000000000, 190070.0000000000000000 25455.0000000000000000, 190048.0001584033016115 22906.0183531817456242, 190500.0000000000000000 23395.0000000000000000, 191919.0000000000000000 25494.0000000000000000, 191940.0000000000000000 25290.0000000000000000, 192011.0000000000000000 25177.0000000000000000, 192236.0000000000000000 25018.0000000000000000, 192508.0000000000000000 24709.0000000000000000, 192585.0000000000000000 24539.0000000000000000, 192599.0000000000000000 24415.0000000000000000, 192914.0000000000000000 24970.0000000000000000, 192900.0000000000000000 23701.0000000000000000)))

---

Le polygone A02A est simple et "propre" (pas de chiffres après la virgule) et pourtant le résultat de cette union introduit des valeurs après la virgule. J'ai du mal à comprendre.
Je ne sais pas si cela vous permets de tester. Sous quelle forme puis je vous soumettre mes données vous (moi je les lis à partir d'une database) ?
Merci

Je ne pense pas qu'il y ait d'erreur sur les exemples que tu me fournis.

4 images à la suite, dans l'ordre :

1. l'enveloppe avant (vert) et le polygone (rouge)
2. petit zoom sur la partie "intéressante"
3. très fort zoom sur cette même partie
4. même endroit dans l'enveloppe résultat

http://plx.free.fr/general.png
http://plx.free.fr/envavant_poly.png
http://plx.free.fr/zoom.png
http://plx.free.fr/envapres.png

Les coordonnées non entières que tu récupères sont celles du point d'intersection d'un segment de l'enveloppe (vert) et d'un segment du polygone (rouge) :

Code:

---

1 2 3 4 5 6 7

from shapely.geometry import LineString >>> segment_enveloppe = LineString(((189950.,22800.),(190500.,23395.))) >>> >>> segment_polygone = LineString(((190070.,25455.),(190048.,22906))) >>> >>> segment_enveloppe.intersection(segment_polygone).wkt 'POINT (190048.0001584033016115 22906.0183531817456242)'

---

Tu retrouves les coordonnées non entières de l'enveloppe résultat.

Si, par contre, lors des unions successives, et selon l'ordre, une exception est levée, c'est un autre problème. Pour y jeter un coup d'oeil, il me faudrait tous les polygones (le format utilisé me convient)

chapeau!!!

bravo pour cette analyse.
Voici les volumes dont l'union doit être faite, je les mets en fichiers car sinon le post est refusé car trop long et en 2 fichiers car je dépasse le max autorisé.
Merci de votre implication

Pouvez vous me dire ce que vous utilisez pour visualiser les volumes ? Ca m'interesse (surtout le zoom) Est ce un free ware ? Un développement perso ? Est ce une interface en python ?
Merci

J'utilise OpenJump, écrit en Java, qui s'appuye sur JTS, Java Topology Suite (la base pour GEOS et Shapely) et qui tourne aussi bien sous Windows, Linux et MacOSX.

Vu les messages d'erreurs de ton premier post, l'ordre des polygones dans le calcul de l'union est lié à une boucle sur les clés d'un dictionnaire (les 'A017', 'A011', etc...) donc pas d'ordre évident.

Tu pourrais m'indiquer l'ordre de ces clés (dans ton script qui plante, à chaque itération, faire afficher la clé avant de calculer l'union, sans oublier la première des clés pour le polygone qui te sert à initialiser l'union) ? Merci

Je viens de l'installer. Pouvez vous me donner le fichier que vous avez utilisé avec mes données afin que je vois le format utilisé.
merci

Citation:

---

Envoyé par plxpy

Vu les messages d'erreurs de ton premier post, l'ordre des polygones dans le calcul de l'union est lié à une boucle sur les clés d'un dictionnaire (les 'A017', 'A011', etc...) donc pas d'ordre évident.

Tu pourrais m'indiquer l'ordre de ces clés (dans ton script qui plante, à chaque itération, faire afficher la clé avant de calculer l'union, sans oublier la première des clés pour le polygone qui te sert à initialiser l'union) ? Merci

---

En fait l'union de mes polygones me sert pour calculer l'enveloppe de l'ensemble :

Code:

---

envelop = envelop.envelope

---

ensuite je fais une boucle pour enlever mes volumes un à un : le résultat doit être vide sinon il y a un trou : c'est mon objectif premier de trouver les trous.
Donc l'ordre de l'union n'était pas un souci pour moi (c'est peut être une erreur de ma part) : Afin que la boucle d'union ne soit pas aléatoire, je vais trier la liste des polygones par le nom des polygones. Je fais ça tout de suite ...
Cordialement

Lorsque mes volumes sont triés par nom, je n'ai plus le cas "null geometry", dans l'ordre inverse c'est pareil. Je vais retrouver un ordre de traitement qui redéclenche l'erreur ...

ça y est avec le code suivant, j'initialise en dur la liste des polygones :

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12

# get the sectors envelope # l_polygon = self.d_polygons.keys() l_polygon = ['A038', 'A039', 'A032', 'A033', 'A030', 'A031', 'A036', 'A037', 'A034', 'A035', 'A020B', 'A020C', 'A020A', 'A012A', 'A025', 'A024D', 'A027', 'A091E', 'A021B', 'A021', 'A020', 'A008D', 'A022', 'A008A', 'A008B', 'A008C', 'A021D', 'A021A', 'A021C', 'A060', 'A051A', 'A077', 'A016J', 'A016I', 'A016H', 'A016G', 'A016F', 'A016E', 'A016D', 'A016C', 'A016B', 'A016A', 'A075A', 'A075B', 'A019A', 'A011B', 'A022A', 'A022B', 'A011A', 'A035A', 'A024F', 'A008', 'A091F', 'A024', 'A091D', 'A026', 'A091B', 'A091C', 'A023', 'A091A', 'A024E', 'A029', 'A028', 'A018A', 'A023C', 'A023B', 'A023A', 'A023G', 'A023F', 'A023E', 'A023D', 'A010A', 'A006B', 'A053A', 'A053B', 'A006A', 'A049C', 'A050', 'A051', 'A052', 'A053', 'A090A', 'A055', 'A056', 'A090B', 'A058', 'A059', 'A057', 'A024B', 'A024C', 'A024A', 'A003A', 'A024G', 'A003C', 'A003B', 'A040A', 'A013A', 'A043', 'A042', 'A041', 'A040', 'A047', 'A046', 'A045', 'A044', 'A049', 'A048', 'A025A', 'A025C', 'A025B', 'A025E', 'A025D', 'A025G', 'A025F', 'A004A', 'A004B', 'A092A', 'A078', 'A079', 'A076', 'A049F', 'A049E', 'A075', 'A072', 'A049B', 'A049A', 'A071', 'A009C', 'A009B', 'A009A', 'A018', 'A019', 'A009D', 'A014', 'A015', 'A016', 'A017', 'A010', 'A011', 'A012', 'A013', 'A090', 'A091', 'A092', 'A093', 'A050A', 'A050C', 'A050B', 'A050D', 'A015B', 'A026A', 'A026B', 'A026C', 'A026D', 'A026E', 'A026F', 'A026G', 'A026H', 'A058A', 'A061', 'A054', 'A063', 'A062', 'A065', 'A064', 'A067', 'A009', 'A015C', 'A007', 'A006', 'A005', 'A004', 'A003', 'A002', 'A001', 'A000', 'A015A', 'A083', 'A082', 'A081', 'A080', 'A089', 'A014A', 'A014B', 'A015D', 'A049D', 'A067B', 'A067A']   envelop = self.d_polygons[ l_polygon[0] ] for p in l_polygon[1:]: # add the current polygon points = self.d_polygons[ p ] try: envelop = envelop.union(points) except TopologicalError, err: print "Warning : with volume %s : %s" %(p,err)

---

C'est le volume A049B qui provoque le "null geometry"
Cordialement

OpenJump reconnait le format WKT (format texte) parmi d'autres. La représentation wkt d'une géométrie shapely s'obtient via l'attribut wkt

Code:

---

1 2 3 4 5

>>> wkt = open('polygones.wkt','w') >>> for poly in polygons: ... print >> wkt, poly.wkt ... >>> wkt.close()

---

Ce n'est pas le plus compact des formats mais, pour tes polygones, ce sera amplement suffisant et il a l'avantage d'être simple. C'est ce format que j'ai utilisé pour les affichages précédents.

Pour le problème de fond, je continue(rai) à chercher. Le fait que l'ordre des polygones fasse parfois planter le calcul n'est pas satisfaisant. Afin d'orienter les recherches, est-il possible de connaître ce que représente les polygones ? des objets gépgraphiques ? et quelle est l'unité des coordonnées ?

est-il normal d'avoir des polygones qui "se touchent presque" mais pas parfaitement ?

ps : pour les fichiers au format wkt, utilise bien l'extension .wkt. Ca te facilitera la vie

Je vais essayer OPenJump demain.
Mes polygones représentent des volumes de l'espace aérien. On raisonne à altitude constante. Le x,y corresponds respectivement à la latitude/longitude converti en secondes d'arc :
latitude 1650N06 => x = 16*3600 + 50*60 + 6 et N détermine un signe - ou +
idem pour longitude.
Logiquement les polygones doivent être joints: pas de trous ni d'overlaps dans l'ensemble (c'est le but de mon travail).
Cordialement

J'ai oublié aussi une contrainte forte sur les volumes: si on a 3 points consécutifs communs entre V1 et V2 : A - B - C, alors V1 et V2 doivent utiliser les segments A-B et B-C (on ne peut pas avoir A-C pour un volume et A-B et B-C pour l'autre même si les 3 points sont alignés). C'est pour des problèmes de projection où la projection de A - C n'est pas égale à la projection de A - B et B - C : ça peut introduire un trou entre V1 et V2 (la terre est courbe et les distances sont assez grandes).

Dans mon algorithme, je détermine les altitudes qui interviennent dans la définition des volumes, et pour chaque tranche d'altitude je traite les volumes concernés par cette altitude.

Pour répondre à votre question : on ne doit pas avoir des polygones qui se touchent presque (le presque est peut être couvert par un polygone ?), mais on peut avoir des polygones avec des formes tordus et des portions fines (1 à 2 secondes d'arc, ce qui fait 1,8 à 3,2 km) sur une surface totale de 500 à 2000 km de coté.

Cordialement

Cordialement

Un petit point sur mes recherches

Pour commencer, j'ai trouvé un certain nombre de polygones qui étaient strictement égaux. Je ne sais pas dire si c'est normal ou pas. Il y a 188 polygones/identifiants mais seulement 108 polygones distincts (question géométrie)

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

from shapely.wkt import loads   def load_polygons(): polygons = dict() for name in('volumes1.txt','volumes2.txt'): with open(name) as f: for line in f: ident = line.split("'")[1] geom = loads(line[line.find('POLYGON'):]) polygons[ident] = geom return polygons   def distinct_polygons(polygons):   """les cles sont des tuples formes des identifiants tries du/des polygone(s)"""   distinct = dict()   for curid, curgeom in polygons.iteritems():   found = False for idents, geometry in distinct.iteritems(): if curgeom.equals(geometry): old_key = idents new_key = tuple(sorted(idents+(curid,))) found = True break   if found: distinct[new_key] = distinct[old_key] del distinct[old_key]   else: distinct[(curid,)] = curgeom   return distinct

---

A l'utilisation :

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

>>> from awalter import * >>> polygons = load_polygons() >>> len(polygons) 188 >>> >>> dist_polygons = distinct_polygons(polygons) >>> len(dist_polygons) 108 >>> # les identifiants des polygones geometriquement identiques ... for k in dist_polygons: ... if len(k) != 1: ... print k ... ('A033', 'A048') ('A016C', 'A016F') ('A028', 'A043') ('A082', 'A083') ('A035A', 'A053A', 'A053B') ('A000', 'A036', 'A054') ('A004B', 'A011B', 'A015D') ('A050', 'A050C') ('A009A', 'A009B', 'A009D') ('A006A', 'A006B') ('A008A', 'A008C') ('A032', 'A047') ('A007', 'A071') ('A030', 'A045') ('A023E', 'A024E') ('A023D', 'A024D') ('A016A', 'A016J') ('A039', 'A057') ('A029', 'A044') ('A014', 'A067') ('A004', 'A011') ('A016E', 'A016I') ('A062', 'A063') ('A020A', 'A021A') ('A014B', 'A067B') ('A038', 'A056') ('A004A', 'A011A') ('A080', 'A081', 'A093') ('A025G', 'A026G') ('A049', 'A049E') ('A023F', 'A024F') ('A014A', 'A067A') ('A049A', 'A049D') ('A031', 'A046') ('A023C', 'A024C') ('A010', 'A010A') ('A025', 'A026') ('A005', 'A012', 'A012A') ('A025F', 'A026F') ('A025D', 'A026D', 'A026H') ('A076', 'A077') ('A035', 'A053') ('A025B', 'A026B') ('A025E', 'A026E') ('A018', 'A019') ('A025C', 'A026C') ('A009', 'A009C') ('A027', 'A042') ('A025A', 'A026A') ('A016B', 'A016G', 'A023G', 'A024G') ('A041', 'A059') ('A023B', 'A024B') ('A002', 'A037', 'A055') ('A013', 'A013A') ('A078', 'A079') ('A040', 'A058') ('A017', 'A034') ('A020', 'A021') ('A020C', 'A021C', 'A021D') ('A050B', 'A050D') ('A049C', 'A049F') ('A018A', 'A019A') ('A023', 'A024') ('A020B', 'A021B') ('A016D', 'A016H') ('A040A', 'A058A') ('A023A', 'A024A') ('A001', 'A061') ('A008B', 'A008D') >>>

---

Sinon, je reproduis bien la levée d'exception avec la liste d'identifiants que tu as fournie au post#14 (la liste l_polygon, ligne 3 du bout de code) à l'itération 124 lors de l'union de l'enveloppe courante avec le polygone 'A049B' (ma liste est dans la variable ordre) :

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

>>> import sys >>> enveloppe = Polygon() >>> for i, ident in enumerate(ordre): ... try: ... enveloppe = enveloppe.union(polygons[ident]) ... print >> sys.stderr, "%s U" % ident, ... except: ... # message d'erreur ... print >> sys.stderr, "\nechec indice %d ident %s" % (i,ident) ... # l'enveloppe en l'etat, avant tentative d'union ... print >> open('enveloppe_avant_echec.wkt','w'), enveloppe.wkt ... # le polygone en cause ... print >> open('%s.wkt' % ident,'w'), polygons[ident].wkt ... break ... A038 U A039 U A032 U A033 U A030 U A031 U A036 U A037 U A034 U A035 U A020B U A020C U A020A U A012A U A025 U A024D U A027 U A091E U A021B U A021 U A020 U A008D U A022 U A008A U A008B U A008C U A021D U A021A U A021C U A060 U A051A U A077 U A016J U A016I U A016H U A016G U A016F U A016E U A016D U A016C U A016B U A016A U A075A U A075B U A019A U A011B U A022A U A022B U A011A U A035A U A024F U A008 U A091F U A024 U A091D U A026 U A091B U A091C U A023 U A091A U A024E U A029 U A028 U A018A U A023C U A023B U A023A U A023G U A023F U A023E U A023D U A010A U A006B U A053A U A053B U A006A U A049C U A050 U A051 U A052 U A053 U A090A U A055 U A056 U A090B U A058 U A059 U A057 U A024B U A024C U A024A U A003A U A024G U A003C U A003B U A040A U A013A U A043 U A042 U A041 U A040 U A047 U A046 U A045 U A044 U A049 U A048 U A025A U A025C U A025B U A025E U A025D U A025G U A025F U A004A U A004B U A092A U A078 U A079 U A076 U A049F U A049E U A075 U A072 U echec indice 124 ident A049B >>>

---

A quoi ressemblent l'enveloppe à ce moment là (en vert) et le polygone 'A049B' (rouge) :

http://plx.free.fr/dvp_images/enveloppe_avant_echec.png

En zoomant :

http://plx.free.fr/dvp_images/zoom_A049B.png


Géométriquement, l'enveloppe est un multi-polygone formés de 2 polygones :

• le premier avec 6 trous (dont un "micro" trou qu'on aperçoit sur le zoom plus haut)
• le second sans trou.

C'est le premier polygone troué qui provoque l'exception, pas le second :

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

>>> type(enveloppe) <class 'shapely.geometry.multipolygon.MultiPolygon'> >>> len(enveloppe.geoms) 2 >>> for g,geom in enumerate(enveloppe.geoms): ... print "polygone indice %d, %d trou(s)" % (g,len(geom.interiors)) ... polygone indice 0, 6 trou(s) polygone indice 1, 0 trou(s) >>> >>> enveloppe.geoms[0].union(polygons['A049B']) Traceback (most recent call last): File "<stdin>", line 1, in <module> File "/usr/local/lib/python2.7/dist-packages/Shapely-1.2.14-py2.7-linux-i686.egg/shapely/geometry/base.py", line 343, in union return geom_factory(self.impl['union'](self, other)) File "/usr/local/lib/python2.7/dist-packages/Shapely-1.2.14-py2.7-linux-i686.egg/shapely/topology.py", line 53, in __call__ "This operation produced a null geometry. Reason: unknown") shapely.geos.TopologicalError: This operation produced a null geometry. Reason: unknown >>> >>> enveloppe.geoms[1].union(polygons['A049B']) <shapely.geometry.multipolygon.MultiPolygon object at 0x880e3cc> >>>

---

En tenant compte du périmètre extérieur du premier polygone de l'enveloppe et d'un trou à la fois, c'est celui d'indice 5 (celui du zoom, ça se trouve sous OpenJump) qui provoque plus particulièrement l'exception :

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

>>> for i in range(6): ... poly_un_trou = Polygon(enveloppe.geoms[0].exterior,(enveloppe.geoms[0].interiors[i],)) ... print "avec trou indice", i, ... try: ... _ = poly_un_trou.union(polygons['A049B']) ... print "ok" ... except: ... print "CA PLANTE" ... avec trou indice 0 ok avec trou indice 1 ok avec trou indice 2 ok avec trou indice 3 ok avec trou indice 4 ok avec trou indice 5 CA PLANTE

---

J'ai même fait toutes les combinaisons possibles, avec 1 trou, 2 trous, ... , les 6 trous : toutes les combinaisons n'incluant pas ce trou d'indice 5 ne lèvent pas d'exception et toutes les autres si :

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

>>> ooops = list() >>> from itertools import permutations >>> exterior = enveloppe[0].exterior >>> interiors = enveloppe[0].interiors >>> for nbtrous in range(1,7): ... for perm in permutations(range(6),nbtrous): ... holes = [ interiors[i] for i in perm ] ... newPoly = Polygon(exterior,holes) ... try: ... _ = newPoly.union(polygons['A049B']) ... if 5 in perm: ... print perm, "fonctionne" ... except: ... ooops.append(perm) ... >>> >>> len(ooops) 1631 >>> >>> for perm in ooops: ... assert 5 in perm ... >>>

---

Qu'a-t-il de particulier, en lui-même et/ou par rapport au polygone A049B ?


Impossible de zoomer sur ces deux éléments et d'avoir une vue d'ensemble car ils ont des dimensions vraiment trop différentes. Un schéma, pas à l'échelle, mais topologiquement juste et qui respecte les positions relatives des points constituant ces 2 éléments :

http://plx.free.fr/dvp_images/schema.png


Le trou a deux points ayant le même X (ici même latitude), séparés de -3.64e-12 (secondes d'arc) en longitude :

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

>>> coords = tuple(enveloppe.geoms[0].interiors[5].coords) >>> for i in range(len(coords)): ... print "x=%.2f y=%.2f" % (coords[i][0],coords[i][1]), ... if i == len(coords)-1: ... print ... else: ... delta_x = coords[i+1][0]-coords[i][0] ... delta_y = coords[i+1][1]-coords[i][1] ... print "delta_x=%10.2e delta_y=%10.2e" % (delta_x,delta_y) ... x=188859.78 y=23443.17 delta_x= 1.46e+02 delta_y= 3.46e+02 x=189006.00 y=23789.00 delta_x= -1.46e+02 delta_y= -3.46e+02 x=188859.78 y=23443.17 delta_x= 0.00e+00 delta_y= -3.64e-12 x=188859.78 y=23443.17 >>>

---

Ce qui fait (on est à une latitude d'environ 189000./3600. = 52.5°), à la louche :

1° <==> 111 km à l'équateur
1° <==> 111 x cos(52.5°) sur la zone d'étude

à 3.64e-12 secondes d'arc, on titille 1e-10 mètres ...


Autre question importante (enfin, qui m'importe quand je trouve des "bizarreries" en travaillant avec Shapely) : est-ce du à Shapely lui-même ? ou est-ce que le problème est déjà présent dans la bibliothèque GEOS ?

Réponse avec ce code (vite fait) en C :

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

#include <stdio.h> #include <stdlib.h> #include <assert.h> #include <geos_c.h>   #define NPOLYS 125 #define BUFSIZE 2500   int main(int argc,char *argv[]) { GEOSGeometry *(polygons[NPOLYS]); GEOSGeometry *union1; GEOSGeometry *union2; GEOSWKTReader *reader; GEOSWKTWriter *writer; FILE *fp; char buffer[BUFSIZE]; int i;   initGEOS((void *)printf, (void *)printf); reader = GEOSWKTReader_create(); writer = GEOSWKTWriter_create();   // load the NPOLYS polygons from polygons.wkt assert(fp = fopen("polygons.wkt","r")); for (i=0; i<NPOLYS; i++) { fgets(buffer,BUFSIZE,fp); polygons[i] = GEOSWKTReader_read(reader,buffer); } fclose(fp);   union1 = GEOSGeom_clone(polygons[0]); for (i=1; i<NPOLYS; i++) { FILE *fp; char name[32];   fprintf(stderr,"union with polygon %d ...", i); assert(union2 = GEOSUnion(union1,polygons[i])); fprintf(stderr,"done\n");   sprintf(name,"union_%03d.wkt",i); fp = fopen(name,"w"); fprintf(fp,"%s",GEOSWKTWriter_write(writer,union2)); fclose(fp);   GEOSGeom_destroy(union1); union1 = union2; } }

---

et l'exécution :

Code:

---

1 2 3 4 5 6 7 8 9

plx@satellite:~/Bureau/awalter/valid_C$ ./awalter union with polygon 1 ...done union with polygon 2 ...done union with polygon 3 ...done ... je saute des lignes union with polygon 122 ...done union with polygon 123 ...done union with polygon 124 ...awalter: awalter.c:64: main: Assertion `union2 = GEOSUnion(union1,polygons[i])' failed. TopologyException: no outgoing dirEdge found at 189006 23789Abandon (core dumped)

---

ps: les coordonnées fournies dans le message d'erreur sont celles du point commun entre le trou et le polygone (cf schéma), j'ai vérifié

Conclusion : le problème survient "à la base", dans GEOS.

Je pense que la topologie des polygones n'a pas d'impact ici et que la cause du plantage est uniquement dû à des problèmes de calcul numérique.

Je continue à cogiter.

Bonjour,
Une analyse toujours aussi précise et complete :ccool:.
Qu'il y ait des polygones identiques, cela est assez courant. Il faut tenir compte de l'altitude associé à chacun d'eux : on part d'un espace en 3D, qu'on réduit à du 2D par tranche d'altitude. Dans ma problèmatique, je pourrais raisonner par tranche d'altitude et donc éviter cette duplication de polygones, mais comme je voulais l'enveloppe de l'ensemble, ce n'était pas strictement nécessaire. Je peux aussi tester avant chaque union si un polygone identique n'a pas été déjà traité.
Sur le problème identifié, étant donné que notre unité est la seconde d'arc, soit 1852m / 60 = 30,86 m, ne peut t'on pas considérer que toutes distances inférieurs sont à rejeter purement et simplement, compte tenu que l'on va au bout de la précision des outils utilisés, ici, geos ?
Merci

préparé avant d'avoir pris connaissance du post précédent - je le mets quand même. Je re-posterai plus tard

C'est bien un problème de précision de calcul numérique.

J'ai fait une fonction qui relance le calul de l'union entre l'enveloppe courante (juste avant le crash) et le polygone A049B en deplaçant de delta_Y le 2ieme point du micro-trou (voir schema post précédent) :

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

def new_try(enveloppe,A049B,delta_Y):   # coordonnees du perimetre exterieur du premier polygone de l'enveloppe exterior = tuple(enveloppe[0].exterior.coords)   # coordonnees des 5 premiers trous (du premier polygone de l'enveloppe) holes = [ tuple(enveloppe[0].interiors[i].coords) for i in range(5) ]   # les coordonnees du "micro-trou" qui plante l'union micro_hole = list(enveloppe[0].interiors[5].coords)   # ajout du delta_Y au point d'indice 2 (voir schema post precedent) micro_hole[2] = (micro_hole[2][0],micro_hole[2][1]+delta_Y)   # ajout de ce nouveau micro-trou a la liste holes.append(micro_hole) new_poly = Polygon(exterior,holes)   try: new_union = new_poly.union(A049B) print "CA MARCHE !!!" return new_poly, new_union   except: print "ca coince toujours ..."

---

J'ai ajouté n fois 1e-12, de 0 a 15 fois : à partir de 10 (soit 1e-11) ça passe (1e-11 secondes d'arc, c'est du, pour fixer les idées, 1e-10 METRES)

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

>>> for i in range(15): ... delta = float(i) * pow(0.1,12) ... print "essai avec + %d x 1e-12 :" % i, ... _ = new_try(enveloppe,A049B,delta) ... essai avec + 0 x 1e-12 : ca coince toujours ... essai avec + 1 x 1e-12 : ca coince toujours ... essai avec + 2 x 1e-12 : ca coince toujours ... essai avec + 3 x 1e-12 : ca coince toujours ... essai avec + 4 x 1e-12 : ca coince toujours ... essai avec + 5 x 1e-12 : ca coince toujours ... essai avec + 6 x 1e-12 : ca coince toujours ... essai avec + 7 x 1e-12 : ca coince toujours ... essai avec + 8 x 1e-12 : ca coince toujours ... essai avec + 9 x 1e-12 : ca coince toujours ... essai avec + 10 x 1e-12 : CA MARCHE !!! essai avec + 11 x 1e-12 : CA MARCHE !!! essai avec + 12 x 1e-12 : CA MARCHE !!! essai avec + 13 x 1e-12 : CA MARCHE !!! essai avec + 14 x 1e-12 : CA MARCHE !!!

---

Verifications visuelles

Code:

---

1 2 3 4

>>> new_poly, new_union = new_try(enveloppe,A049B,1e-11) CA MARCHE !!! >>> print >> open('new_poly.wkt','w'), new_poly.wkt >>> print >> open('new_union.wkt','w'), new_union.wkt

---

Avec le même très gros zoom

le micro-trou avant :

http://plx.free.fr/dvp_images/micro_trou_avant.png

et après :

http://plx.free.fr/dvp_images/micro_trou_apres.png

L'union résultat, zoomée sur la zone d'étude :

http://plx.free.fr/dvp_images/nouvel...sans_A049B.png

avec le polygone A049B, on voit bien que le trou a été, en partie, rebouché

http://plx.free.fr/dvp_images/nouvelle_union.png


Et maintenant le pourquoi. En fait, c'est sur mon schéma d'avant que j'ai, ici, surchargé :

http://plx.free.fr/dvp_images/schema_annote.png

Si les deux points de même X (même latitude) du micro-trou sont très proches (en 1e-12), les points d'intersection entre :

• les segments rouge et bleu
• les segments vert et bleu

le sont encore plus. Je ne suis pas rentré dans le détail de la bibliothèque GEOS mais ces points sont obligatoirement calculés à un moment ou à un autre lors de l'union.

Et si ces points, vu les valeurs numériques, avaient, par calcul, les mêmes coordonnées, ca mettrait forcément une belle pagaille. Eh bien, c'est ce qui se passe :

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

>>> micro_trou_coords = tuple(enveloppe.geoms[0].interiors[5].coords) >>> >>> rouge = LineString(micro_trou_coords[:2]) >>> vert = LineString(micro_trou_coords[1:3]) >>> >>> poly_049B_coords = tuple(polygons['A049B'].exterior.coords) >>> >>> bleu = LineString(poly_049B_coords[1:3]) >>> >>> inter_rouge = rouge.intersection(bleu) >>> inter_vert = vert.intersection(bleu) >>> >>> inter_rouge, inter_vert (<shapely.geometry.point.Point object at 0x6743d0>, <shapely.geometry.point.Point object at 0x674210>) >>> >>> inter_rouge.equals(inter_vert) True >>> >>> inter_rouge.coords[0] == inter_vert.coords[0] True >>>

---

Conclusion : certaines des erreurs que tu souhaites mettre en évidence (vu le domaine, ce micro-trou est assurément un artefact) sont tellement petites qu'elles mettent en échec la méthode utilisée (unions successives des polygones).


Sinon, il y a bien, dans shapely, la fonction "cascaded_union" (principe décrit ici) qui calcule l'union d'une liste de polygones. Ici, elle ne plante pas et produit un polygone avec, quand même, un autre micro-trou :

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13

>>> from shapely.ops import cascaded_union >>> enveloppe = cascaded_union(polygons.values()) >>> len(enveloppe.interiors) 1 >>> for xy in enveloppe.interiors[0].coords: ... print xy ... (184951.04588681934, 13094.414817494975) (184951.04588681937, 13094.414817494995) (184951.04588681934, 13094.414817495006) (184951.04588681934, 13094.414817494991) (184951.04588681934, 13094.414817494975) >>>

---

MAIS, avec cascaded_union :

1. tu ne maitrises pas l'ordre des unions : tu n'es pas à l'abri d'avoir, un jour, selon les polygones en entrée, un plantage comme celui que tu as rencontré
2. ici, le micro-trou qui nous pose problème n'existe plus dans le résultat (un ou plusieurs polygones le bouchent) et je ne suis pas sur que c'est ce que tu souhaites (pas de trou, pas de recouvrement) car il met en évidence un problème dans tes polygones en entrées (problème masqué avec cascaded_union !)

Bref, cascaded_union me semble être une fausse bonne idée ici.

Une première question pour être bien sur de comprendre ce à quoi tu veux arriver.

Tu travailles à la base sur des volumes (de l'espace aérien) puis tu en fais différentes coupes à altitude constante ce qui fournit les différents polygones, chaque polygone étant donc attaché à une altitude (tu l'indiques quelque part).

Quel est le but de la démarche qui consiste à faire, ensuite, l'union des polygones, toutes altitudes de coupe confondues ? Est-ce pour déterminer l'emprise au sol de l'espace aérien, de la zone de terrain survolée/survolable ?

Si c'est cela, tu as tout intêret à :

• d'abord travailler à altitude constante en faisant l'union des polygones attachés à cette altitude, uniquement ceux-là,
• faire, à la fin, l'union de ces unions intermédiaires

Je pense que ça peut t'éviter des cas de figures "dangereux", numériquement parlant (volumes avec des faces quasiment verticales, tranchées selon des z très proches)

En plus, c'est à altitude constante qu'il faut rechercher les éventuels trous mais, surtout, les éventuels recouvrements : selon les volumes, ca peut se recouvrir dans tous les sens sans qu'il y ait d'intersection entre tes volumes, tout simplement parce que tu traites tous les polygones en même temps.

Par exemple, un seul volume formé par un 'L' en 3D : si tu le découpes tout en haut et tout en bas, tu as 2 polygones qui s'intersectent (l'un est même contenu dans l'autre) sans qu'il y ait de problème.

Enfin, pour être "couvrant" (dans ta recherche de trous et de recouvrements, j'ai l'impression qu'il faut "trancher" tes volumes à toutes les altitudes (différentes) des points définissant les volumes. S'il n'y en a pas trop (j'ai l'impression que c'est le cas) j'ai parlé trop vite


Ensuite, concernant les précisions de calcul et ton souhait de travailler à une tolérance près.

Le seul endroit où Shapely prend en compte une tolérance est la méthode almost_equals qui indique si deux géométries sont égales à N décimales près. On ne peut pas en faire grand-chose ici, il me semble.

Mais même si tu travailles en secondes d'arc entières (je comprends que tu définis les différents points/sommets de tes volumes avec des (x,y) entiers, des réels, certes, mais valeurs entieres), je ne vois pas trop où on pourrait prendre en compte une quelconque tolérance.

Même des petits trous (ou des petits recouvrements, voir un de mes premiers messages) peuvent être révélateur d'un problème bien réel. Donc éliminer des trous, des recouvrements ou retoucher des coordonnées sur de tels critères (précision, arrondis à la seconde d'arc entière la plus proche, ...) ne me semble pas une bonne idée et plutôt contraire à la vérification des contraintes dont tu parlais au post#17)


Tout dernier point : tu parles enveloppe, trous, recouvrements mais je ne vois pas trop la logique des choses ni l'ordre dans lequel tu fais tout ça. Tu pourrais préciser ?

Aurais-tu également l'altitude des différents polygones dont tu as fournis la représentation wkt plus haut ?

Bonjour,
Je vais tenter d'expliquer sur quoi je travaille : tout d'abord cela concerne de la simulation, il y a donc des simplifications opérationnelles au cas ou tu serais de la partie.
Un espace aérien de simulation est défini par une zone rectangulaire d'altitude 0 à 999 (centaines de pieds) (en fait on parle de niveau de vol, mais disons altitude pour simplifier).
Notre espacé aérien de simulation est défini par des volumes. Chaque volume est un polygone + une altitude min et une max. L'ensemble des altitudes min et max définit les différentes altitudes sur lesquels on va raisonner. On a donc des volumes qui vont concerner une ou plusieurs tranches d'altitude.
Pour le simulateur : cet espace aérien ne doit comporter ni trous, ni overlaps, ni segments communs (on va écarter ce dernier cas qui n'est pas traité avec shapely).

Mon premier traitement a consisté à déterminer l'espace aérien de simulation :
1er temps : j'ai donc fait l'union de tous les volumes avec le pb que tu connais. J'ai changé mon algo, et aujourd'hui je prends un objet multipoint avec tous les points qui définissent mes volumes, puis je fait

Code:

---

1 2	multipoint = MultiPoint(self.l_envelop) envelop = multipoint.envelope

---

j'ai donc ma zone rectangulaire, l'altitude n'a pas été prise en compte ici (car il n'y en a pas au niveau des points)
2ème temps : recherche des trous

Code:

---

1 2 3 4 5 6 7

pour chaque altitude pour chaque volume si volume concerné par l altitude on fait la somme des volumes (avec union()) si la somme des volumes != enveloppe : il y a au moins un trou holes = envelop.difference(volumes) gestion du contenu des trous obtenus : traitement des différents cas : polygons, multipolygons, polylines ...

---

3ème temps: recherche des overlaps

Code:

---

1 2 3 4 5 6

pour chaque altitude pour chaque couple de volumes vol_i, vol_j si volumes concernés par l altitude inter = vol_i.intersection(vol_j) si intersection non vide gestion du contenu des intersections obtenues : traitement des différents cas : polygons, multipolygons, polylines ...

---

Voilà pour ce qui est de l'algorithmie.
Sur le fond, j'ai l'impression que ce que je dois faire est assez simple : union ou intersection, par rapport à toutes les fonctions de géométrie possibles. Mes polygones sont définis par des valeurs entières alors que la précision va jusqu'à 10**-12. Si j'ai rencontré ce cas, d'autres vont le rencontrer à coup sur.
Ce qui est bloquant, ce n'est pas les erreurs de précision qui font apparaitre des micro trous (car je pourrais les exclure en fonction d'une taille fournie par mon client), mais le fait que l'union ne se fait pas et qu'on part en erreur : je pense que ça relève d'un bug de shapely (ou geos). Cette erreur ne se reproduit plus pour l'instant car j'ai changé la façon de calculer l'enveloppe, mais j'ai toujours une union plus loin, donc on ne peut pas savoir si cela va se reproduire ou pas. Potentiellement c'est possible. Si ça se reproduit, je pourrais capter l'exception, mais pour en faire quoi ? si l'union n'est pas faite, comment faire autrement ?
Voilà ou j'en suis.
Cordialement