De codesnippers op deze pagina hebben de volgende import nodig als u buiten de console van PyQGIS bent:

 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. Afhandeling van geometrie

Naar punten, lijnen en polygonen die een ruimtelijk object weergeven wordt gewoonlijk verwezen als geometrieën. In QGIS worden zij weergegeven door de klasse QgsGeometry.

Soms is één geometrie in feite een verzameling van enkele (ééndelige) geometrieën. Een dergelijke geometrie wordt een geometrie met meerdere delen genoemd. Als het slechts één type eenvoudige geometrie bevat, noemen we het multi-punt, multi-lijn of multi-polygoon. Een land dat bijvoorbeeld bestaat uit meerdere eilanden kan worden weergegeven als een multi-polygoon.

De coördinaten van geometrieën kunnen in elk coördinaten referentiesysteem (CRS) staan. Bij het ophalen van objecten vanaf een laag, zullen de geassocieerde geometrieën in coördinaten in het CRS van de laag staan.

Beschrijving en specificaties van alle mogelijke constructies van geometrieën en relaties zijn beschikbaar in de OGC Simple Feature Access Standards voor uitgebreide details.

7.1. Construeren van geometrie

PyQGIS verschaft verscheidene opties voor het maken van een geometrie:

  • uit coördinaten

    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)

    Coördinaten worden opgegeven met behulp van de klassen QgsPoint of QgsPointXY. Het verschil tussen deze klassen is dat QgsPoint dimensies M en Z ondersteunt.

    Een Polylijn (Lijn) wordt weergegeven door een lijst met punten.

    Een polygoon wordt weergegeven als een lijst van lineaire ringen (d.i. gesloten lijnen). De eerste ring is de buitenste ring (grens), optionele volgende ringen zijn gaten in de polygoon. Onthoud dat, anders dan andere programma’s, QGIS de ring voor u zal sluiten dus is er geen reden om het eerste punt als laatste te dupliceren.

    Geometrieën die bestaan uit meerdere delen gaan een niveau verder: multi-punt is een lijst van punten, multi-lijnen zijn een lijst van lijnen en multi-polygoon is een lijst van polygonen.

  • uit bekende tekst (WKT)

    geom = QgsGeometry.fromWkt("POINT(3 4)")
print(geom)

  • uit bekende binaire (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. Toegang tot geometrie

Als eerste zou u het type geometrie moeten zoeken, de methode wkbType() is die om te gebruiken. Het geeft een waarde uit de enumeratie QgsWkbTypes.Type terug.

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

Als alternatief kan men de methode type() gebruiken die een waarde teruggeeft uit de enumeratie van de klasse QgsWkbTypes.GeometryType.

U kunt de functie displayString() gebruiken om een voor mensen leesbaar type geometrie te verkrijgen.

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

Er is ook een hulpfunctie isMultipart() om uit te zoeken of een geometrie meerdelig is of niet.

Voor elk type vector zijn er functies voor toegang om informatie uit de geometrie op te halen. Hier is een voorbeeld hoe deze functies te gebruiken:

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)>]]

Notitie

De tuples (x,y) zijn geen echte tuples, zij zijn objecten QgsPoint, de waarden zijn toegankelijk met de methoden x() en y().

Voor meerdelige geometrieën zijn er soortgelijke functies voor toegang: asMultiPoint(), asMultiPolyline() en asMultiPolygon().

7.3. Predicaten en bewerking voor geometrieën

QGIS gebruikt de bibliotheek GEOS voor geavanceerde bewerkingen met geometrieën, zoals de predicaten voor geometrieën (contains(), intersects(), …) en het instellen van bewerkingen (combine(), difference(), …). Het kan ook geometrische eigenschappen van geometrieën berekenen, zoals gebied (in het geval van polygonen) of lengten (voor polygonen en lijnen).

Laten we een voorbeeld bekijken dat het doorlopen van de objecten op een laag combineert met het uitvoeren van enkele geometrische berekeningen, gebaseerd op hun geometrieën. De onderstaande code zal het gebied en de perimeter van elk land op de laag countries in ons project in de handleiding voor QGIS berekenen en afdrukken.

De volgende code gaat ervan uit dat layer een object QgsVectorLayer is.

 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

Nu hebt u de gebieden en perimeters van de geometrieën berekend en afgedrukt. Het zal u echter snel opvallen dat de waarden vreemd zijn. Dat komt omdat gebieden en perimeters geen rekening houden met het CRS bij het berekenen met behulp van de methoden area() en length() uit de klasse QgsGeometry. Voor een meer krachtiger berekening van gebied en afstand kan de klasse QgsDistanceArea worden gebruikt, die op ellipsoïde gebaseerde berekeningen kan uitvoeren:

De volgende code gaat ervan uit dat layer een object QgsVectorLayer is.

 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

Als alternatief zou u misschien de afstand en richting tussen twee punten willen weten.

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

U kunt zoeken naar vele voorbeelden van algoritmen die zijn opgenomen in QGIS en die methoden gebruiken om vectorgegevens te analyseren en te transformeren. Hier zijn enkele koppelingen naar de code van sommige ervan.