7. Manejo de Geometría

Consejo

Los fragmentos de código en esta página necesitan las siguientes adiciones si está fuera de la consola de pyqgis:

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

Los puntos, las cadenas de líneas y los polígonos que representan una característica espacial se conocen comúnmente como geometrías. En QGIS están representados con la clase QgsGeometry.

A veces una geometría es realmente una colección simple (partes simples) geométricas. Tal geometría se llama geometría de múltiples partes. Si contiene un tipo de geometría simple, lo llamamos un punto múltiple, lineas múltiples o polígonos múltiples. Por ejemplo, un país consiste en múltiples islas que se pueden representar como un polígono múltiple.

Las coordenadas de las geometrías pueden estar en cualquier sistema de referencia de coordenadas (SRC). Cuando extrae características de una capa, las geometrías asociadas tendrán sus coordenadas en el SRC de la capa.

La descripción y las especificaciones de todas las geometrías posibles, la construcción y las especificaciones están disponibles en los Estándares de Acceso a Funciones Simples de OGC para detalles avanzados.

7.1. Construcción de Geometría

PyQGIS proporciona varias opciones para crear una geometría:

• desde coordenadas

  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)
  

  Las coordenadas son dadas usando la clase QgsPoint o la clase QgsPointXY. La diferencia entre estas clases es que soporten QgsPoint dimensione M y Z.

  Una polilínea (cadena lineal) es representada por una lista de puntos.

  Un polígono está representado por una lista de anillos lineales (es decir, cadenas de líneas cerradas). El primer anillo es el anillo exterior (límite), los anillos subsiguientes opcionales son agujeros en el polígono. Tenga en cuenta que, a diferencia de algunos programas, QGIS cerrará el anillo por usted, por lo que no es necesario duplicar el primer punto como el último.

  Las geometrías multi-parte van un nivel más allá: multi-punto es una lista de puntos, multi-linea es una lista de polilíneas y multi-polígono es una lista de polígonos.

• desde well-known text (WKT)

  geom = QgsGeometry.fromWkt("POINT(3 4)")
  print(geom)
  

• desde well-known binary (WKB)

  g = QgsGeometry()
  wkb = bytes.fromhex("010100000000000000000045400000000000001440")
  g.fromWkb(wkb)
  
  # print WKT representation of the geometry
  print(g.asWkt())
  

7.2. Acceso a Geometría

Primero, debe averiguar el tipo de geometría. El método wkbType() es el que se va a utilizar. Devuelve un valor de la enumeración 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

Como alternativa, se puede utilizar el método type() que devuelve un valor de la enumeración QgsWkbTypes.GeometryType.

Puede usar la función displayString() para obtener un tipo de geometría legible por humanos.

print(QgsWkbTypes.displayString(gPnt.wkbType()))
# output: 'Point'
print(QgsWkbTypes.displayString(gLine.wkbType()))
# output: 'LineString'
print(QgsWkbTypes.displayString(gPolygon.wkbType()))
# output: 'Polygon'

Point
LineString
Polygon

También hay una función de ayuda isMultipart() para saber si una geometría es multiparte o no.

Para extraer información de una geometría, existen funciones de acceso para cada tipo de vector. A continuación, se muestra un ejemplo sobre cómo utilizar estos accesos:

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

Las tuplas (x, y) no son tuplas reales, son objetos QgsPoint, los valores son accesibles con los métodos x() y y().

Para geometrías multiparte, existen funciones de acceso similares: asMultiPoint(), asMultiPolyline() y asMultiPolygon().

Es posible iterar sobre todas las partes de una geometría, independientemente del tipo de geometría. P.ej.

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)

También es posible modificar cada parte de la geometría usando el método 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. Geometría predicados y Operaciones

QGIS utiliza la biblioteca GEOS para operaciones de geometría avanzadas como predicados de geometría (contains(), intersects(), …) y hacer operaciones (combine(), difference(), …). También puede calcular propiedades geométricas de geometrías, como área (en el caso de polígonos) o longitudes (para polígonos y líneas).

Veamos un ejemplo que combina iterar sobre las entidades en una capa determinada y realizar algunos cálculos geométricos basados en sus geometrías. El siguiente código calculará e imprimirá el área y el perímetro de cada país en la capa de países dentro de nuestro proyecto tutorial QGIS.

El siguiente código asume que layer es un obleto QgsVectorLayer que tiene el tipo de entidad Polígono.

# 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

Ahora ha calculado e impreso las áreas y perímetros de las geometrías. Sin embargo, puede notar rápidamente que los valores son extraños. Esto se debe a que las áreas y los perímetros no tienen en cuenta el CRS cuando se calculan con los métodos area() y length() desde la clase QgsGeometry. Para un cálculo de área y distancia más potente, la clase QgsDistanceArea se puede utilizar, que puede realizar cálculos basados en elipsoides:

El siguiente código asume que layer es un obleto QgsVectorLayer que tiene el tipo de entidad Polígono.

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

Alternativamente, puede querer saber la distancia entre dos puntos.

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

Puede encontrar muchos ejemplos de algoritmos que se incluyen en QGIS y utilizan estos métodos para analizar y transformar los datos vectoriales. Aquí hay algunos enlaces al código de algunos de ellos.

• Distancia y área usando QgsDistanceArea class: Distance matrix algorithm

• Lines to polygons algorithm