19. 网络分析库

提示

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

from qgis.core import (
  QgsVectorLayer,
  QgsPointXY,
)

网络分析库可用于:

  • 从地理数据(多段线矢量图层)创建数学图形

  • 从图论实现基本方法(目前只有Dijkstra算法)

网络分析库是通过从RoadGraph核心插件导出基本函数创建的,现在您可以在插件中使用它的方法,也可以直接从Python控制台使用它的方法。

19.1. 一般信息

简而言之,一个典型的用例可以描述为:

  1. 从地理数据(通常为多段线向量图层)创建图形

  2. 运行图分析

  3. 使用分析结果(例如,将其可视化)

19.2. 构建图表

你需要做的第一件事-是准备输入数据,也就是把一个矢量层转换成一个图形。所有进一步的操作都将使用此图形,而不是层。

作为源,我们可以使用任何折线矢量层。多段线的节点成为图形顶点,多段线的线段成为图形边。如果几个节点具有相同的坐标,则它们是相同的图形顶点。因此,具有公共节点的两条线路相互连接。

此外,在图形创建过程中,可以将任意数量的附加点“固定”(“连接”)到输入向量层。对于每个附加点,将在最近的图顶点或最接近的图边中找到匹配项。在后一种情况下,将分割边并添加新顶点。

边的矢量层属性和长度可用作边的属性。

从矢量层到图形的转换是使用 Builder 编程模式。图是使用所谓的导向器构造的。目前只有一位导演: QgsVectorLayerDirector 。导向器设置将用于从线向量层构建图形的基本设置,构建器使用该线向量层来创建图形。目前,与导演的情况一样,只有一个建筑商: QgsGraphBuilder ,这创造了 QgsGraph 物体。您可能希望实现自己的构建器,以构建与以下库兼容的图形 BGLNetworkX

要计算边缘属性,请使用编程模式 strategy 使用的是。只是暂时的 QgsNetworkDistanceStrategy 策略(考虑到路线的长度)和 QgsNetworkSpeedStrategy (这也考虑了速度)是可用的。您可以实现自己的策略,该策略将使用所有必要的参数。例如,RoadGraph插件使用一种策略,该策略使用属性的边长和速度值来计算行程时间。

现在是时候投入到这个过程中去了。

首先,要使用这个库,我们应该导入分析模块

from qgis.analysis import *

然后是一些创建导演的示例

 1# don't use information about road direction from layer attributes,
 2# all roads are treated as two-way
 3director = QgsVectorLayerDirector(vectorLayer, -1, '', '', '', QgsVectorLayerDirector.DirectionBoth)
 4
 5# use field with index 5 as source of information about road direction.
 6# one-way roads with direct direction have attribute value "yes",
 7# one-way roads with reverse direction have the value "1", and accordingly
 8# bidirectional roads have "no". By default roads are treated as two-way.
 9# This scheme can be used with OpenStreetMap data
10director = QgsVectorLayerDirector(vectorLayer, 5, 'yes', '1', 'no', QgsVectorLayerDirector.DirectionBoth)

要构建导向器,我们应该传递一个矢量层,该矢量层将用作图形结构和每个路段上允许移动的信息(单向或双向移动、正向或反向移动)的源。呼叫如下所示

1director = QgsVectorLayerDirector(vectorLayer,
2                                  directionFieldId,
3                                  directDirectionValue,
4                                  reverseDirectionValue,
5                                  bothDirectionValue,
6                                  defaultDirection)

以下是这些参数含义的完整列表:

  • vectorLayer -用于构建图形的矢量层

  • directionFieldId -属性表字段的索引,存储道路方向信息。如果 -1 ,那么根本就不要使用这些信息。一个整数。

  • directDirectionValue -具有直接方向(从第一个线点移动到最后一个线点)的道路的字段值。一根绳子。

  • reverseDirectionValue -方向相反(从最后一个线点移动到第一个线点)的道路的字段值。一根绳子。

  • bothDirectionValue -双向道路的字段值(对于这种道路,我们可以从第一个点移动到最后一个,从最后一个移动到第一个)。一根绳子。

  • defaultDirection -默认道路方向。该值将用于那些道路位置字段 directionFieldId 未设置或具有与上面指定的三个值中的任何一个不同的值。可能的值包括:

因此,有必要创建用于计算边属性的策略

1# The index of the field that contains information about the edge speed
2attributeId = 1
3# Default speed value
4defaultValue = 50
5# Conversion from speed to metric units ('1' means no conversion)
6toMetricFactor = 1
7strategy = QgsNetworkSpeedStrategy(attributeId, defaultValue, toMetricFactor)

把这个策略告诉主管

director = QgsVectorLayerDirector(vectorLayer, -1, '', '', '', 3)
director.addStrategy(strategy)

现在我们可以使用构建器,它将创建图形。这个 QgsGraphBuilder 类构造函数有几个参数:

  • crs -要使用的坐标系。必选参数。

  • otfEnabled -是否使用“On the Fly”重新投影。默认情况下 True (使用OTF)。

  • topologyTolerance -拓扑容差。默认值为0。

  • ellipsoidID -使用椭球体。默认情况下为“WGS84”。

# only CRS is set, all other values are defaults
builder = QgsGraphBuilder(vectorLayer.crs())

我们还可以定义几个点,这些点将用于分析。例如

startPoint = QgsPointXY(1179720.1871, 5419067.3507)
endPoint = QgsPointXY(1180616.0205, 5419745.7839)

现在一切都已就绪,我们可以构建图表并将这些点“绑定”到它上

tiedPoints = director.makeGraph(builder, [startPoint, endPoint])

构建图表可能需要一些时间(这取决于层中要素的数量和层的大小)。 tiedPoints 是一个带有坐标的列表,这些点被“绑在一起”。构建操作完成后,我们可以获得图形并使用它进行分析

graph = builder.graph()

使用下一段代码,我们可以获得点的顶点索引

startId = graph.findVertex(tiedPoints[0])
endId = graph.findVertex(tiedPoints[1])

19.3. 图形分析

网络分析被用来寻找两个问题的答案:哪些顶点是相连的,以及如何找到最短路径。为了解决这些问题,网络分析库提供了Dijkstra算法。

Dijkstra的算法寻找从图的一个顶点到所有其他顶点的最短路径以及优化参数值。结果可以表示为最短路径树。

最短路径树是具有以下属性的有向加权图(或者更准确地说是树):

  • 只有一个折点没有传入边,即树的根

  • 所有其他折点只有一条传入边

  • 如果折点B可以从折点A到达,则从A到B的路径是唯一可用路径,并且它在该图上是最优(最短)的

要获得最短路径树,请使用以下方法 shortestTree()dijkstra()QgsGraphAnalyzer 班级。建议您使用 dijkstra() 方法,因为它的运行速度更快,使用内存也更高效。

这个 shortestTree() 当您要遍历最短路径树时,方法非常有用。它始终创建一个新的图形对象(QgsGraph)并接受三个变量:

  • source -输入图

  • startVertexIdx -树上点的索引(树的根)

  • criterionNum -要使用的边缘属性数(从0开始)。

tree = QgsGraphAnalyzer.shortestTree(graph, startId, 0)

这个 dijkstra() 方法具有相同的参数,但返回两个数组。在第一个数组元素 n 包含传入边的索引,如果没有传入边,则为-1。在第二数组元素中 n 包含从树根到顶点的距离 n 如果是顶点,则为Double_Max n 是从根本上无法到达的。

(tree, cost) = QgsGraphAnalyzer.dijkstra(graph, startId, 0)

下面是一些非常简单的代码,用于使用使用 shortestTree() 方法(在中选择线串图层 Layers 面板,并用您自己的坐标替换)。

警告

仅以此代码为例,它创建了许多 QgsRubberBand 对象,在大型数据集上可能会很慢。

 1from qgis.core import *
 2from qgis.gui import *
 3from qgis.analysis import *
 4from qgis.PyQt.QtCore import *
 5from qgis.PyQt.QtGui import *
 6
 7vectorLayer = QgsVectorLayer('testdata/network.gpkg|layername=network_lines', 'lines')
 8director = QgsVectorLayerDirector(vectorLayer, -1, '', '', '', QgsVectorLayerDirector.DirectionBoth)
 9strategy = QgsNetworkDistanceStrategy()
10director.addStrategy(strategy)
11builder = QgsGraphBuilder(vectorLayer.crs())
12
13pStart = QgsPointXY(1179661.925139,5419188.074362)
14tiedPoint = director.makeGraph(builder, [pStart])
15pStart = tiedPoint[0]
16
17graph = builder.graph()
18
19idStart = graph.findVertex(pStart)
20
21tree = QgsGraphAnalyzer.shortestTree(graph, idStart, 0)
22
23i = 0
24while (i < tree.edgeCount()):
25  rb = QgsRubberBand(iface.mapCanvas())
26  rb.setColor (Qt.red)
27  rb.addPoint (tree.vertex(tree.edge(i).fromVertex()).point())
28  rb.addPoint (tree.vertex(tree.edge(i).toVertex()).point())
29  i = i + 1

相同的事情,但使用 dijkstra() 方法

 1from qgis.core import *
 2from qgis.gui import *
 3from qgis.analysis import *
 4from qgis.PyQt.QtCore import *
 5from qgis.PyQt.QtGui import *
 6
 7vectorLayer = QgsVectorLayer('testdata/network.gpkg|layername=network_lines', 'lines')
 8
 9director = QgsVectorLayerDirector(vectorLayer, -1, '', '', '', QgsVectorLayerDirector.DirectionBoth)
10strategy = QgsNetworkDistanceStrategy()
11director.addStrategy(strategy)
12builder = QgsGraphBuilder(vectorLayer.crs())
13
14pStart = QgsPointXY(1179661.925139,5419188.074362)
15tiedPoint = director.makeGraph(builder, [pStart])
16pStart = tiedPoint[0]
17
18graph = builder.graph()
19
20idStart = graph.findVertex(pStart)
21
22(tree, costs) = QgsGraphAnalyzer.dijkstra(graph, idStart, 0)
23
24for edgeId in tree:
25  if edgeId == -1:
26    continue
27  rb = QgsRubberBand(iface.mapCanvas())
28  rb.setColor (Qt.red)
29  rb.addPoint (graph.vertex(graph.edge(edgeId).fromVertex()).point())
30  rb.addPoint (graph.vertex(graph.edge(edgeId).toVertex()).point())

19.3.1. 查找最短路径

为了找到两点之间的最佳路径,使用了以下方法。在构建图形时,这两个点(起点A和终点B)都“绑定”到图形。然后使用 shortestTree()dijkstra() 方法建立以起点A为根的最短路径树,并在同一棵树中找到终点B,开始遍历从B点到A点的最短路径树。整个算法可以写成:

1assign T = B
2while T != B
3    add point T to path
4    get incoming edge for point T
5    look for point TT, that is start point of this edge
6    assign T = TT
7add point A to path

在这一点上,我们有了路径,其形式为倒置的顶点列表(顶点从终点到起点以相反的顺序列出),在通过该路径行进期间将被访问。

以下是QGIS Python控制台的示例代码(您可能需要在TOC中加载并选择一个线串图层,并将代码中的坐标替换为您的),它使用 shortestTree() 方法

 1from qgis.core import *
 2from qgis.gui import *
 3from qgis.analysis import *
 4
 5from qgis.PyQt.QtCore import *
 6from qgis.PyQt.QtGui import *
 7
 8vectorLayer = QgsVectorLayer('testdata/network.gpkg|layername=network_lines', 'lines')
 9builder = QgsGraphBuilder(vectorLayer.sourceCrs())
10director = QgsVectorLayerDirector(vectorLayer, -1, '', '', '', QgsVectorLayerDirector.DirectionBoth)
11
12startPoint = QgsPointXY(1179661.925139,5419188.074362)
13endPoint = QgsPointXY(1180942.970617,5420040.097560)
14
15tiedPoints = director.makeGraph(builder, [startPoint, endPoint])
16tStart, tStop = tiedPoints
17
18graph = builder.graph()
19idxStart = graph.findVertex(tStart)
20
21tree = QgsGraphAnalyzer.shortestTree(graph, idxStart, 0)
22
23idxStart = tree.findVertex(tStart)
24idxEnd = tree.findVertex(tStop)
25
26if idxEnd == -1:
27    raise Exception('No route!')
28
29# Add last point
30route = [tree.vertex(idxEnd).point()]
31
32# Iterate the graph
33while idxEnd != idxStart:
34    edgeIds = tree.vertex(idxEnd).incomingEdges()
35    if len(edgeIds) == 0:
36        break
37    edge = tree.edge(edgeIds[0])
38    route.insert(0, tree.vertex(edge.fromVertex()).point())
39    idxEnd = edge.fromVertex()
40
41# Display
42rb = QgsRubberBand(iface.mapCanvas())
43rb.setColor(Qt.green)
44
45# This may require coordinate transformation if project's CRS
46# is different than layer's CRS
47for p in route:
48    rb.addPoint(p)

下面是相同的示例,但使用了 dijkstra() 方法

 1from qgis.core import *
 2from qgis.gui import *
 3from qgis.analysis import *
 4
 5from qgis.PyQt.QtCore import *
 6from qgis.PyQt.QtGui import *
 7
 8vectorLayer = QgsVectorLayer('testdata/network.gpkg|layername=network_lines', 'lines')
 9director = QgsVectorLayerDirector(vectorLayer, -1, '', '', '', QgsVectorLayerDirector.DirectionBoth)
10strategy = QgsNetworkDistanceStrategy()
11director.addStrategy(strategy)
12
13builder = QgsGraphBuilder(vectorLayer.sourceCrs())
14
15startPoint = QgsPointXY(1179661.925139,5419188.074362)
16endPoint = QgsPointXY(1180942.970617,5420040.097560)
17
18tiedPoints = director.makeGraph(builder, [startPoint, endPoint])
19tStart, tStop = tiedPoints
20
21graph = builder.graph()
22idxStart = graph.findVertex(tStart)
23idxEnd = graph.findVertex(tStop)
24
25(tree, costs) = QgsGraphAnalyzer.dijkstra(graph, idxStart, 0)
26
27if tree[idxEnd] == -1:
28    raise Exception('No route!')
29
30# Total cost
31cost = costs[idxEnd]
32
33# Add last point
34route = [graph.vertex(idxEnd).point()]
35
36# Iterate the graph
37while idxEnd != idxStart:
38    idxEnd = graph.edge(tree[idxEnd]).fromVertex()
39    route.insert(0, graph.vertex(idxEnd).point())
40
41# Display
42rb = QgsRubberBand(iface.mapCanvas())
43rb.setColor(Qt.red)
44
45# This may require coordinate transformation if project's CRS
46# is different than layer's CRS
47for p in route:
48    rb.addPoint(p)

19.3.2. 可供使用的领域

顶点A的可用区域是可从顶点A访问的图顶点的子集,并且从A到这些顶点的路径的成本不大于某个值。

下面的例子可以更清楚地说明这一点:“有一个消防站。消防车可以在5分钟内到达城市的哪些地方?10分钟?15分钟?”这些问题的答案是消防站的可用区域。

为了找到可用的区域,我们可以使用 dijkstra() 的方法。 QgsGraphAnalyzer 班级。将成本数组的元素与预定义的值进行比较就足够了。如果成本 [i] 小于或等于预定义的值,则顶点i在可用区域内,否则在可用区域外。

一个更困难的问题是确定可用区域的边界。下边界是仍可访问的顶点集,上边界是不可访问的顶点集。事实上,这很简单:它是基于最短路径树的边的可用边界,边的源顶点是可访问的,而边的目标顶点是不可访问的。

下面是一个例子

 1director = QgsVectorLayerDirector(vectorLayer, -1, '', '', '', QgsVectorLayerDirector.DirectionBoth)
 2strategy = QgsNetworkDistanceStrategy()
 3director.addStrategy(strategy)
 4builder = QgsGraphBuilder(vectorLayer.crs())
 5
 6
 7pStart = QgsPointXY(1179661.925139, 5419188.074362)
 8delta = iface.mapCanvas().getCoordinateTransform().mapUnitsPerPixel() * 1
 9
10rb = QgsRubberBand(iface.mapCanvas())
11rb.setColor(Qt.green)
12rb.addPoint(QgsPointXY(pStart.x() - delta, pStart.y() - delta))
13rb.addPoint(QgsPointXY(pStart.x() + delta, pStart.y() - delta))
14rb.addPoint(QgsPointXY(pStart.x() + delta, pStart.y() + delta))
15rb.addPoint(QgsPointXY(pStart.x() - delta, pStart.y() + delta))
16
17tiedPoints = director.makeGraph(builder, [pStart])
18graph = builder.graph()
19tStart = tiedPoints[0]
20
21idStart = graph.findVertex(tStart)
22
23(tree, cost) = QgsGraphAnalyzer.dijkstra(graph, idStart, 0)
24
25upperBound = []
26r = 1500.0
27i = 0
28tree.reverse()
29
30while i < len(cost):
31    if cost[i] > r and tree[i] != -1:
32        outVertexId = graph.edge(tree [i]).toVertex()
33        if cost[outVertexId] < r:
34            upperBound.append(i)
35    i = i + 1
36
37for i in upperBound:
38    centerPoint = graph.vertex(i).point()
39    rb = QgsRubberBand(iface.mapCanvas())
40    rb.setColor(Qt.red)
41    rb.addPoint(QgsPointXY(centerPoint.x() - delta, centerPoint.y() - delta))
42    rb.addPoint(QgsPointXY(centerPoint.x() + delta, centerPoint.y() - delta))
43    rb.addPoint(QgsPointXY(centerPoint.x() + delta, centerPoint.y() + delta))
44    rb.addPoint(QgsPointXY(centerPoint.x() - delta, centerPoint.y() + delta))