Discussion :

[awalter1]

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

Code :

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

1
2
3
4
polygon=Polygon(((0,0),(0,4),(4,4),(4,0),(0,0)))
print polygon.is_ring
 
False

Merci

[awalter1]

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

[plxpy]

Bonjour

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

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

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

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

doc constructeur de Polygon

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

et doc constructeur de LinearRing :

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

Code :

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

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

[awalter1]

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

[plxpy]

mais ce n'est pas très efficace.

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

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

[plxpy]

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

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

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
>>> from shapely.geometry import Polygon
>>> p1=Polygon(((0,0),(0,1),(1,1),(1,0)))
>>> p2=Polygon(((1,0),(1,1),(2,1),(2,0)))
>>>
>>> p1.intersects(p2)
True
>>> p1.touches(p2)
True
>>> 
>>> p3=Polygon(((0.5,0),(0.5,1),(1.5,1),(1.5,0)))
>>> p1.intersects(p3)
True
>>> p1.touches(p3)
False

La "logique" Shapely est héritée de la bibliothèque C/C++ GEOS, portage de JTS (Java Topology Suite, voir les spécifications techniques), elle-même calquée sur les standards OGC. Elle peut parfois sembler déroutante au premier abord mais, en fait, c'est la façon de traiter les cas "finauds" (tordus ?) qu'on rencontre plus souvent qu'on ne le pense (et je sais de quoi je parle !)

[awalter1]

oui, en creusant un peu on trouve ce qui convient.
au final je fais

Code :

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

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

[plxpy]

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

Code :

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

1
2
3
4
5
6
7
8
9
>>> from shapely.geometry import Polygon
>>> 
>>> polygone_en_U = Polygon(((0,0),(0,5),(1,5),(1,1),(2,1),(2,5),(3,5),(3,0)))
>>> polygone_horiz = Polygon(((0,3),(0,4),(4,4),(4,3)))
>>> 
>>> inter = polygone_en_U.intersection(polygone_horiz)
>>> 
>>> type(inter)
<class 'shapely.geometry.multipolygon.MultiPolygon'>

Mais tout dépend des polygones que tu traites. Peut-être que ce cas de figure ne se présente (mais, perso, j'ai appris à me méfier, dans ce domaine, des cas de figure qui n'arrivent "jamais" ...)

[awalter1]

oui, en creusant un peu on trouve ce qui convient.
au final je fais

Code :

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

1
2
3
4
typ_polygon=type(poly1)
inter = poly1.intersection(poly2)
if type(inter) == typ_polygon:
   print "c'est bien une intersection de surface"

Maintenant, comment détecter un polygone dont au moins 2 segments se croisent (à moins de tester tous les couples de segments soi même). A priori, ce n'est pas détecter lors de la création du polygone.

Merci

[awalter1]

j'ai trouvé : on fait poly1.is_valid

[plxpy]

j'ai trouvé : on fait poly1.is_valid

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

Par contre, ce que tu fais ne marche pas :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
>>> pu = polygone_en_U
>>> ph = polygone_horiz
>>> 
>>> pu.intersects(ph) and not pu.touches(ph)
True
>>> 
>>> typ_polygon=type(pu)
>>> inter = pu.intersection(ph)
>>> if type(inter) == typ_polygon:
...     print "c'est bien une intersection de surface"
... 
>>>

et mes deux polygones sont tout ce qu'il y a de plus valides :

Code :

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

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

[awalter1]

J'ai une intersection qui me renvoie une linestring et un polygon :

Code :

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

1
2
print ">>> ",poly1.intersection(poly2)
>>> GEOMETRYCOLLECTION (LINESTRING (213551.0 19113.0 ....), POLYGON ((213878.0 19114.0 ....)))

Comment puis je extraire ces deux entités ?
Merci

[plxpy]

A partir du moment où tu as une "multi-géométrie" (MultiPoint, MultiLineString, MultiPolygon et, à fortiori, une GeometryCollection), les géométries sont accessibles via l'attribut geoms

[awalter1]

Bonjour,
Je vois que si je fais

Code :

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

1
2
p=Polygon([[0,0][0,5],[5,5][5,0]])
print p.geoms

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

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

[plxpy]

Shapely est il bien déterministe et fiable ?

Sans aucune hésitation : oui !

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

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

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

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

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

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

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

[plxpy]

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

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

Code :

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
>>> from shapely.geometry import Point, LineString, Polygon, GeometryCollection
>>>
>>> for couleur in polygons:
...     inter = vert.intersection(polygons[couleur])
...     if inter.is_empty:
...             print "pas d intersection avec le polygone", couleur
...     elif isinstance(inter,Point):
...             print "intersection ponctuelle avec le polygone", couleur
...     elif isinstance(inter,LineString):
...             print "intersection lineaire avec le polygone", couleur
...     elif isinstance(inter,Polygon):
...             print "intersection surfacique avec le polygone", couleur
...     elif isinstance(inter,GeometryCollection):
...             print "intersection compliquee avec le polygone", couleur
...
intersection ponctuelle avec le polygone bleu
intersection surfacique avec le polygone orange
pas d intersection avec le polygone violet
intersection compliquee avec le polygone marron
intersection lineaire avec le polygone rouge
>>>
>>>>>> vert_marron = vert.intersection(polygons['marron'])
>>>
>>> for geom in vert_marron.geoms:
...     if isinstance(geom,Point):
...             print "un point"
...     elif isinstance(geom, LineString):
...             print "un arc"
...     elif isinstance(geom,Polygon):
...             print "un polygone"
...     else:
...             print "autre chose"
...
un point
un arc
un arc
un polygone
>>>

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

[plxpy]

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

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

[awalter1]

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

Code :

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
# l_polygon contains a list of polygons
#  l_polygon[i][0] = polygon name
#  l_polygon[i][1] = upper level
#  l_polygon[i][2] = lower level
#  l_polygon[i][3] = list of points coordinates
#
for l in range(1,len(l_level)):
       selected_level = l_level[l]
 
       for i in range(len(l_polygon)):
               if selected_level <= l_polygon[i][1] and selected_level > l_polygon[i][2]:
                       for j in range(i+1,len(l_polygon)):
                               if selected_level <= l_polygon[j][1] and selected_level > l_polygon[j][2]:
                                       inter = l_polygon[i][3].intersection(l_polygon[j][3])
                                       # if intersection
                                       if inter.is_empty == False:
                                               l_inter = []
                                               # check if collection or polygon, exclude point and linestring
                                               if isinstance(inter,Polygon):
                                                       l_inter.append(inter)
                                               elif isinstance(inter,GeometryCollection):
                                                       for geom in inter.geoms:
                                                               if isinstance(geom,Polygon) == True:
                                                                       l_inter.append(geom)
                                               for poly in l_inter:
                                                       print "Intersection between volumes <%s> and <%s> between level
%d and %d" % (l_polygon[i][0],l_polygon[j][0],min_level,selected_level)
                                                       # find the points name corresponding to the overlapping polygon
                                                       l_points_overlap = []
                                                       point_overlap = ""
                                                       for coord in list(geom.exterior.coords)[:-1]:
                                                               for pt in l_points:
                                                                       if (str(coord[0]) == "%s.0" % pt[1]) and (str(coord[1]) ==
"%s.0" % pt[2]):
                                                                               if pt[0] not in l_points_overlap:
                                                                                       point_overlap = pt[0]
                                                               # the point name is not found : get only its coordinates
                                                               if point_overlap == "":
                                                                       l_points_overlap.append((self.ApplicationConvertNumlatlong(coord[0],coord[1])))
                                                               else:
                                                                       l_points_overlap.append(point_overlap)
 
                                                       ch_points_overlap = ""
                                                       for pt in l_points_overlap:
                                                               if type(pt) == type(""):
                                                                       ch_points_overlap = "%s %s" % (ch_points_overlap,pt)
                                                               else:
                                                                       ch_points_overlap = "%s %s/%s" % (ch_points_overlap,pt[0],pt[1])
 
                                                       print "   intersection definition : %s" % (ch_points_overlap)

[plxpy]

Quelques petites remarques

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

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

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

Sinon, tes "levels" correspondent aux différents z ?

[awalter1]

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