PyQGisCookbook--几何体处理(八)


几何体处理

如果您在pyqgis控制台之外,则此以下代码段需要以下导入:

from qgis.core import (
  QgsGeometry,
  QgsPoint,
  QgsPointXY,
  QgsWkbTypes,
  QgsProject,
  QgsFeatureRequest,
  QgsDistanceArea
)

表示空间要素的点,线和多边形通常称为几何体。在QGIS中,它们以QgsGeometry类表示 。

有时,一个几何体实际上是一组简单(单个部分)几何体的集合。这种几何形状称为多部分几何体。如果仅包含一种类型的简单几何体,则将其称为多点,多线或多多边形。例如,由多个岛屿组成的国家可以表示为多面。

几何的坐标可以在任何坐标参考系统(CRS)中。从图层获取Feature时,关联的几何体的坐标在图层的CRS中。

OGC简单要素访问标准的高级描述信息中提供了所有可能的几何体构造和关系的描述和规范。

几何体构造

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 class或QgsPointXY class 给出的。这些类之间的区别是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())
    

访问几何体

首先,您应该找出几何体类型。wkbType() 方法是一种可使用的方法。它返回一个QgsWkbTypes.Type 枚举的值。

if gPnt.wkbType() == QgsWkbTypes.Point:
  print(gPnt.wkbType())
  # output: 1 for Point
if gLine.wkbType() == QgsWkbTypes.LineString:
  print(gLine.wkbType())
if gPolygon.wkbType() == QgsWkbTypes.Polygon:
  print(gPolygon.wkbType())
  # output: 3 for Polygon

作为一种替代方法,可以使用type() 方法获得QgsWkbTypes.GeometryType 枚举类型的返回值。

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

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

几何体处理和运算

QGIS使用GEOS库先进的几何操作,如几何处理(contains()intersects(),...)和设置操作(combine()difference(),...)。它还可以计算几何体形状的几何属性,例如面积(对于多边形)或长度(对于多边形和直线)。

让我们看一个示例,该示例将遍历给定图层中的Feature并根据其几何形状执行合并在一起且执行一些几何计算。以下代码将计算和打印在countries层每个国家的面积和周长。

以下代码假定layerQgsVectorLayer具有Polygon要素类型的对象。

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

现在,您已经计算并打印了几何图形的面积和周长。但是,您可能很快注意到这些值很奇怪。这是因为使用QgsGeometry类中的area()length() 方法计算时,面积和周长未考虑CRS 。为了进行更精确的面积和距离计算,可以使用QgsDistanceArea 类,该类可以执行基于椭球的计算:

以下代码假定layerQgsVectorLayer类型的对象,具有Polygon Feature。

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

或者,您可能想知道两点之间的距离和方位。

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中包含的许多算法示例,并使用这些方法来分析和转换矢量数据。这里是其中一些代码的一些链接。

下一个   前一个

©版权所有2002-现在,QGIS项目 最近更新于2020年4月3日09:14。

使用Sphinx使用Read the Docs提供的主题构建。

  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值