7. Manipulation de la géométrie

Indication

Les extraits de code sur cette page nécessitent les importations suivantes si vous êtes en dehors de la console pyqgis :

 1from qgis.core import (
 2  QgsGeometry,
 3  QgsGeometryCollection,
 4  QgsPoint,
 5  QgsPointXY,
 6  QgsWkbTypes,
 7  QgsProject,
 8  QgsFeatureRequest,
 9  QgsVectorLayer,
10  QgsDistanceArea,
11  QgsUnitTypes,
12  QgsCoordinateTransform,
13  QgsCoordinateReferenceSystem
14)

Les points, les linestring et les polygons qui représentent une entité spatiale sont communément appelés géométries. Dans QGIS, ils sont représentés par la classe QgsGeometry

Parfois, une entité correspond à une collection d’éléments géométriques simples (d’un seul tenant). Une telle géométrie est appelée multi-parties. Si elle ne contient qu’un seul type de géométrie, il s’agit de multi-points, de multi-lignes ou de multi-polygones. Par exemple, un pays constitué de plusieurs îles peut être représenté par un multi-polygone.

Les coordonnées des géométries peuvent être dans n’importe quel système de coordonnées de référence (SCR). Lorsqu’on accède aux entités d’une couche, les géométries correspondantes auront leurs coordonnées dans le SCR de la couche.

La description et les spécifications liées à la construction et aux relations entre géométries sont disponibles dans OGC Simple Feature Access Standards pour plus de détails.

7.1. Construction de géométrie

PyQGIS fournit plusieurs options pour créer une géométrie:

  • à partir des coordonnées

    1gPnt = QgsGeometry.fromPointXY(QgsPointXY(1,1))
    2print(gPnt)
    3gLine = QgsGeometry.fromPolyline([QgsPoint(1, 1), QgsPoint(2, 2)])
    4print(gLine)
    5gPolygon = QgsGeometry.fromPolygonXY([[QgsPointXY(1, 1),
    6    QgsPointXY(2, 2), QgsPointXY(2, 1)]])
    7print(gPolygon)
    

    Les coordonnées sont données en utilisant la classe QgsPoint ou la classe QgsPointXY. La différence entre ces classes est que la classe QgsPoint supporte les dimensions M et Z.

    Une polyligne (Linestring) est représentée par une liste de points.

    Un Polygone est représenté par une liste d’anneaux linéaires (c’est-à-dire des chaînes de caractères fermées). Le premier anneau est l’anneau extérieur (limite), les anneaux suivants facultatifs sont des trous dans le polygone. Notez que contrairement à certains programmes, QGIS fermera l’anneau pour vous, il n’est donc pas nécessaire de dupliquer le premier point comme le dernier.

    Les géométries multi-parties sont d’un niveau plus complexe: les multipoints sont une succession de points, les multilignes une succession de lignes et les multipolygones une succession de polygones.

  • depuis un Well-Known-Text (WKT)

    geom = QgsGeometry.fromWkt("POINT(3 4)")
    print(geom)
    
  • depuis un Well-Known-Binary (WKB)

    1g = QgsGeometry()
    2wkb = bytes.fromhex("010100000000000000000045400000000000001440")
    3g.fromWkb(wkb)
    4
    5# print WKT representation of the geometry
    6print(g.asWkt())
    

7.2. Accéder à la Géométrie

Tout d’abord, vous devez trouver le type de géométrie. La méthode wkbType() est celle à utiliser. Elle renvoie une valeur provenant de l’énumération QgsWkbTypes.Type.

1if gPnt.wkbType() == QgsWkbTypes.Point:
2  print(gPnt.wkbType())
3  # output: 1 for Point
4if gLine.wkbType() == QgsWkbTypes.LineString:
5  print(gLine.wkbType())
6  # output: 2 for LineString
7if gPolygon.wkbType() == QgsWkbTypes.Polygon:
8  print(gPolygon.wkbType())
9  # output: 3 for Polygon

Comme alternative, on peut utiliser la méthode type() qui retourne une valeur de l’énumération QgsWkbTypes.GeometryType.

Vous pouvez utiliser la fonction displayString() pour obtenir un type de géométrie lisible par l’homme.

1print(QgsWkbTypes.displayString(gPnt.wkbType()))
2# output: 'Point'
3print(QgsWkbTypes.displayString(gLine.wkbType()))
4# output: 'LineString'
5print(QgsWkbTypes.displayString(gPolygon.wkbType()))
6# output: 'Polygon'
Point
LineString
Polygon

Il existe également une fonction d’aide isMultipart() pour savoir si une géométrie est multipartite ou non.

Pour extraire des informations d’une géométrie, il existe des fonctions d’accès pour chaque type de vecteur. Voici un exemple d’utilisation de ces accesseurs :

1print(gPnt.asPoint())
2# output: <QgsPointXY: POINT(1 1)>
3print(gLine.asPolyline())
4# output: [<QgsPointXY: POINT(1 1)>, <QgsPointXY: POINT(2 2)>]
5print(gPolygon.asPolygon())
6# output: [[<QgsPointXY: POINT(1 1)>, <QgsPointXY: POINT(2 2)>, <QgsPointXY: POINT(2 1)>, <QgsPointXY: POINT(1 1)>]]

Note

Les tuples (x,y) ne sont pas de vrais tuples, ce sont des objets QgsPoint, les valeurs sont accessibles avec les méthodes x() et y().

Pour les géométries en plusieurs parties, il existe des fonctions d’accès similaires : asMultiPoint(), asMultiPolyline() et asMultiPolygon().

Il est possible d’itérer sur les parties d’une géométrie, peu importe son type. Par exemple:

geom = QgsGeometry.fromWkt( 'MultiPoint( 0 0, 1 1, 2 2)' )
for part in geom.parts():
  print(part.asWkt())
Point (0 0)
Point (1 1)
Point (2 2)
geom = QgsGeometry.fromWkt( 'LineString( 0 0, 10 10 )' )
for part in geom.parts():
  print(part.asWkt())
LineString (0 0, 10 10)
gc = QgsGeometryCollection()
gc.fromWkt('GeometryCollection( Point(1 2), Point(11 12), LineString(33 34, 44 45))')
print(gc[1].asWkt())
Point (11 12)

Il est aussi possible de modifier chaque partie de la géométrie à l’aide de la méthode QgsGeometry.parts().

1geom = QgsGeometry.fromWkt( 'MultiPoint( 0 0, 1 1, 2 2)' )
2for part in geom.parts():
3  part.transform(QgsCoordinateTransform(
4    QgsCoordinateReferenceSystem("EPSG:4326"),
5    QgsCoordinateReferenceSystem("EPSG:3111"),
6    QgsProject.instance())
7  )
8
9print(geom.asWkt())
MultiPoint ((-10334728.12541878595948219 -5360106.25905461423099041),(-10462135.16126426123082638 -5217485.4735023295506835),(-10589399.84444035589694977 -5072021.45942386891692877))

7.3. Prédicats et opérations géométriques

QGIS utilise la bibliothèque GEOS pour des opérations de géométrie avancées telles que les prédicats de géométrie (contains(), intersects(), …) et les opérations de prédicats (combine(), difference(), …). Il peut également calculer les propriétés géométriques des géométries, telles que l’aire (dans le cas des polygones) ou les longueurs (pour les polygones et les lignes).

Voyons un exemple qui combine l’itération sur les entités d’une couche donnée et l’exécution de certains calculs géométriques basés sur leurs géométries. Le code ci-dessous permet de calculer et d’imprimer la superficie et le périmètre de chaque pays dans la couche pays dans le cadre de notre projet tutoriel QGIS.

Le code suivant suppose que layer est un objet QgsVectorLayer qui a le type d’entité Polygon.

 1# let's access the 'countries' layer
 2layer = QgsProject.instance().mapLayersByName('countries')[0]
 3
 4# let's filter for countries that begin with Z, then get their features
 5query = '"name" LIKE \'Z%\''
 6features = layer.getFeatures(QgsFeatureRequest().setFilterExpression(query))
 7
 8# now loop through the features, perform geometry computation and print the results
 9for f in features:
10  geom = f.geometry()
11  name = f.attribute('NAME')
12  print(name)
13  print('Area: ', geom.area())
14  print('Perimeter: ', geom.length())
1Zambia
2Area:  62.822790653431205
3Perimeter:  50.65232014052552
4Zimbabwe
5Area:  33.41113559136521
6Perimeter:  26.608288555013935

Vous avez maintenant calculé et imprimé les surfaces et les périmètres des géométries. Vous pouvez cependant rapidement remarquer que les valeurs sont étranges. C’est parce que les surfaces et les périmètres ne prennent pas en compte le CRS lorsqu’ils sont calculés à l’aide des méthodes area() et length() de la classe QgsGeometry. Pour un calcul de surface et de distance plus puissant, la classe QgsDistanceArea peut être utilisée, ce qui permet d’effectuer des calculs basés sur des ellipsoïdes :

Le code suivant suppose que layer est un objet QgsVectorLayer qui a le type d’entité Polygon.

 1d = QgsDistanceArea()
 2d.setEllipsoid('WGS84')
 3
 4layer = QgsProject.instance().mapLayersByName('countries')[0]
 5
 6# let's filter for countries that begin with Z, then get their features
 7query = '"name" LIKE \'Z%\''
 8features = layer.getFeatures(QgsFeatureRequest().setFilterExpression(query))
 9
10for f in features:
11  geom = f.geometry()
12  name = f.attribute('NAME')
13  print(name)
14  print("Perimeter (m):", d.measurePerimeter(geom))
15  print("Area (m2):", d.measureArea(geom))
16
17  # let's calculate and print the area again, but this time in square kilometers
18  print("Area (km2):", d.convertAreaMeasurement(d.measureArea(geom), QgsUnitTypes.AreaSquareKilometers))
1Zambia
2Perimeter (m): 5539361.250294601
3Area (m2): 751989035032.9031
4Area (km2): 751989.0350329031
5Zimbabwe
6Perimeter (m): 2865021.3325076113
7Area (m2): 389267821381.6008
8Area (km2): 389267.8213816008

Vous pouvez également vouloir connaître la distance et le bearing entre deux points.

 1d = QgsDistanceArea()
 2d.setEllipsoid('WGS84')
 3
 4# Let's create two points.
 5# Santa claus is a workaholic and needs a summer break,
 6# lets see how far is Tenerife from his home
 7santa = QgsPointXY(25.847899, 66.543456)
 8tenerife = QgsPointXY(-16.5735, 28.0443)
 9
10print("Distance in meters: ", d.measureLine(santa, tenerife))

Vous trouverez de nombreux exemples d’algorithmes inclus dans QGIS et utiliser ces méthodes pour analyser et modifier les données vectorielles. Voici des liens vers le code de quelques-uns.