本节代码片段需导入以下模块:
| from qgis.core import (
QgsVectorLayer,
QgsPointXY,
)
|
19 网络分析库
网络分析库可用于:
- 从地理数据(矢量线折线图层)创建图
- 实现图论中的基本算法(目前只有Dijkstra的算法)
网络分析库是通过从RoadGraph
核心插件导出基本函数创建的,现在你可以在插件中使用它的方法,也可以直接从Python控制台使用它。
19.1 一般信息
简而言之,一个典型用例可以描述为:
- 从地理数据创建图(矢量折线图层)
- 运行图算法
- 使用分析结果
19.2 构建一个图
你需要做的第一件事是准备输入数据,也就是将矢量图层转换为图。所有进一步的操作都将使用这个图,而不是图层。
作为数据源,我们可以使用任何折线矢量图层。折线的节点成为图的顶点,折线的线段成为图的边。如果几个节点具有相同的坐标,那么它们就是相同的图顶点。因此,具有公共节点的两条线就会相互连接。
此外,在图创建过程中,可以将任意数量的附加点“固定”(“绑”)到输入矢量图层。对于每个附加点,将找到一个匹配——最近的顶点或最近的边。在后一种情况下,边将被分割并添加一个新顶点。
矢量图层属性和边的长度可以用作边的属性。
使用Builder编程模式完成从矢量图层到图的转换。图是使用所谓的控制器构造的。目前只有一个控制器:QgsVectorLayerDirector
。控制器设置了基本的设置——这些设置将用于从线矢量图层构造图,构建器用来创建图。目前,控制器一样,只有一个构建器存在:QgsGraphBuilder
,它可以创建QgsGraph
对象。你可能希望实现自己的构建器,以建立一个与BGL或NetworkX等库兼容的图。
为了计算边属性,使用编程模式策略。目前只有QGSNetworkDistanceTreatgy
策略(考虑到路线的长度)和QgsNetworkSpeedStrategy
(也考虑到速度)可用。你可以实现自己的策略,使用所有必要的参数。例如,RoadGraph插件使用的策略是使用边长度和属性中的速度值来计算行程时间。
是时候深入研究这个过程了。
首先,要使用这个库,我们应该导入分析模块
| from qgis.analysis import *
|
然后是一些创建控制器的示例:
| # 不要使用图层属性中有关道路方向的信息,所有道路都是双向的
director = QgsVectorLayerDirector(vectorLayer, -1, '', '', '', QgsVectorLayerDirector.DirectionBoth)
# 使用第5个字段作为道路的方向信息.
# 正向单方向道路使用 "yes",
# 反向单方向道路使用 "1"
# 因此,双向道路使用 “no”。默认情况道路是双向路。
# 此方案可用于OpenStreetMap数据
director = QgsVectorLayerDirector(vectorLayer, 5, 'yes', '1', 'no', QgsVectorLayerDirector.DirectionBoth)
|
为了构造一个控制器,我们应该传递一个矢量图层,该矢量图层作为图结构的源,以及关于每个路段上允许移动的信息(单向或双向移动、直接或反向),像这样:
| director = QgsVectorLayerDirector(vectorLayer,
directionFieldId,
directDirectionValue,
reverseDirectionValue,
bothDirectionValue,
defaultDirection)
|
以下是这些参数的全部含义:
然后有必要创建用于计算边属性的策略:
| # 包含速度信息的字段索引值
attributeId = 1
# 默认速度
defaultValue = 50
# 转化速度到米制单位 ('1' 表示不转化)
toMetricFactor = 1
strategy = QgsNetworkSpeedStrategy(attributeId, defaultValue, toMetricFactor)
|
告诉控制器这个策略
| director = QgsVectorLayerDirector(vectorLayer, -1, '', '', '', 3)
director.addStrategy(strategy)
|
现在我们可以使用构造器来创建图。QgsGraphBuilder
类构造函数接受几个参数:
crs
——坐标参考系统。必需。
otfEnabled
——使用“on the fly” 重投影. 默认为True
(use OTF).
topologyTolerance
——拓扑容差。默认为0。
ellipsoidID
——参考椭球。默认为“WGS84”。
| # 只设置CRS,,其它值默认
builder = QgsGraphBuilder(vectorLayer.crs())
|
我们还可以定义几个点,这些点将用于分析。例如:
| startPoint = QgsPointXY(1179720.1871, 5419067.3507)
endPoint = QgsPointXY(1180616.0205, 5419745.7839)
|
现在一切就绪,我们可以构建图表并将这些点“连接”到它
| tiedPoints = director.makeGraph(builder, [startPoint, endPoint])
|
构建图可能需要一些时间(这取决于图层中的要素数量和图层大小)。tiedPoints
是一个包含“绑定”点坐标的列表。当构建操作完成时,我们可以得到图并将其用于分析
通过下面的代码,我们可以得到点的顶点索引
| startId = graph.findVertex(tiedPoints[0])
endId = graph.findVertex(tiedPoints[1])
|
19.3 图分析
网络分析用于找到两个问题的答案:哪些顶点是相连的,以及如何找到最短路径。为了解决这些问题,网络分析库提供了Dijkstra算法。
Dijkstra算法找到从图的一个顶点到所有其他顶点的最短路径以及优化参数的值。结果可以表示为最短路径树。
最短路径树是一个有向加权图(或更准确地说是一棵树),具有以下特性:
- 只有一个顶点没有入射边——树的根
- 所有其他顶点只有一条入射边
- 如果顶点B可以从顶点A到达,那么从A到B的路径是唯一可用的路径,并且它在该图上是最优的(最短的)路径
要获得最短路径树,可以使用QgsGraphAnalyzer
类的shortestTree()
和dijkstra()
方法。建议使用dijkstra()
方法,因为它更快,而且更有效地使用内存。
shortestTree()
方法在你想在最短路径树上行走时很有用。它总是创建一个新的图对象(QgsGraph)并接受三个变量:
source
——输入的图
startVertexIdx
——树上的点的索引(树的根)
criterionNum
——使用的边属性的数量(从0开始)
| tree = QgsGraphAnalyzer.shortestTree(graph, startId, 0)
|
dijkstra()
方法有相同的参数,但返回两个数组。在第一个数组中,n元素包含传入边的索引,如果没有传入边则为-1。在第二个数组中,n元素包含从树的根到顶点n的距离,如果顶点n从根部无法到达,则为DOUBLE_MAX。
| (tree, cost) = QgsGraphAnalyzer.dijkstra(graph, startId, 0)
|
下面是一些非常简单的代码,使用shortestTree()
方法创建的图显示最短路径树(在Layers
面板中选择linestring
图层,用你自己的坐标替换)。
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
29 | from qgis.core import *
from qgis.gui import *
from qgis.analysis import *
from qgis.PyQt.QtCore import *
from qgis.PyQt.QtGui import *
vectorLayer = QgsVectorLayer('testdata/network.gpkg|layername=network_lines', 'lines')
director = QgsVectorLayerDirector(vectorLayer, -1, '', '', '', QgsVectorLayerDirector.DirectionBoth)
strategy = QgsNetworkDistanceStrategy()
director.addStrategy(strategy)
builder = QgsGraphBuilder(vectorLayer.crs())
pStart = QgsPointXY(1179661.925139,5419188.074362)
tiedPoint = director.makeGraph(builder, [pStart])
pStart = tiedPoint[0]
graph = builder.graph()
idStart = graph.findVertex(pStart)
tree = QgsGraphAnalyzer.shortestTree(graph, idStart, 0)
i = 0
while (i < tree.edgeCount()):
rb = QgsRubberBand(iface.mapCanvas())
rb.setColor (Qt.red)
rb.addPoint (tree.vertex(tree.edge(i).fromVertex()).point())
rb.addPoint (tree.vertex(tree.edge(i).toVertex()).point())
i = i + 1
|
使用dijkstra()
方法
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
29
30 | from qgis.core import *
from qgis.gui import *
from qgis.analysis import *
from qgis.PyQt.QtCore import *
from qgis.PyQt.QtGui import *
vectorLayer = QgsVectorLayer('testdata/network.gpkg|layername=network_lines', 'lines')
director = QgsVectorLayerDirector(vectorLayer, -1, '', '', '', QgsVectorLayerDirector.DirectionBoth)
strategy = QgsNetworkDistanceStrategy()
director.addStrategy(strategy)
builder = QgsGraphBuilder(vectorLayer.crs())
pStart = QgsPointXY(1179661.925139,5419188.074362)
tiedPoint = director.makeGraph(builder, [pStart])
pStart = tiedPoint[0]
graph = builder.graph()
idStart = graph.findVertex(pStart)
(tree, costs) = QgsGraphAnalyzer.dijkstra(graph, idStart, 0)
for edgeId in tree:
if edgeId == -1:
continue
rb = QgsRubberBand(iface.mapCanvas())
rb.setColor (Qt.red)
rb.addPoint (graph.vertex(graph.edge(edgeId).fromVertex()).point())
rb.addPoint (graph.vertex(graph.edge(edgeId).toVertex()).point())
|
19.3.1 查找最短路径
为了找到两点之间的最佳路径,采用了以下方法。两点(起点A和终点B)在图建立时都被 "捆绑 "在一起。然后使用shortestTree()
或dijkstra()
方法,我们建立以起点A为根的最短路径树。在同一棵树上,我们也找到了终点B,并开始从B点走到A点,整个算法可以写成这样:
| assign T = B
while T != B
add point T to path
get incoming edge for point T
look for point TT, that is start point of this edge
assign T = TT
add point A to path
|
在这一点上,我们有一个路径,其形式是顶点的倒置列表(顶点是按照从终点到起点的相反顺序排列的),这些顶点将在这条路径的行程中被访问。
以下是使用shortestTree()
方法的QGIS Python控制台示例代码(你可能需要在图层目录树中中加载并选择一个线图层,并将代码中的坐标替换为你的坐标)。
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
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49 | from qgis.core import *
from qgis.gui import *
from qgis.analysis import *
from qgis.PyQt.QtCore import *
from qgis.PyQt.QtGui import *
vectorLayer = QgsVectorLayer('testdata/network.gpkg|layername=network_lines', 'lines')
builder = QgsGraphBuilder(vectorLayer.sourceCrs())
director = QgsVectorLayerDirector(vectorLayer, -1, '', '', '', QgsVectorLayerDirector.DirectionBoth)
strategy = QgsNetworkDistanceStrategy()
director.addStrategy(strategy)
startPoint = QgsPointXY(1179661.925139,5419188.074362)
endPoint = QgsPointXY(1180942.970617,5420040.097560)
tiedPoints = director.makeGraph(builder, [startPoint, endPoint])
tStart, tStop = tiedPoints
graph = builder.graph()
idxStart = graph.findVertex(tStart)
tree = QgsGraphAnalyzer.shortestTree(graph, idxStart, 0)
idxStart = tree.findVertex(tStart)
idxEnd = tree.findVertex(tStop)
if idxEnd == -1:
raise Exception('No route!')
# 添加最后一个点
route = [tree.vertex(idxEnd).point()]
# 遍历图
while idxEnd != idxStart:
edgeIds = tree.vertex(idxEnd).incomingEdges()
if len(edgeIds) == 0:
break
edge = tree.edge(edgeIds[0])
route.insert(0, tree.vertex(edge.fromVertex()).point())
idxEnd = edge.fromVertex()
# 显示
rb = QgsRubberBand(iface.mapCanvas())
rb.setColor(Qt.green)
# 如果项目的坐标系和图层的坐标系不一样,则需要坐标转换
for p in route:
rb.addPoint(p)
|
使用dijkstra()
方法
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
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47 | from qgis.core import *
from qgis.gui import *
from qgis.analysis import *
from qgis.PyQt.QtCore import *
from qgis.PyQt.QtGui import *
vectorLayer = QgsVectorLayer('testdata/network.gpkg|layername=network_lines', 'lines')
director = QgsVectorLayerDirector(vectorLayer, -1, '', '', '', QgsVectorLayerDirector.DirectionBoth)
strategy = QgsNetworkDistanceStrategy()
director.addStrategy(strategy)
builder = QgsGraphBuilder(vectorLayer.sourceCrs())
startPoint = QgsPointXY(1179661.925139,5419188.074362)
endPoint = QgsPointXY(1180942.970617,5420040.097560)
tiedPoints = director.makeGraph(builder, [startPoint, endPoint])
tStart, tStop = tiedPoints
graph = builder.graph()
idxStart = graph.findVertex(tStart)
idxEnd = graph.findVertex(tStop)
(tree, costs) = QgsGraphAnalyzer.dijkstra(graph, idxStart, 0)
if tree[idxEnd] == -1:
raise Exception('No route!')
# Total cost
cost = costs[idxEnd]
# 添加最后一个点
route = [graph.vertex(idxEnd).point()]
# 遍历图
while idxEnd != idxStart:
idxEnd = graph.edge(tree[idxEnd]).fromVertex()
route.insert(0, graph.vertex(idxEnd).point())
# 显示
rb = QgsRubberBand(iface.mapCanvas())
rb.setColor(Qt.red)
# 如果项目的坐标系和图层的坐标系不一样,则需要坐标转换
for p in route:
rb.addPoint(p)
|
19.3.2 可达区域
顶点A的可达区域是指可以从顶点A进入的图形顶点的子集,并且从A到这些顶点的路径成本不超过某个值。
这一点可以通过以下例子更清楚地表明。"有一个消防站,消防车可以在5分钟内到达城市的哪些地方?10分钟?15分钟?"。这些问题的答案就是消防站的可达区域。
为了找到可达的区域,我们可以使用QgsGraphAnalyzer
类的dijkstra()
方法。只需将cost数组的元素与预定义的值进行比较。如果cost[i]小于或等于一个预定义的值,那么顶点i就在可达区域内,否则它就在外面。
一个更困难的问题是要得到可达区域的边界。底部边界是仍可访问的顶点集合,而顶部边界是不可访问的顶点集合。事实上这很简单:它是基于最短路径树的边的可达性边界,对于这些边的源顶点是可访问的,而边的目标顶点则不是。
下面是一个例子:
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
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44 | director = QgsVectorLayerDirector(vectorLayer, -1, '', '', '', QgsVectorLayerDirector.DirectionBoth)
strategy = QgsNetworkDistanceStrategy()
director.addStrategy(strategy)
builder = QgsGraphBuilder(vectorLayer.crs())
pStart = QgsPointXY(1179661.925139, 5419188.074362)
delta = iface.mapCanvas().getCoordinateTransform().mapUnitsPerPixel() * 1
rb = QgsRubberBand(iface.mapCanvas())
rb.setColor(Qt.green)
rb.addPoint(QgsPointXY(pStart.x() - delta, pStart.y() - delta))
rb.addPoint(QgsPointXY(pStart.x() + delta, pStart.y() - delta))
rb.addPoint(QgsPointXY(pStart.x() + delta, pStart.y() + delta))
rb.addPoint(QgsPointXY(pStart.x() - delta, pStart.y() + delta))
tiedPoints = director.makeGraph(builder, [pStart])
graph = builder.graph()
tStart = tiedPoints[0]
idStart = graph.findVertex(tStart)
(tree, cost) = QgsGraphAnalyzer.dijkstra(graph, idStart, 0)
upperBound = []
r = 1500.0
i = 0
tree.reverse()
while i < len(cost):
if cost[i] > r and tree[i] != -1:
outVertexId = graph.edge(tree [i]).toVertex()
if cost[outVertexId] < r:
upperBound.append(i)
i = i + 1
for i in upperBound:
centerPoint = graph.vertex(i).point()
rb = QgsRubberBand(iface.mapCanvas())
rb.setColor(Qt.red)
rb.addPoint(QgsPointXY(centerPoint.x() - delta, centerPoint.y() - delta))
rb.addPoint(QgsPointXY(centerPoint.x() + delta, centerPoint.y() - delta))
rb.addPoint(QgsPointXY(centerPoint.x() + delta, centerPoint.y() + delta))
rb.addPoint(QgsPointXY(centerPoint.x() - delta, centerPoint.y() + delta))
|