7. Gestione della Geometria

Suggerimento

I frammenti di codice di questa pagina necessitano delle seguenti importazioni se sei è al di fuori della console pyqgis:

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

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 per dettagli avanzati nel documento OGC Simple Feature Access Standards .

7.1. Costruzione della Geometria

PyQGIS offre diverse opzioni per la creazione di una geometria:

• dalle coordinate

  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)
  

  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)

  g = QgsGeometry()
  wkb = bytes.fromhex("010100000000000000000045400000000000001440")
  g.fromWkb(wkb)
  
  # print WKT representation of the geometry
  print(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.

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

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.

print(QgsWkbTypes.displayString(gPnt.wkbType()))
# output: 'Point'
print(QgsWkbTypes.displayString(gLine.wkbType()))
# output: 'LineString'
print(QgsWkbTypes.displayString(gPolygon.wkbType()))
# 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:

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

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)

gc = QgsGeometryCollection()
gc.fromWkt('GeometryCollection( Point(1 2), Point(11 12), LineString(33 34, 44 45))')
print(gc[1].asWkt())

Point (11 12)

È anche possibile modificare ciascuna parte della geometria usando il metodo 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. 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.

# 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

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.

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

In alternativa, potresti voler conoscere la distanza tra due punti.

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

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.

• Distanza e area utilizzando la classe QgsDistanceArea: Distance matrix algorithm.

• Lines to polygons algorithm