7. Gestione della Geometria

Suggerimento

I frammenti di codice di questa pagina necessitano delle seguenti importazioni se sei è al di fuori della console 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)

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.

La descrizione e le specifiche di tutte le possibili geometrie e relazioni sono disponibili nel documento OGC Simple Feature Access Standards per dettagli avanzati.

7.1. Costruzione della Geometria

PyQGIS offre diverse opzioni per la creazione di una geometria:

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

    Le coordinate vengono fornite utilizzando la classe QgsPoint o la classe QgsPointXY. La differenza tra queste classi è che QgsPoint supporta le dimensioni M e Z.

    Una Polilinea (Linestring) è rappresentata da un elenco di punti.

    Un Poligono è rappresentato da un elenco di anelli lineari (cioè di linee chiuse). Il primo anello è l’anello esterno (confine), mentre eventuali anelli successivi sono buchi nel poligono. Da notare che, a differenza di alcuni programmi, QGIS chiude l’anello per voi, quindi non è necessario duplicare il primo punto come ultimo.

    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

Per prima cosa, devi trovare il tipo di geometria. Il metodo wkbType() è quello da utilizzare. Restituisce un valore dell’enumerazione QgsWkbTypes.Type.

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

In alternativa, si può usare il metodo type() che restituisce un valore dell’enumerazione QgsWkbTypes.GeometryType.

Puoi usare la funzione displayString() per ottenere un tipo di geometria leggibile in chiaro.

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

Esiste anche una funzione di aiuto isMultipart() per scoprire se una geometria è multipart o meno.

Per estrarre informazioni da una geometria esistono funzioni di selezione per ogni tipologia di vettore. Ecco un esempio di utilizzo di queste funzioni di selezione:

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

Le tuple (x,y) non sono tuple reali, ma oggetti QgsPoint, i cui valori sono accessibili con i metodi x() e y().

Per le geometrie multiparte esistono funzioni di accesso simili: asMultiPoint(), asMultiPolyline() e asMultiPolygon().

È possibile iterare su tutte le parti di una geometria, indipendentemente dal tipo di geometria. Ad esempio

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)

È anche possibile modificare ciascuna parte della geometria usando il metodo QgsGeometry.parts().

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 utilizza la libreria GEOS per operazioni geometriche avanzate come i predicati geometrici (contains(), intersects(), …) e operazioni di insieme (combine(), difference(), …). Può anche calcolare le proprietà geometriche delle geometrie, come l’area (nel caso dei poligoni) o le lunghezze (per poligoni e linee).

Vediamo un esempio che combina l’iterazione degli elementi in un dato layer e l’esecuzione di alcuni calcoli geometrici basati sulle loro geometrie. Il codice seguente calcolerà e stamperà l’area e il perimetro di ogni countrie nel layer ``countries”” del nostro progetto QGIS.

Il codice seguente presuppone che layer sia un oggetto QgsVectorLayer che ha una tipologia Poligono.

 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

Ora hai calcolato e stampato le aree e i perimetri delle geometrie. Tuttavia, puoi subito notare che i valori sono strani. Questo perché le aree e i perimetri non tengono conto dei SR quando vengono calcolati con i metodi area() e length() della classe QgsGeometry. Per un calcolo più preciso di aree e distanze, occorre utilizzare la classe QgsDistanceArea, che può eseguire calcoli basati su ellissoidi:

Il codice seguente presuppone che layer sia un oggetto QgsVectorLayer che ha una tipologia Poligono.

 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

In aggiunta, potresti voler conoscere la distanza e la direzione tra due punti.

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

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