Python: le module Shapely, géométries, prédicats spatiaux, analyse spatiale, matrices de Clementini (DE-9IM), traitements de shapefiles ou autres, et représentation avec matplotlib

Parmi les normes de l'Open Geospatial Consortium, ou OGC figure Simple Features.  Celle-ci définit le modèle géométrique à utiliser pour les entités vectorielles SIG ( Simple Feature Access - Part 1: Common Architecture: www.opengeospatial.org/standards/sfa, Simple Feature Access - Part 2: SQL Option: www.opengeospatial.org/standards/sfs) sous la forme de:

• un modèle de géométrie 2D pour les objets vectoriels (points, lignes, polygones...);
• le traitement des projections avec des identifiants de systèmes de référence spatiale ou SRID;  
• un format de données "well-known" (WKT, WKB) pour le stockage de ces objets;
• les opérations spatiales (intersections, buffers, etc.) ;
• un schéma de base de données spatiale et un modèle d'accès.

Quelques exemples de représentation de géométries OGC au format WKT (sans SRID):
POINT (180 270)
LINESTRING (250 370, 180 330, 180 190)
POLYGON ((180 370, 130 310, 160 230, 160 230, 250 260, 180 370))
LINEARRING (180 370, 130 310, 160 230, 160 230, 250 260, 180 370)

Il est aussi nécessaire de bien comprendre les concepts d'intérieur, de limite et d'extérieur d'un objet géométrique, ainsi que sa dimension:

• l'intérieur d'une géométrie se compose de toute la géométrie à l'exception des limites 
•  la limite d'une géométrie est une géométrie de dimension inférieure
•  l'extérieur d'une géométrie se compose de tous les points qui ne sont ni dans l'intérieur ni dans la limite;

Vu la composition de l’OGC (ESRI, Oracle et al.), cette norme devrait, en théorie, être respectée par tous les grands acteurs du SIG, mais ce n'est pas toujours le cas, comme toujours... (voir benjamin.chartier.free.fr/pro/). Dans le monde du libre, par contre, c'est le modèle utilisé, entre autres, par PostGIS qui se base sur la librairie GEOS, un élément de la Java Topology Suite (JTS).

En Java, il est possible d'utiliser cette Java Topology Suite seule grâce au programme JTS Test Builder (voir www.perrygeo.net/wordpress/ ) qui permet d'« explorer » les géométries OGC et les opérations spatiales.

En Python, Sean Gillies a créé le module Shapely pour permettre d'effectuer, en simplifiant, tout ce qu'il est possible de faire avec PostGIS, sans utiliser SQL ni Java, au niveau des géométries. Il est aussi basé sur la librairie GEOS, qui doit être installée auparavant (pour Windows, des distributions complètes de Shapely sont disponibles).  Le module ne traite que les géométries dans un plan cartésien xy, sans projection (SRID). Sean Gillies part du principe qu'il vaut mieux travailler en coordonnées cartésiennes et une fois les démarches effectuées, les projeter ou les reprojeter avec d'autres outils. Je développerai dans la suite les points suivants:

1. Les objets géométriques (points, lignes, polygones, multipoints, multilignes, multipolygones, multigéométriques);
2. Les "constructeurs" ou comment créer des nouvelles géométries à partir d'un objet géométrique existant (buffer...);
3. Les requêtes spatiales binaires (True, False) et les prédicats;
4. les requêtes spatiales D9-IM de l'OGC (matrice de Clementini);
5. L'analyse spatiale ou comment créer de nouveaux objets à partir de deux objets géométriques existants (union, intersection...) avec application à la matrice de Clementini;
6. Les applications pratiques (shapefiles, PostGIS, projections, autre traitement mathématique)
7. Comment dessiner les géométries avec Matplotlib et ma classe Plot_Shapely;
8. Conclusions

Les objets géométriques (points, lignes, polygones, multipoints, multilignes, multipolygones, multigéométriques)

Les principaux objets géométriques (classes) disponibles dans Shapely sont illustrés par le diagramme suivant créé avec Yuml (yuml.me/edit/25878f1a ):

La création d'un objet géométrique se fait en important sa classe géométrie:

>>> from shapely.geometry import Polygon
>>> poly2 = Polygon([(1, 2), (2, 5), (6, 5),(7,2),(5 1),(1 2)])
>>> poly2.wkt
'POLYGON ((1 2, 2 5, 6 5, 7 2, 5 1, 1 2))'

Shapely permet de visualiser l'objet au format WKT ou d'importer un objet défini en WKT. C'est cette forme qui sera utilisée dans la suite pour faciliter la compréhension.

>>> from shapely.wkt import loads
>>> poly1 = loads('POLYGON ((1 2, 2 5, 6 5, 7 2, 5 1, 1 2))')

Toutes les figures ont été faites avec le module matplotlib. Nous détaillerons la démarche par la suite.

Un point est construit à partir de ses coordonnées xy (z). Son aire est par définition nulle. L'objet MultiPoints est une collection de points

# point seul
 
>>> point1 = loads('POINT (4 5)')
 
# multipoints
 
>>> multipt = loads('MULTIPOINT ((2 2), (3 3), (3 4))')

Un LineString est une ligne, définie par l'OGC comme une séquence ordonnée de deux à plusieurs Points. Un MultiLineString est une collection de lignes (LineString)

# ligne simple et quelques propriétés
>>> ligne1 = loads('LINESTRING (3 1, 4 4, 5 5, 5 6)')
>>> ligne1.bounds
(3.0, 1.0, 5.0, 6.0)
>>> ligne1.length
5.5764912225414749
 
# multiligne
>>> multiligne = loads('MULTILINESTRING ((1 2, 1 2, 2 1, 3 2),(2 4, 3 4, 5 3))')

# les coordonnées des points constitutifs de la ligne
>>> list(ligne1.coords)  
[(3.0, 1.0), (4.0, 4.0), (5.0, 5.0), (5.0, 6.0)] 
# les noeuds limites (boundary)                
>>> ligne1.boundary.wkt
'MULTIPOINT (3 1, 5 6)' 

Un polygone s.s. peut être défini comme une aire simple ou en « anneau » (« troué »). Il est défini par une ou des séquences ordonnées de points.

# polygone simple et quelques propriétés
 
>>> poly1 = loads('POLYGON ((1 2, 2 5, 6 5, 7 2, 5 1, 1 2))')
>>> poly1.area
18.0
# limites (bounds)
>>> poly1.bounds
(1.0, 1.0, 7.0, 5.0)
 
# polygone troué et quelques propriétés
 
>>> poly2 = loads('POLYGON ((1 2, 2 5, 6 5, 7 2, 5 1, 1 2), (3 4, 2 2, 5 2, 6 4, 3 4))')
>>> poly2.is_valid
True
 
# création d'un polygone par buffer
 
>>> point1 = loads('POINT (4 5)')
>>> poly2 = point1.buffer(2)
poly2.wkt
'POLYGON ((6.00 5.00,[...]))'
>>> line1 = loads('LINESTRING (3 1, 4 4, 5 5, 5 6)')
>>> poly3 = line1.buffer(2)
'POLYGON ((2.1026334038989725 4.6324555320336760,,[...]))' 
 
# création avec une liste de points (MultiPoints) et la fonction convex_hull
>>> multi = loads('MULTIPOINT (2 5, 5 1, 1 2, 7 2, 1 2, 6 5)')
>>> poly = multi.convex_hull
>>> poly.wkt
'POLYGON ((5 1, 1 2, 2 5, 6 5, 7 2, 5 1))'

Les limites extérieures et intérieures d'un polygone sont des lignes fermées. C'est pourquoi le concept de LinearRing a été introduit par l'OGC pour les différencier des LineString.  Dans Shapely, Sean Gillies en fait bien une sous-classe de Polygon. 

>>> poly1.exterior.wkt
'LINEARRING (1 2, 2 5, 6 5, 7 2, 5 1, 1 2)'
 
>>> poly1.boundary.wkt
'LINESTRING (1 2, 2 5, 6 5, 7 2, 5 1, 1 2)'
 
# prédicat unaire
>>> poly2.exterior.is_ring
True
>>> poly2.boundary.is_ring
False

C'est à partir de ces LinearRing que les coordonnées des séquences de points peuvent être retrouvées

>>> list(poly1.exterior.coords)
[(1.0, 2.0), (2.0, 5.0), (6.0, 5.0), (7.0, 2.0), (5.0, 1.0), (1.0, 2.0)]
>>> list(poly2.exterior.coords)
[(1.0, 2.0), (2.0, 5.0), (6.0, 5.0), (7.0, 2.0), (5.0, 1.0), (1.0, 2.0)]
>>> list(poly2.interiors[0].coords)
[(3.0, 4.0), (2.0, 2.0), (5.0, 2.0), (6.0, 4.0), (3.0, 4.0)]

et un MultiPolygon qui est une collection de polygones

multipoly = loads('MULTIPOLYGON (((2 4, 1 3, 2 2, 4 3, 2 4)),((5 2, 4 2, 4 1, 4 1, 5 1, 5 2)))')

Le résultat des diverses opérations peut conduire à l'obtention d'objets hétérogènes associant diverses géométries. Cette classe de l'OGC est aussi implantée dans Shapely: GeometryCollection

GEOMETRYCOLLECTION (LINESTRING (210 380, 202.45283018867926 334.7169811320755),
  LINESTRING (235 283.57142857142856, 260 280, 250 370, 210 380),
  POLYGON ((235 283.57142857142856, 250 260, 160 230, 130 310, 180 370, 202.45283018867926 334.7169811320755, 235 283.57142857142856)))

Les "constructeurs" ou comment créer des nouvelles géométries à partir d'un objet géométrique existant (buffer, ...) 

Shapely possède les contructeurs suivants:

• objet.buffer(distance, [résolution]): permet de créer des buffers positifs ou négatifs. Le terme résolution peut être utilisé pour préciser le nombre de segments de la géométrie résultante (voir script "polygone");
• objet.convex.hull (voir script "polygone")
• objet.enveloppe
• objet.simplify
• objet.boundary (voir script linearring - polygone)
• (objet.exterior ( réservé aux polygones, voir script linearring - polygone))
• objet.centroid

Les requêtes spatiales binaires (True, False) et les prédicats

Les relations spatiales entre géométries sont aussi normalisées par l'OGC. Elles se traduisent dans Shapely par des prédicats (fonctions renvoyant une valeur binaire, True ou False). Ils sont amplement expliqués dans le manuel de Shapely.

1)  prédicats unaires, car ils s'appliquent à un seul objet géométrique

• objet.has_z¶
• objet.is_empty
• objet.is_ring  (voir script "linearring -polygone")
• objet.is_simple 
• objet.is_valid (voir script "polygones")

2) prédicats binaires, car ils s’appliquent à deux objets géométriques

• objet1.crosses(objet2) signifie que l'intérieur de objet1 recoupe l'intérieur de objet2, que objet1 ne contient pas objet2 et dont la dimension résultante doit être inférieure à celle des 2 objets;
• objet1.intersects(objet2) signifie qu'il y a une intersection entre les limites et les intérieurs des deux objets;
• objet1.touches(objet2) signifie que la limite de objet1 recoupe uniquement celle de objet2 mais qu'il n'y a pas d'intersection entre leurs intérieurs;
• objet1.within(objet2) signifie que l'intérieur et la limite de objet1 recoupent uniquement l'intérieur de objet2 (pas la limite);
• objet1.contains(objet2) signifie que l'intérieur et la limite de objet2 recoupe uniquement l'intérieur de objet1 et que leurs limites ne se touchent (touches)
• objet1.disjoint(objet2)
• objet1.equals(objet2)
• objet1.almost_equals(objet2) 

Il en résulte, par exemple, qu'un polygone ne contient pas sa limite...(benjamin.chartier.free.fr/pro/)

Les requêtes spatiales D9-IM de l'OGC (matrice de Clementini)

Les résultats sont envoyés sous forme de matrice. C'est la matrice de Clementini ou notation DE-9IM (Dimensionally Extended Nine-Intersection Model de l' OGC, voir benjamin.chartier.free.fr/pro/ ,  postgis.org/documentation/manual-svn/ch04.html#DE-9IM )

• objet1.relate(objet2)

Nous allons construire cette matrice dans le point suivant. Notons que Sean Gillies a aussi développé un module d9im pour traiter ces relations.

L'analyse spatiale ou comment créer de nouveaux objets à partir de deux objets géométriques existants (union, intersection...) avec application à la matrice de Clementini

Shapely fournit aussi les méthodes permettant de visualiser ces relations spatiales en créant de nouvelles géométries (résultats de l'union, de l'intersection, de la différence ou de la différence symétrique entre 2 objets)

• object1.intersection(objet2)
• objet1.difference(objet2)
• objet1.symmetric_difference(objet2)
• objet1.union(objet2)

Pour illustrer toutes ces commandes, nous allons construire la matrice de Clementini pour caractériser l'intersection entre les 2 polygones suivants:

Le résultat obtenu est

>>> polygon1.relate(polygon2)
212101212
#ou matrice
# 2 1 2
# 1 0 1
# 2 1 2

Mais qu'est ce que cela signifie ? En fait, les chiffres obtenus correspondent à des dimensions d'objets. La matrice est en effet construite de la manière suivante:

                                             Intérieur(Polygone 2)    Limite(Polygone 2)    Extérieur(Polygone 2)
        ______________________________________________________________________
       Intérieur(Polygone 1)     dimension(inters.)       dimension(inters.).         dimension(inters.)
       Limite(Polygone 1)                 "                                    "                                      "
       Extérieur(Polygone 1)            "                                    "                                      "    

Cette matrice est construite ici avec les méthodes d'analyse de Shapely:

                       dim =   2                                            dim = 1                                  dim = 2

                 dim =   1                                           dim = 0                                  dim = 1

                  dim =   2                                            dim = 1                                dim=2

 Le résultat est donc:

2 1 2
1 0 1
2 1 2

ou 212101212 conforme à ce qui a été obtenu par polygon1.relate(polygon2)

Quel est l'intérêt d'une telle démarche ? Celle de pouvoir créer des filtres spatiaux très précis pour rechercher des relations impossibles à obtenir par les méthodes traditionnelles. Par exemple, si je veux obtenir toutes les intersections de type "ligne dans polygone" (ligne 1, colonne 2 et ligne 2, colonne 1 de la matrice), je vais utiliser un filtre du genre  "¨*1*1***** ". Cette méthode et surtout utilisée pour les les recherches dans les bases de données spatiales:

-- Identify road segments that cross on a line
SELECT a.id
FROM roads a, roads b
WHERE a.id != b.id
AND a.geom && b.geom
AND ST_Relate(a.geom, b.geom, '1*1***1**');
(tiré de  postgis.org/documentation/manual-svn/ch04.html#DE-9IM )

Chritian Strobl explique tout cela dans gis.hsr.ch/wiki/images/3/3d/9dem_springer.pdf avec des matrices pour de nombreux cas.

Les applications pratiques (shapefiles, PostGIS, projections , autre traitement mathématique)

Tout ça, c'est bien beau, me direz-vous, mais cela reste apparemment purement théorique. Et non, car en pratique Shapely peut être utilisé pour traiter beaucoup de choses, des shapefiles, des résultats obtenus par des requêtes spatiales provenant des bases de données ou des traitements mathématiques divers.

En pratique, la démarche consiste toujours à transformer les objets en géométries Shapely pour pouvoir faire ensuite les traitements.

Pour un shapefile, la démarche a déjà été partiellement expliquée partiellement dans www.portailsig.org/content/python-visualiser-et-traiter-des-donnees-spatiales-de-type-xyz-shapefiles-ou-mnt-srtm. Sachant que la commande geom_type permet d'obtenir le type de géométrie d'un objet:

>>> poly1.geom_type
'Polygon'
>>> ligne1.geom_type
'LineString'
etc.

on peut effectuer le traitement avec le module ogr:

>>> from osgeo import ogr
>>> from shapely.wkb import loads
>>> source = ogr.Open("testpoly.shp")
>>> couche = source.GetLayerByName("testpoly")
>>> for element in couche:
...     geom = loads(element.GetGeometryRef().ExportToWkb())
...     if geom.geom_type == 'Point':
...        print geom.type
...        print geom
...     if geom.geom_type == 'LineString':   
...        print geom.type
...        print geom
...     if geom.geom_type == 'MultiLineString': 
...        print geom.type
...        print geom
...     if geom.geom_type == 'MultiPolygon':
...        print geom.type
...        print geom
...     if geom.geom_type == 'Polygon':
...        print geom.type
...        print geom
... 
Polygon
POLYGON ((0.0909447004608295 0.8075576036866359, 0.1416359447004608 0.8014746543778801, 0.3606221198156682 0.7021198156682027, 0.2409907834101383 0.5480184331797234, 0.0868894009216590 0.5845161290322580, 0.0544470046082949 0.7426728110599077, 0.0909447004608295 0.8075576036866359))
Polygon
POLYGON ((0.2754608294930876 0.7933640552995391, 0.5370276497695853 0.8136405529953916, 0.4173963133640554 0.5926267281105990, 0.2267972350230415 0.6777880184331797, 0.2267972350230415 0.7447004608294931, 0.2754608294930876 0.7933640552995391))

ou avec le module shpUtils (voir www.portailsig.org/content/python-et-les-shapefiles et www.aaronland.info/weblog/2009/04/12/hammock/):

>>> import shpUtils
>>> from shapely.geometry import Polygon
>>> shp = 'testpoly.shp'
>>> polys = []
>>> for element in shpUtils.loadShapefile(shp):
...        for part in element['shp_data']['parts'] :
...            poly=[]
...            for pt in part['points'] :
...                if pt.has_key('x') and pt.has_key('y') :
...                    poly.append((pt['x'], pt['y']))
...            poly = tuple(poly)
...            p = Polygon(poly)
...            polys.append(p)
 
>>> for pol in polys:
...     print pol.wkt
... 
POLYGON ((0.0909447004608295 0.8075576036866359, 0.1416359447004608 0.8014746543778801, 0.3606221198156682 0.7021198156682027, 0.2409907834101383 0.5480184331797234, 0.0868894009216590 0.5845161290322580, 0.0544470046082949 0.7426728110599077, 0.0909447004608295 0.8075576036866359))
POLYGON ((0.2754608294930876 0.7933640552995391, 0.5370276497695853 0.8136405529953916, 0.4173963133640554 0.5926267281105990, 0.2267972350230415 0.6777880184331797, 0.2267972350230415 0.7447004608294931, 0.2754608294930876 0.7933640552995391))

Avec PostGIS les démarches ont été expliquées dans sgillies.net/blog/531/the-shapely-alchemist/ ou ici avec une adaptation du script "requête SELECT avec SQLObject" de www.portailsig.org/content/python-les-bases-de-donnees-geospatiales-2-mapping-objet-relationnel-orm-sqlalchemy-sqlobjec

résultats obtenus directement avec PostGIS en SQL:

>>> polys = testpoly._connection.queryAll("""SELECT Srid(the_geom), AsText(the_geom) FROM testpoly""")
>>> polys
[(4326, 'POLYGON((0.09094470046083 0.807557603686636,0.141635944700461 0.80147465437788,0.360622119815668 0.702119815668203,0.240990783410138 0.548018433179723,0.086889400921659 0.584516129032258,0.054447004608295 0.742672811059908,0.09094470046083 0.807557603686636))'),[...]')]
>>> polys[0][1]
'POLYGON((0.09094470046083 0.807557603686636,0.141635944700461 0.80147465437788,0.360622119815668 0.702119815668203,0.240990783410138 0.548018433179723,0.086889400921659 0.584516129032258,0.054447004608295 0.742672811059908,0.09094470046083 0.807557603686636))'
# et donc
b = loads(polys[0][1]]

résultats obtenus avec SQLObject et Shapely:

>>> from sqlobject import *
>>> from shapely.wkb import loads
>>> db='postgres://localhost:5432/testpostgis'
>>> connection = connectionForURI(db)
>>> class testpoly(SQLObject): 
...     _connection = connection
...     _fromDatabase = True
...     the_geom = StringCol()
...     def _get_geom(self):
...         value = self._SO_get_the_geom()
...         test = loads(value.decode('hex')).wkt
...         return test
... 
>>> a = testpoly.get(1)
>>> a.geom
'POLYGON ((0.0909447004608300 0.8075576036866360, 0.1416359447004610 0.8014746543778800, 0.3606221198156680 0.7021198156682030, 0.2409907834101380 0.5480184331797230, 0.0868894009216590 0.5845161290322580, 0.0544470046082950 0.7426728110599080, 0.0909447004608300 0.8075576036866360))'
>>> 

Le SRID n'est pas pris en compte avec Shapely, mais hackmap.blogspot.com/2008/03/ogr-python-projection.html montre comment utiliser le module ogr pour le faire:

>>> from shapely.geometry import LineString
>>> from shapely.wkb import loads
>>> from osgeo import ogr
>>> def project(geom, from_epsg, to_epsg):
...      to_srs = ogr.osr.SpatialReference()
...      to_srs.ImportFromEPSG(to_epsg)
...      from_srs = ogr.osr.SpatialReference()
...      from_srs.ImportFromEPSG(from_epsg)
...      ogr_geom = ogr.CreateGeometryFromWkb(geom.wkb)
...      ogr_geom.AssignSpatialReference(from_srs)
...      ogr_geom.TransformTo(to_srs)
...      return loads(ogr_geom.ExportToWkb())
>>> ligne = LineString([[4,50], [5,51]])
>>> ligne.wkt
LINESTRING (4 50, 5 51)
>>> ligneproj = project(ligne, from_epsg=4326, to_epsg=31370)
>>> ligneproj.wkt
'LINESTRING (123742.1348178901971551 76638.2337631648406386, 194490.5125554261030629 188001.2492922432720661)'

Le manuel et la distribution fournissent de nombreux autres exemples comme la procédure pour découper des lignes en tronçons suivant une distance déterminée (gispython.org/shapely/docs/1.2/manual.html#linear-referencing-methods).

Mais Shapely peut aussi être utilisé à toute autre chose que le domaine spatial comme l'illustre cette réponse de plxpy sur www.developpez.net/forums/d755595/autres-langages/python-zope/general-python/points-intersection-courbes/ suite à une question pratique:

"Une petite démo pour calculer les intersections entre sinus et cosinus dans [0,2*pi]" (script original légèrement modifié)

>>> import math
>>> from shapely.geometry import LineString
>>> list_x = [ 2.*math.pi/5000*i for i in range(5001) ]
>>> courbe_sinus = [ (list_x[i],math.sin(list_x[i])) for i in range(5001) ]
>>> courbe_cosinus = [ (list_x[i],math.cos(list_x[i])) for i in range(5001) ]
>>> sinus = LineString(courbe_sinus)
>>> cosinus = LineString(courbe_cosinus)
>>> intersections = sinus.intersection(cosinus)
>>> intersections.geom_type
'MultiPoint'
>>> intersections.wkt
'MULTIPOINT (0.7853981633974484 0.7071067811865475, 3.9269908169872414 -0.7071067811865476)'

Comment dessiner les géométries avec matplotlib (classe Plot_Shapely)

En pratique, Shapely utilise des tableaux (array) numpy pour représenter les objets:

>>> line1.xy
(array('d', [3.0, 4.0, 5.0, 5.0]), array('d', [1.0, 4.0, 5.0, 6.0]))

Ces arrays numpy sont directement utilisables par matplotlib:

import pylab
from shapely.wkt import loads
 
def plot_coords(ax, ob):
     x, y = ob.xy
     ax.plot(x, y, 'o', color='r', ms=20) 
    # 'o', rond, 'r', couleur rouge, 20, taille)
def plot_line(ax, ob):
    x, y = ob.xy
    ax.plot(x, y, color='b')
def plot_multipt(ax, ob):
    for pt in ob:
        plot_coords(ax,pt)
def plot_multilignes(ax, ob):       
    for ligne in ob:
        plot_line(ax,ligne)
 
ax = pylab.gca() #ouverture d'une "fenêtre" pylab
 
ligne = loads('LINESTRING (3 1, 4 4, 5 5, 5 6)')
point  = loads('POINT (4 5)')
# ensuite suivant les cas 
plot_coords(ax,point)  #-> point
plot_line(ax, ligne)   #-> ligne
plot_coords(ax, ligne) #-> points
 
pylab.show() 

Cela marche aussi pour les limites des polygones, mais pas pour les aires (hormis en utilisant la commande pylab.fill). Matplotlib offre d'autres commandes pour figurer des surfaces, bien illustrées dans sgillies.net/blog/1013/painting-punctured-polygons-with-matplotlib/. Sur cette base Seans Gillies a développé un module spécifique pour le faire, descartes.

from descartes.patch import PolygonPatch
 
def plot_poly(ax,ob):
    patch = PolygonPatch(ob, facecolor='#6699cc', edgecolor='#6699cc')
    ax.add_patch(patch)
 
def plot_multipol(ax,ob):
    for polygon in ob:
        plot_poly(ax,polygon)
poly = loads('POLYGON ((1 2, 2 5, 6 5, 7 2, 5 1, 1 2))')
multipoly = loads('MULTIPOLYGON (((2 4, 1 3, 2 2, 4 3, 2 4)),((5 2, 4 2, 4 1, 4 1, 5 1, 5 2)))')
 
plot_poly(ax,poly)
plot_multipol(ax, multipoly)

Je me suis donc créé sur ces bases une classe Plot_shapely qui me permet de dessiner toutes les géométries (diagramme obtenu avec Yuml - yuml.me/edit/37c88fdd )

#!/usr/bin/env python
# encoding: utf-8
"""
Plot_shapely.py
Martin Laloux 2010
"""
import pylab
from shapely.wkt import loads
from descartes.patch import PolygonPatch
class Plot_shapely(object):
   """dessin de géométries Shapely avec matplotlib - pylab"""
 
   def __init__(self,obj,ax, coul=None, alph=1):
       """
        Paramètres
        ----------------------------------------------------------------
        ax de pylab, obj = objet géometrique, coul = couleur matplotlib, alph = transparence 
 
        Exemple
        --------------------------------------------------------------------
        >>> from shapely.wkt import loads
        >>> ax = pylab.gca()
        >>> ligne = loads('LINESTRING (3 1, 4 4, 5 5, 5 6)')
        >>> a  = Plot_shapely(ligne,ax,'r', 0.5)
        >>> a.plot
        ou directement
        >>> Plot_shapely(ligne,ax,'#FFEC00'').plot
        >>> pylab.show()
       """
       self.obj = obj
       self.type = obj.geom_type
       self.ax = ax
       self.coul= coul
       self.alph=alph
 
   def plot_coords(self):
       """points"""
       x, y = self.obj.xy
       self.ax.plot(x, y, 'o', color=self.coul) 
 
   def plot_ligne(self):
       """lignes"""
       x, y = self.obj.xy
       self.ax.plot(x, y, color=self.coul, alpha=self.alph, linewidth=3)
 
   def plot_polygone(self):
       """polygones"""
       patch = PolygonPatch(self.obj, facecolor=self.coul,edgecolor=self.coul, alpha=self.alph)
       self.ax.add_patch(patch)
 
   def plot_multi(self):
       """multipoints, multilignes,multipolygones + GeometryCollection"""
       for elem in self.obj:
           Plot_shapely(elem, self.ax, self.coul,self.alph).plot
 
 
   @property
   def plot(self):
       """dessine en fonction du type géometrique"""
       if self.type == 'Point':
           self.plot_coords()
       elif self.type == 'Polygon':
           self.plot_polygone()
       elif self.type == 'LineString':   
           self.plot_ligne()
       elif "Multi" in self.type:
           """ex. MultiPolygon"""
           self.plot_multi()
       elif self.type == 'GeometryCollection':
           self.plot_multi()
       elif self.type == 'LinearRing':
           self.plot_line()
       else:
           raise ValueError("inconnu au bataillon: %s" % self.type)
 
if __name__ == '__main__':
   import doctest
   doctest.testmod()

Conclusions

Shapely est, pour moi, un magnifique outil pour celui qui aime Python et veut comprendre ce que l'on fait réellement avec un SIG ou une base de données spatiale qui respectent les normes de l' OGC.
Bien sur, il ne gère pas les données alphanumériques associées aux géométries, de même que les règles topologiques propres à certains logiciels. Ce devrait être possible en lui associant des modules comme ogr mais, cela reste purement théorique et ce n'est pas le but du module.

J'espère que cette introduction vous donnera les clés pour « jouer  » avec Python et Shapely, ce qui reste pour moi le plus grand plaisir.  Sean Gillies montre surtout comment une bonne connaissance de Python permet d'optimiser n'importe quel script pour n'importe quel logiciel.

Tous les traitements ont été effectués sur Mac OS X avec Python 2.5.4, diverses versions de Shapely suivant l'évolution du module (dernière 1.2), matplotlib, descartes et autres modules Python, Yuml, directement à partir de Python,  et Inkscape pour la seule figure non matplotlib.