Discussion :

[awalter1]

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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

[awalter1]

J'utilise la version 3.3.0 de geos.

[plxpy]

Bonjour 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

[awalter1]

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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

[plxpy]

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

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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)

[fred1599]

chapeau!!!

[awalter1]

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

[awalter1]

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

[plxpy]

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.

[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

[awalter1]

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

[awalter1]

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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

[awalter1]

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 ...

[awalter1]

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

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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

[plxpy]

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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

[awalter1]

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

[awalter1]

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

[plxpy]

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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) :



En zoomant :




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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :




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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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.

[awalter1]

Bonjour,
Une analyse toujours aussi précise et complete .
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

[plxpy]

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :



et après :



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



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




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



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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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.