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

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

7. Manipulation de la géométrie

Les points, les linestring et les polygons qui représentent une entite 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 de toutes les constructions et relations géométriques possibles sont disponibles dans le document « OGC Simple Feature Access Standards » <https://www.opengeospatial.org/standards/sfa>`_ pour des détails avancés.

7.1. Construction de géométrie

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

  • à partir des coordonnées

    1
2
3
4
5
6
7
    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)

    1
2
3
4
5
6
    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.

1
2
3
4
5
6
7
8
9
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.

1
2
3
4
5
6
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 :

1
2
3
4
5
6
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().

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 entites 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’entite Polygon.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
# 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 \'Zu%\''
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())

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
Zubin Potok
Area:  0.040717371293465573
Perimeter:  0.9406133328077781
Zulia
Area:  3.708060762610232
Perimeter:  17.172123598311487
Zuid-Holland
Area:  0.4204687950359031
Perimeter:  4.098878517120812
Zug
Area:  0.027573510374275363
Perimeter:  0.7756605461489624

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’entite Polygon.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
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 \'Zu%\''
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))

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
Zubin Potok
Perimeter (m): 87581.40256396442
Area (m2): 369302069.18814206
Area (km2): 369.30206918814207
Zulia
Perimeter (m): 1891227.0945423362
Area (m2): 44973645460.19726
Area (km2): 44973.64546019726
Zuid-Holland
Perimeter (m): 331941.8000214341
Area (m2): 3217213408.4100943
Area (km2): 3217.213408410094
Zug
Perimeter (m): 67440.22483063207
Area (m2): 232457391.52097562
Area (km2): 232.45739152097562

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
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.