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:

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

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

    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)
    

    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)

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

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

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.

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

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:

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

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

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

 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

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.

 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.250294601
3Area (m2): 751989035032.9031
4Area (km2): 751989.0350329031
5Zimbabwe
6Perimeter (m): 2865021.3325076113
7Area (m2): 389267821381.6008
8Area (km2): 389267.8213816008

Alternativamente, es posible que desee conocer la distancia y el rumbo entre dos puntos.

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

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.