Skip to Content

Nous adressons toutes nos pensées à la famille de notre ami Jérôme !

http://www.forumsig.org/showthread.php/43488-Disparition-de-Phoenix

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


python

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

  1. >>> from shapely.geometry import Polygon
  2. >>> poly2 = Polygon([(1, 2), (2, 5), (6, 5),(7,2),(5 1),(1 2)])
  3. >>> poly2.wkt
  4. '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.

wkt

  1. >>> from shapely.wkt import loads
  2. >>> 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

 

points - multipoint

  1. # point seul
  2.  
  3. >>> point1 = loads('POINT (4 5)')
  4.  
  5. # multipoints
  6.  
  7. >>> 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)

 

lignes - multilignes

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

traitement lignes

  1. # les coordonnées des points constitutifs de la ligne
  2. >>> list(ligne1.coords)
  3. [(3.0, 1.0), (4.0, 4.0), (5.0, 5.0), (5.0, 6.0)]
  4. # les noeuds limites (boundary)
  5. >>> ligne1.boundary.wkt
  6. '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.

polygones

  1. # polygone simple et quelques propriétés
  2.  
  3. >>> poly1 = loads('POLYGON ((1 2, 2 5, 6 5, 7 2, 5 1, 1 2))')
  4. >>> poly1.area
  5. 18.0
  6. # limites (bounds)
  7. >>> poly1.bounds
  8. (1.0, 1.0, 7.0, 5.0)
  9.  
  10. # polygone troué et quelques propriétés
  11.  
  12. >>> 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))')
  13. >>> poly2.is_valid
  14. True
  15.  
  16. # création d'un polygone par buffer
  17.  
  18. >>> point1 = loads('POINT (4 5)')
  19. >>> poly2 = point1.buffer(2)
  20. poly2.wkt
  21. 'POLYGON ((6.00 5.00,[...]))'
  22. >>> line1 = loads('LINESTRING (3 1, 4 4, 5 5, 5 6)')
  23. >>> poly3 = line1.buffer(2)
  24. 'POLYGON ((2.1026334038989725 4.6324555320336760,,[...]))'
  25.  
  26. # création avec une liste de points (MultiPoints) et la fonction convex_hull
  27. >>> multi = loads('MULTIPOINT (2 5, 5 1, 1 2, 7 2, 1 2, 6 5)')
  28. >>> poly = multi.convex_hull
  29. >>> poly.wkt
  30. '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. 

 

linearring - polygone

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

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

 

coord. points

  1. >>> list(poly1.exterior.coords)
  2. [(1.0, 2.0), (2.0, 5.0), (6.0, 5.0), (7.0, 2.0), (5.0, 1.0), (1.0, 2.0)]
  3. >>> list(poly2.exterior.coords)
  4. [(1.0, 2.0), (2.0, 5.0), (6.0, 5.0), (7.0, 2.0), (5.0, 1.0), (1.0, 2.0)]
  5. >>> list(poly2.interiors[0].coords)
  6. [(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

multipolygone

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

d9im

  1. >>> polygon1.relate(polygon2)
  2. 212101212
  3. #ou matrice
  4. # 2 1 2
  5. # 1 0 1
  6. # 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:

geom_type

  1. >>> poly1.geom_type
  2. 'Polygon'
  3. >>> ligne1.geom_type
  4. 'LineString'
  5. etc.

on peut effectuer le traitement avec le module ogr:

 

Shapely-ogr

  1. >>> from osgeo import ogr
  2. >>> from shapely.wkb import loads
  3. >>> source = ogr.Open("testpoly.shp")
  4. >>> couche = source.GetLayerByName("testpoly")
  5. >>> for element in couche:
  6. ... geom = loads(element.GetGeometryRef().ExportToWkb())
  7. ... if geom.geom_type == 'Point':
  8. ... print geom.type
  9. ... print geom
  10. ... if geom.geom_type == 'LineString':
  11. ... print geom.type
  12. ... print geom
  13. ... if geom.geom_type == 'MultiLineString':
  14. ... print geom.type
  15. ... print geom
  16. ... if geom.geom_type == 'MultiPolygon':
  17. ... print geom.type
  18. ... print geom
  19. ... if geom.geom_type == 'Polygon':
  20. ... print geom.type
  21. ... print geom
  22. ...
  23. Polygon
  24. 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))
  25. Polygon
  26. 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/):

 

Shapely-shpUtils

  1. >>> import shpUtils
  2. >>> from shapely.geometry import Polygon
  3. >>> shp = 'testpoly.shp'
  4. >>> polys = []
  5. >>> for element in shpUtils.loadShapefile(shp):
  6. ... for part in element['shp_data']['parts'] :
  7. ... poly=[]
  8. ... for pt in part['points'] :
  9. ... if pt.has_key('x') and pt.has_key('y') :
  10. ... poly.append((pt['x'], pt['y']))
  11. ... poly = tuple(poly)
  12. ... p = Polygon(poly)
  13. ... polys.append(p)
  14.  
  15. >>> for pol in polys:
  16. ... print pol.wkt
  17. ...
  18. 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))
  19. 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:

directement avec PostGIS

  1. >>> polys = testpoly._connection.queryAll("""SELECT Srid(the_geom), AsText(the_geom) FROM testpoly""")
  2. >>> polys
  3. [(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))'),[...]')]
  4. >>> polys[0][1]
  5. '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))'
  6. # et donc
  7. b = loads(polys[0][1]]

résultats obtenus avec SQLObject et Shapely:

Shapely - SQLObject - PostGIS

  1. >>> from sqlobject import *
  2. >>> from shapely.wkb import loads
  3. >>> db='postgres://localhost:5432/testpostgis'
  4. >>> connection = connectionForURI(db)
  5. >>> class testpoly(SQLObject):
  6. ... _connection = connection
  7. ... _fromDatabase = True
  8. ... the_geom = StringCol()
  9. ... def _get_geom(self):
  10. ... value = self._SO_get_the_geom()
  11. ... test = loads(value.decode('hex')).wkt
  12. ... return test
  13. ...
  14. >>> a = testpoly.get(1)
  15. >>> a.geom
  16. '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))'
  17. >>>

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:

Shapely - ogr- projections

  1. >>> from shapely.geometry import LineString
  2. >>> from shapely.wkb import loads
  3. >>> from osgeo import ogr
  4. >>> def project(geom, from_epsg, to_epsg):
  5. ... to_srs = ogr.osr.SpatialReference()
  6. ... to_srs.ImportFromEPSG(to_epsg)
  7. ... from_srs = ogr.osr.SpatialReference()
  8. ... from_srs.ImportFromEPSG(from_epsg)
  9. ... ogr_geom = ogr.CreateGeometryFromWkb(geom.wkb)
  10. ... ogr_geom.AssignSpatialReference(from_srs)
  11. ... ogr_geom.TransformTo(to_srs)
  12. ... return loads(ogr_geom.ExportToWkb())
  13. >>> ligne = LineString([[4,50], [5,51]])
  14. >>> ligne.wkt
  15. LINESTRING (4 50, 5 51)
  16. >>> ligneproj = project(ligne, from_epsg=4326, to_epsg=31370)
  17. >>> ligneproj.wkt
  18. '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é)

traitement

  1. >>> import math
  2. >>> from shapely.geometry import LineString
  3. >>> list_x = [ 2.*math.pi/5000*i for i in range(5001) ]
  4. >>> courbe_sinus = [ (list_x[i],math.sin(list_x[i])) for i in range(5001) ]
  5. >>> courbe_cosinus = [ (list_x[i],math.cos(list_x[i])) for i in range(5001) ]
  6. >>> sinus = LineString(courbe_sinus)
  7. >>> cosinus = LineString(courbe_cosinus)
  8. >>> intersections = sinus.intersection(cosinus)
  9. >>> intersections.geom_type
  10. 'MultiPoint'
  11. >>> intersections.wkt
  12. '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:

Shapely - array numpy

  1. >>> line1.xy
  2. (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:

dessine points, lignes

  1. import pylab
  2. from shapely.wkt import loads
  3.  
  4. def plot_coords(ax, ob):
  5. x, y = ob.xy
  6. ax.plot(x, y, 'o', color='r', ms=20)
  7. # 'o', rond, 'r', couleur rouge, 20, taille)
  8. def plot_line(ax, ob):
  9. x, y = ob.xy
  10. ax.plot(x, y, color='b')
  11. def plot_multipt(ax, ob):
  12. for pt in ob:
  13. plot_coords(ax,pt)
  14. def plot_multilignes(ax, ob):
  15. for ligne in ob:
  16. plot_line(ax,ligne)
  17.  
  18. ax = pylab.gca() #ouverture d'une "fenêtre" pylab
  19.  
  20. ligne = loads('LINESTRING (3 1, 4 4, 5 5, 5 6)')
  21. point = loads('POINT (4 5)')
  22. # ensuite suivant les cas
  23. plot_coords(ax,point) #-> point
  24. plot_line(ax, ligne) #-> ligne
  25. plot_coords(ax, ligne) #-> points
  26.  
  27. 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.

dessine aire

  1. from descartes.patch import PolygonPatch
  2.  
  3. def plot_poly(ax,ob):
  4. patch = PolygonPatch(ob, facecolor='#6699cc', edgecolor='#6699cc')
  5. ax.add_patch(patch)
  6.  
  7. def plot_multipol(ax,ob):
  8. for polygon in ob:
  9. plot_poly(ax,polygon)
  10. poly = loads('POLYGON ((1 2, 2 5, 6 5, 7 2, 5 1, 1 2))')
  11. 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)))')
  12.  
  13. plot_poly(ax,poly)
  14. 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 Yumlyuml.me/edit/37c88fdd )

 

Classe

  1. #!/usr/bin/env python
  2. # encoding: utf-8
  3. """
  4. Plot_shapely.py
  5. Martin Laloux 2010
  6. """
  7. import pylab
  8. from shapely.wkt import loads
  9. from descartes.patch import PolygonPatch
  10. class Plot_shapely(object):
  11. """dessin de géométries Shapely avec matplotlib - pylab"""
  12.  
  13. def __init__(self,obj,ax, coul=None, alph=1):
  14. """
  15. Paramètres
  16. ----------------------------------------------------------------
  17. ax de pylab, obj = objet géometrique, coul = couleur matplotlib, alph = transparence
  18.  
  19. Exemple
  20. --------------------------------------------------------------------
  21. >>> from shapely.wkt import loads
  22. >>> ax = pylab.gca()
  23. >>> ligne = loads('LINESTRING (3 1, 4 4, 5 5, 5 6)')
  24. >>> a = Plot_shapely(ligne,ax,'r', 0.5)
  25. >>> a.plot
  26. ou directement
  27. >>> Plot_shapely(ligne,ax,'#FFEC00'').plot
  28. >>> pylab.show()
  29. """
  30. self.obj = obj
  31. self.type = obj.geom_type
  32. self.ax = ax
  33. self.coul= coul
  34. self.alph=alph
  35.  
  36. def plot_coords(self):
  37. """points"""
  38. x, y = self.obj.xy
  39. self.ax.plot(x, y, 'o', color=self.coul)
  40.  
  41. def plot_ligne(self):
  42. """lignes"""
  43. x, y = self.obj.xy
  44. self.ax.plot(x, y, color=self.coul, alpha=self.alph, linewidth=3)
  45.  
  46. def plot_polygone(self):
  47. """polygones"""
  48. patch = PolygonPatch(self.obj, facecolor=self.coul,edgecolor=self.coul, alpha=self.alph)
  49. self.ax.add_patch(patch)
  50.  
  51. def plot_multi(self):
  52. """multipoints, multilignes,multipolygones + GeometryCollection"""
  53. for elem in self.obj:
  54. Plot_shapely(elem, self.ax, self.coul,self.alph).plot
  55.  
  56.  
  57. @property
  58. def plot(self):
  59. """dessine en fonction du type géometrique"""
  60. if self.type == 'Point':
  61. self.plot_coords()
  62. elif self.type == 'Polygon':
  63. self.plot_polygone()
  64. elif self.type == 'LineString':
  65. self.plot_ligne()
  66. elif "Multi" in self.type:
  67. """ex. MultiPolygon"""
  68. self.plot_multi()
  69. elif self.type == 'GeometryCollection':
  70. self.plot_multi()
  71. elif self.type == 'LinearRing':
  72. self.plot_line()
  73. else:
  74. raise ValueError("inconnu au bataillon: %s" % self.type)
  75.  
  76. if __name__ == '__main__':
  77. import doctest
  78. 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.


Site officiel : Shapely
Site officiel : matplotlib
Site officiel : descartes
Site officiel : d9im
Site officiel : GEOS
Autres Liens : manuel Shapely 1.05
Autres Liens : manuel Shapely 1.2

Commentaires

Une petite erreur

Hello,

Merci pour cette présentation, mais il y a une petite erreurs:
LinearRing n’étends pas Polygon mais LineString.

Je trouve également que les illustrations en dessus et en dessous de la phrase « 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. » supère que les polygones sont plein et les mutipolygones sont troué, mais les deux peuvent être troué !

CU
Stéph

Vous n'avez pas bien lu  et

Vous n'avez pas bien lu  et ce n'est pas une petite erreur

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

C'est de la topologie; Dans les shapefiles, les géométries de PostGIS ou d'Oracle Spatial, lignes et polygones (LinearStrig et LinearRing) sont 2 géométries/topologies différentes. Une surface est constituée d'un contour et d'une aire, un contour (LinearRing) délimite toujours une aire, et une aire est toujours limitée par un contour, et ceci n'a rien à voir avec une ligne (LineString).
C'est pourquoi Sean Gillies l'a bien inclus dans sa classe Polygon, même si son implémentation pratique utilise la classe Point et la classe LineString.

Or le diagramme de classe s'adresse au module Shapely.

Poster un nouveau commentaire

Le contenu de ce champ sera maintenu privé et ne sera pas affiché publiquement.