The code snippets on this page need the following imports if you’re outside the pyqgis console:
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)
7. Gestione della Geometria
Punti, linee e poligoni che rappresentano un elemento spaziale sono comunemente indicati come geometrie. In QGIS sono rappresentati con la QgsGeometry
class.
Alcune volte una geometria é effettivamente una collezione di geometrie (parti singole) piú semplici. Se contiene un tipo di geometria semplice, la chiameremo punti multipli, string multi linea o poligoni multipli. Ad esempio, un Paese formato da piú isole puó essere rappresentato come un poligono multiplo.
Le coordinate delle geometrie possono essere in qualsiasi sistema di riferimento delle coordinate (CRS). Quando si estraggono delle caratteristiche da un vettore, le geometrie associate avranno le coordinate nel CRS del vettore.
Description and specifications of all possible geometries construction and relationships are available in the OGC Simple Feature Access Standards for advanced details.
7.1. Costruzione della Geometria
PyQGIS provides several options for creating a geometry:
dalle coordinate
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)
Coordinates are given using
QgsPoint
class orQgsPointXY
class. The difference between these classes is thatQgsPoint
supports M and Z dimensions.A Polyline (Linestring) is represented by a list of points.
A Polygon is represented by a list of linear rings (i.e. closed linestrings). The first ring is the outer ring (boundary), optional subsequent rings are holes in the polygon. Note that unlike some programs, QGIS will close the ring for you so there is no need to duplicate the first point as the last.
Le geometrie a parti multiple vanno ad un livello successivo: punti multipli é una lista di punti, una stringa multi linea é una linea di linee ed un poligono multiplo é una lista di poligoni.
da well-known text (WKT)
geom = QgsGeometry.fromWkt("POINT(3 4)") print(geom)
da well-known binary (WKB)
1g = QgsGeometry() 2wkb = bytes.fromhex("010100000000000000000045400000000000001440") 3g.fromWkb(wkb) 4 5# print WKT representation of the geometry 6print(g.asWkt())
7.2. Accedere alla 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. Predicati ed Operazioni delle Geometrie
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.250080013
3Area (m2): 752000605894.2937
4Area (km2): 752000.6058942936
5Zimbabwe
6Perimeter (m): 2865021.3323912495
7Area (m2): 389250992553.95465
8Area (km2): 389250.99255395465
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))
É possibile trovare molti esempi di algoritmi che sono inclusi in QGIS ed utilizzare questi metodi per analizzare e trasformare i dati vettoriali. Di seguito i link al codice di alcuni di questi.
Distance and area using the
QgsDistanceArea
class: Distance matrix algorithm