Discussion :

[awalter1]

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 :

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
		# 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]

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

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

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

[awalter1]

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 :

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

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 :

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

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

[awalter1]

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

[plxpy]

je regarde tout ça (voir tes MP)

[plxpy]

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 :

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

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

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.

[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