7. Manipulação Geométrica

Dica

Os trechos de código desta página precisam das seguintes importações se você estiver fora do console do pyqgis:

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

Pontos, linhas e polígonos que representam um feição espacial são comumente referidos como geometrias. No QGIS, eles são representados com a classe QgsGeometry.

Às vezes, uma geometria é realmente uma coleção dex simples geometrias (single-part). Tal geometria é chamada de geometria de várias partes. Se ele contém apenas um tipo de simples geometria, podemos chamar de multi-ponto, multi-cadeia linear ou multi-polígono. Por exemplo, um país que consiste de múltiplas ilhas pode ser representado como um sistema multi-polígono.

As coordenadas de geometrias podem estar em qualquer sistema de referência de coordenadas (SRC). Ao buscar feições a partir de uma camada, geometrias associadas terão coordenadas no SRC da camada.

A descrição e as especificações de todas as possíveis construções e relações de geometrias estão disponíveis no OGC Simple Feature Access Standards para detalhes mais avançados.

7.1. Construção de Geométria

O PyQGIS fornece várias opções para criar uma geometria:

  • a partir das coordenadas

    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)
    

    As coordinadas são dadas usando a classe QgsPoint ou a classe QgsPointXY. A diferençae entre estas classes é que a classe QgsPoint suporta M e Z dimensões.

    A Polyline (Linestring) is represented by a list of points.

    Um polígono é representado por uma lista de anéis lineares (ou seja, fformas de linhas fechadas). O primeiro anel é o anel externo (limite); os anéis subsequentes opcionais são orifícios no polígono. Observe que, diferentemente de alguns programas, o QGIS fechará o anel para você, não sendo necessário duplicar o primeiro ponto como o último.

    Geometrias multi-parte passam para um nível maior: multi-ponto é uma lista de pontos, multi-cadeia linear é uma lista de cadeias lineares e multi-polígono é uma lista de polígonos.

  • a partir de textos conhecidos (WKT)

    geom = QgsGeometry.fromWkt("POINT(3 4)")
    print(geom)
    
  • a partir de binários conhecidos (WKB)

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

7.2. Acesso a Geometria

First, you should find out the geometry type. The wkbType() method is the one to use. It returns a value from the QgsWkbTypes.Type enumeration.

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

As an alternative, one can use the type() method which returns a value from the QgsWkbTypes.GeometryType enumeration.

You can use the displayString() function to get a human readable geometry type.

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

There is also a helper function isMultipart() to find out whether a geometry is multipart or not.

To extract information from a geometry there are accessor functions for every vector type. Here’s an example on how to use these accessors:

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

Nota

The tuples (x,y) are not real tuples, they are QgsPoint objects, the values are accessible with x() and y() methods.

For multipart geometries there are similar accessor functions: asMultiPoint(), asMultiPolyline() and asMultiPolygon().

It is possible to iterate over all the parts of a geometry, regardless of the geometry’s type. E.g.

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)

It’s also possible to modify each part of the geometry using QgsGeometry.parts() method.

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 ((-10334726.79314761981368065 -5360105.10101188533008099),(-10462133.82917750626802444 -5217484.34365727473050356),(-10589398.51346865110099316 -5072020.358805269934237))

7.3. Operações e Predicados Geométricos

QGIS uses GEOS library for advanced geometry operations such as geometry predicates (contains(), intersects(), …) and set operations (combine(), difference(), …). It can also compute geometric properties of geometries, such as area (in the case of polygons) or lengths (for polygons and lines).

Let’s see an example that combines iterating over the features in a given layer and performing some geometric computations based on their geometries. The below code will compute and print the area and perimeter of each country in the countries layer within our tutorial QGIS project.

The following code assumes layer is a QgsVectorLayer object that has Polygon feature type.

 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

Now you have calculated and printed the areas and perimeters of the geometries. You may however quickly notice that the values are strange. That is because areas and perimeters don’t take CRS into account when computed using the area() and length() methods from the QgsGeometry class. For a more powerful area and distance calculation, the QgsDistanceArea class can be used, which can perform ellipsoid based calculations:

The following code assumes layer is a QgsVectorLayer object that has Polygon feature type.

 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.250294596
3Area (m2): 751989035032.9031
4Area (km2): 751989.0350329031
5Zimbabwe
6Perimeter (m): 2865021.332507607
7Area (m2): 389267821381.6009
8Area (km2): 389267.82138160086

Alternatively, you may want to know the distance and bearing between two 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))

Você pode encontrar muitos exemplo de algoritmos que estão incluídos no QGIS e usar esses métodos para analisar e transformar dados vetoriais. Aqui estão alguns links para o código de alguns deles.