7. 几何体处理

提示

如果您在pyqgis控制台之外，则此页面上的代码片段需要以下导入：

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

表示空间要素的点、线串和多边形通常称为几何图形。在QGIS中，它们用 QgsGeometry 班级。

有时，一个几何图形实际上是简单(单个部分)几何图形的集合。这样的几何图形称为多零件几何图形。如果它只包含一种类型的简单几何图形，我们称之为多点、多线串或多多边形。例如，由多个岛屿组成的国家可以表示为多个多边形。

几何图形的坐标可以在任何坐标参考系(CRS)中。从层中提取要素时，关联的几何图形将在该层的CRS中具有坐标。

所有可能的几何结构和关系的描述和规范可在 OGC Simple Feature Access Standards 了解更高级的详细信息。

7.1. 几何构造

PyQGIS提供了几个用于创建几何体的选项：

• 从坐标开始

  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)
  

  坐标是使用 QgsPoint 类或 QgsPointXY 班级。这些类之间的区别在于 QgsPoint 支持M和Z尺寸。

  多段线(直线串)由点列表表示。

  一个多边形由一列线性环(即闭合的线串)表示。第一个环是外环(边界)，可选的后续环是多边形中的孔。请注意，与某些程序不同的是，QGIS将为您关闭环，因此没有必要将第一个点复制为最后一个点。

  多零件几何图形更进一步：多点是点的列表，多线串是线串的列表，多多边形是多边形的列表。

• 从熟知文本(WKT)

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

• 来自熟知二进制(WKB)

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

7.2. 访问几何体

首先，您应该找出几何图形类型。这个 wkbType() 方法是要使用的方法。它返回一个来自 QgsWkbTypes.Type 枚举。

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

作为另一种选择，可以使用 type() 方法，该方法从 QgsWkbTypes.GeometryType 枚举。

print(gLine.type())
# output: 'GeometryType.Line'

您可以使用 displayString() 函数以获取人类可读的几何类型。

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

还有一个帮助器功能 isMultipart() 以确定几何体是否为多部分。

为了从几何图形中提取信息，每个向量类型都有访问器函数。下面是一个关于如何使用这些访问器的示例：

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

备注

元组(x，y)不是真正的元组，它们是 QgsPoint 对象，这些值可以通过 x() 和 y() 方法：研究方法。

对于多部分几何图形，有类似的存取器函数： asMultiPoint() ， asMultiPolyline() 和 asMultiPolygon() 。

可以遍历几何体的所有部分，而不考虑几何体的类型。例如。

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)

还可以使用修改几何的每个部分 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. 几何谓词和运算

QGIS使用GEOS库进行高级几何操作，如几何谓词 (contains() ， intersects() ，…)和集合运算 (combine() ， difference() ，…)。它还可以计算几何图形的几何特性，例如面积(对于多边形)或长度(对于多边形和直线)。

让我们来看一个例子，它结合了迭代给定层中的要素并根据它们的几何形状执行一些几何计算。下面的代码将计算和打印每个国家/地区在 countries 在我们的教程QGIS项目中的层。

以下代码假定 layer 是一种 QgsVectorLayer 具有多边形要素类型的对象。

# 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.82279065343119
Perimeter:  50.65232014052552
Zimbabwe
Area:  33.41113559136517
Perimeter:  26.608288555013935

现在，您已经计算并打印了几何图形的面积和周长。然而，您可能很快就会注意到这些值很奇怪。这是因为在使用 area() 和 length() 方法，方法来自 QgsGeometry 班级。对于更强大的面积和距离计算， QgsDistanceArea 类，它可以执行基于椭球体的计算：

以下代码假定 layer 是一种 QgsVectorLayer 具有多边形要素类型的对象。

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

或者，您可能想知道两点之间的距离。

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

您可以找到QGIS中包含的许多算法示例，并使用这些方法来分析和转换矢量数据。这里有一些指向其中几个代码的链接。

• 距离和面积使用 QgsDistanceArea 班级： Distance matrix algorithm

• Lines to polygons algorithm