跳转至

7 几何处理⚓︎

此页面上的代码片段需要导入以下模块：

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

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

有时，一种几何实际上是简单（单部件）几何的集合。另一种几何形状称为多部件几何。如果它只包含一种类型的简单几何，我们称之为多点、多线或多多边形。例如，由多个岛组成的国家可以表示为多多边形。

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

有关所有可访问的几何结构和关系说明，请参阅OGC简单要素访问标准，以获取更详细的信息。

7.1 几何构造⚓︎

PyQGIS提供了几种创建几何的选项：

  • 坐标
1
2
3
4
5
6
7
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维度。

折线（Linestring）由一系列点表示。

多边形由线环列表（即闭合的线段）表示。第一个环是外环（边界），第二个可选项线环是多边形中的孔。请注意，与某些程序不同，QGIS会为你闭合环，因此无需将第一个点复制为最后一个。

多部件几何图形更进一步：多点是一个点列表，多线是一个线列表，多多边形是一个多边形列表。

  • WKT
1
2
gem = QgsGeometry.fromWkt("POINT(3 4)")
print(geom)
  • WKB
1
2
3
4
5
6
g = QgsGeometry()
wkb = bytes.fromhex("010100000000000000000045400000000000001440")
g.fromWkb(wkb)

#使用WKT打印几何
print(g.asWkt())

7.2 访问几何⚓︎

首先，你应该找出几何类型。wkbType() 方法是其中一种方法。它从QgsWkbTypes.Type 枚举中返回一个值。

1
2
3
4
5
6
7
8
9
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

作为替代方案，可以使用type() ，从QgsWkbTypes.GeometryType 枚举中返回值。

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

1
2
3
4
5
6
print(QgsWkbTypes.displayString(gPnt.wkbType()))
# output: 'Point'
print(QgsWkbTypes.displayString(gLine.wkbType()))
# output: 'LineString'
print(QgsWkbTypes.displayString(gPolygon.wkbType()))
# output: 'Polygon'

还有一个辅助函数 isMultipart()可以确定几何是否是多部件 。

从几何中提取信息，每种矢量类型都有访问器函数。以下是如何使用这些访问器的示例：

1
2
3
4
5
6
gPnt.asPoint()
# output: <QgsPointXY: POINT(1 1)>
gLine.asPolyline()
# output: [<QgsPointXY: POINT(1 1)>, <QgsPointXY: POINT(2 2)>]
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()。

不用管几何的类型，直接遍历所有几何是可以的，例如：

1
2
3
4
5
6
7
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)

1
2
3
4
5
geom = QgsGeometry.fromWkt( 'LineString( 0 0, 10 10 )' )
for part in geom.parts():
  print(part.asWkt())

# LineString (0 0, 10 10)

1
2
3
4
5
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() 方法修改所有几何。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
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 图层中每个国家的面积和周长。

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# 访问'countries'图层
layer = QgsProject.instance().mapLayersByName('countries')[0]

# 过滤以Z开头的国家，然后获取其要素
query = '"name" LIKE \'Z%\''
features = layer.getFeatures(QgsFeatureRequest().setFilterExpression(query))


# 现在遍历要素，执行几何计算并打印结果
for f in features:
  geom = f.geometry()
  name = f.attribute('NAME')
  print(name)
  print('Area: ', geom.area())
  print('Perimeter: ', geom.length())


# Zambia
# Area:  62.822790653431205
# Perimeter:  50.65232014052552
# Zimbabwe
# Area:  33.41113559136521
# Perimeter:  26.608288555013935

现在，你已经计算并打印了几何图形的面积和周长。但是，你可能会很快注意到这些值很奇怪。这是因为当使用 QgsGeometry类中的area()和length() 方法计算时，面积和周长不会考虑CRS。可以使用更强大的QgsDistanceArea 类计算面积和周长，它可以执行基于椭球的计算：

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
d = QgsDistanceArea()
d.setEllipsoid('WGS84')

layer = QgsProject.instance().mapLayersByName('countries')[0]

# 过滤以Z开头的国家，然后获取其要素
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)) ））
  # 打印（“面积（m2）：” ， d 。measureArea （GEOM ））

  # 计算并重新打印面积，单位为平方公里
  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

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

1
2
3
4
5
6
7
8
9
d = QgsDistanceArea()
d.setEllipsoid('WGS84')

# 让我们创造两个点
# 圣诞老人是一个工作狂，他需要放个暑假，让我们看一下他家离特内里费有多远
santa = QgsPointXY(25.847899, 66.543456)
tenerife = QgsPointXY(-16.5735, 28.0443)

print("Distance in meters: ", d.measureLine(santa, tenerife))

你可以在QGIS中找到许多算法示例，并使用这些方法来分析和转换矢量数据。以下是一些代码的链接。