7. Afhandeling van geometrie

Hint

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

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

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

  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)

  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.

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.

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:

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

Het is mogelijk door alle delen van een geometrie te gaan, ongeacht het type geometrie. Bijv.

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)

Het is ook mogelijk elk deel van de geometrie aan te passen met de methode 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. 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.

# 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.822790653431014
Perimeter:  50.65232014052552
Zimbabwe
Area:  33.41113559136511
Perimeter:  26.608288555013935

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.

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

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

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

• Afstand en gebied gebruiken de klasse QgsDistanceArea: Algoritme Afstandsmatrix

• Algoritme Lijnen naar polygonen