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 :

from qgis.core import (
  QgsGeometry,
  QgsGeometryCollection,
  QgsPoint,
  QgsPointXY,
  QgsWkbTypes,
  QgsProject,
  QgsFeatureRequest,
  QgsVectorLayer,
  QgsDistanceArea,
  QgsUnitTypes,
  QgsCoordinateTransform,
  QgsCoordinateReferenceSystem
)

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

  gPnt = QgsGeometry.fromPointXY(QgsPointXY(1,1))
  print(gPnt)
  gLine = QgsGeometry.fromPolyline([QgsPoint(1, 1), QgsPoint(2, 2)])
  print(gLine)
  gPolygon = QgsGeometry.fromPolygonXY([[QgsPointXY(1, 1),
      QgsPointXY(2, 2), QgsPointXY(2, 1)]])
  print(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)

  g = QgsGeometry()
  wkb = bytes.fromhex("010100000000000000000045400000000000001440")
  g.fromWkb(wkb)
  
  # print WKT representation of the geometry
  print(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.

if gPnt.wkbType() == QgsWkbTypes.Point:
  print(gPnt.wkbType())
  # output: 1 for Point
if gLine.wkbType() == QgsWkbTypes.LineString:
  print(gLine.wkbType())
  # output: 2 for LineString
if gPolygon.wkbType() == QgsWkbTypes.Polygon:
  print(gPolygon.wkbType())
  # 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.

print(QgsWkbTypes.displayString(gPnt.wkbType()))
# output: 'Point'
print(QgsWkbTypes.displayString(gLine.wkbType()))
# output: 'LineString'
print(QgsWkbTypes.displayString(gPolygon.wkbType()))
# 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 :

print(gPnt.asPoint())
# output: <QgsPointXY: POINT(1 1)>
print(gLine.asPolyline())
# output: [<QgsPointXY: POINT(1 1)>, <QgsPointXY: POINT(2 2)>]
print(gPolygon.asPolygon())
# 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().

geom = QgsGeometry.fromWkt( 'MultiPoint( 0 0, 1 1, 2 2)' )
for part in geom.parts():
  part.transform(QgsCoordinateTransform(
    QgsCoordinateReferenceSystem("EPSG:4326"),
    QgsCoordinateReferenceSystem("EPSG:3111"),
    QgsProject.instance())
  )

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

# let's access the 'countries' layer
layer = QgsProject.instance().mapLayersByName('countries')[0]

# let's filter for countries that begin with Z, then get their features
query = '"name" LIKE \'Z%\''
features = layer.getFeatures(QgsFeatureRequest().setFilterExpression(query))

# now loop through the features, perform geometry computation and print the results
for f in features:
  geom = f.geometry()
  name = f.attribute('NAME')
  print(name)
  print('Area: ', geom.area())
  print('Perimeter: ', geom.length())

Zambia
Area:  62.822790653431205
Perimeter:  50.65232014052552
Zimbabwe
Area:  33.41113559136521
Perimeter:  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.

d = QgsDistanceArea()
d.setEllipsoid('WGS84')

layer = QgsProject.instance().mapLayersByName('countries')[0]

# let's filter for countries that begin with Z, then get their features
query = '"name" LIKE \'Z%\''
features = layer.getFeatures(QgsFeatureRequest().setFilterExpression(query))

for f in features:
  geom = f.geometry()
  name = f.attribute('NAME')
  print(name)
  print("Perimeter (m):", d.measurePerimeter(geom))
  print("Area (m2):", d.measureArea(geom))

  # let's calculate and print the area again, but this time in square kilometers
  print("Area (km2):", d.convertAreaMeasurement(d.measureArea(geom), QgsUnitTypes.AreaSquareKilometers))

Zambia
Perimeter (m): 5539361.250294601
Area (m2): 751989035032.9031
Area (km2): 751989.0350329031
Zimbabwe
Perimeter (m): 2865021.3325076113
Area (m2): 389267821381.6008
Area (km2): 389267.8213816008

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

d = QgsDistanceArea()
d.setEllipsoid('WGS84')

# Let's create two points.
# Santa claus is a workaholic and needs a summer break,
# lets see how far is Tenerife from his home
santa = QgsPointXY(25.847899, 66.543456)
tenerife = QgsPointXY(-16.5735, 28.0443)

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

• Distance et surface à l’aide de la classe QgsDistanceArea: algorithme Matrice des distances

• algorithme Lignes vers polygones