19. ネットワーク分析ライブラリ

ヒント

pyqgisコンソールを使わない場合、このページにあるコードスニペットは次のインポートが必要です:

from qgis.core import (
  QgsVectorLayer,
  QgsPointXY,
)

ネットワーク分析ライブラリは次のことに使われます:

  • 地理データ（ポリライン・ベクタレイヤ）から数学的グラフを作る

  • グラフ理論の基本的な手法を実装する（現在はダイクストラ法のみ）

ネットワーク分析ライブラリはRoadGraphコアプラグインから基本機能をエクスポートすることによって作成されました。今はそのメソッドをプラグインで、またはPythonのコンソールから直接使用できます。

19.1. 一般情報

手短に言えば、一般的なユースケースは、次のように記述できます。

  1. 地理データから地理的なグラフ（たいていはポリラインベクタレイヤ）を作成する

  2. グラフ分析を実行する

  3. 分析結果を利用する（例えば、その可視化）

19.2. グラフを構築する

最初にする必要がある事は --- 入力データを準備すること、つまりベクタレイヤをグラフに変換することです。これからのすべての操作は、レイヤではなく、このグラフを使用します。

ソースとしてどんなポリラインベクタレイヤも使用できます。ポリラインの頂点は、グラフの頂点となり、ポリラインのセグメントは、グラフの辺です。いくつかのノードが同じ座標を持っている場合、それらは同じグラフの頂点です。だから共通のノードを持つ2つの線は互いに接続しています。

さらに、グラフの作成時には、入力ベクタレイヤに好きな数だけ追加の点を「固定する」（「結びつける」）ことが可能です。追加の点それぞれについて、対応箇所---最も近いグラフの頂点または最も近いグラフの辺、が探し出されます。後者の場合、辺は分割されて新しい頂点が追加されるでしょう。

ベクタレイヤ属性と辺の長さは、辺のプロパティとして使用できます。

ベクタレイヤからグラフへの変換は、Builder プログラミングパターンを使って行われます。グラフはいわゆるDirectorを使って構築されます。今のところDirectorは1つしかありません: QgsVectorLayerDirector。ディレクタは、ビルダがラインベクタレイヤからグラフを作るのに使う基本的な設定を行います。現在のところ、ディレクタと同様に、ビルダは1つしかありません: QgsGraphBuilder で、 QgsGraph オブジェクトを作成します。BGL <https://www.boost.org/doc/libs/1_48_0/libs/graph/doc/index.html>`_ や NetworkX などのライブラリと互換性のあるグラフを構築する独自のビルダを実装することもできます。

辺のプロパティを計算するために、プログラミングパターン strategy が使用されます。今のところ（経路の長さを考慮する） QgsNetworkDistanceStrategy 戦略と（速度も考慮する）:class:QgsNetworkSpeedStrategy <qgis.analysis.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 が設定されていない、または上記の3つの値のいずれとも異なる値を持つ道路に使用されます。設定可能な値は以下の通りです：

それから、辺のプロパティを計算するための戦略を作成することが必要です

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 --- 「オンザフライ」再投影を使用するかどうか。デフォルトは 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. グラフ分析

ネットワーク分析はこの二つの質問に対する答えを見つけるために使用されます：どの頂点が接続されているか、どのように最短経路を検索するか。これらの問題を解決するため、ネットワーク解析ライブラリではダイクストラのアルゴリズムを提供しています。

ダイクストラ法は、グラフの１つの頂点から他のすべての頂点への最短ルートおよび最適化パラメーターの値を見つけます。結果は、最短経路木として表現できます。

最短経路木は、以下の性質を持つ有向加重グラフ（より正確には木）です：

  • 流入する辺がない頂点が1つだけあります - 木の根

  • 他のすべての頂点には流入する辺が１つだけあります

  • 頂点Bが頂点Aから到達可能である場合、このグラフ上のAからBへの経路は、単一の利用可能な経路であり、それは最適（最短）です

最短経路木を得るには QgsGraphAnalyzer クラスの shortestTree() と dijkstra() メソッドを使用します。dijkstra() メソッドを使用することを推奨します。なぜなら、このメソッドはより速く動作し、より効率的にメモリを使用するからです。

shortestTree() メソッドは最短経路木を様々な角度から考えたいときに便利です。それは常に新しいグラフオブジェクト(QgsGraph)を生成し、3つの変数を受け取ります：

  • source --- 入力グラフ

  • startVertexIdx --- 木上の点のインデックス（木の根）

  • criterionNum --- 使用する辺のプロパティの数（0から始まる）。

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

dijkstra() メソッドは同じ引数を持ちますが、2つの配列を返します。最初の配列要素 n には入力辺のインデックス、または、入力辺がないときは-1が格納されます。2 番目の配列の要素 n には、木の根から頂点 n までの距離、または頂点 n が根から到達できない場合は DOUBLE_MAX が格納されます。

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

以下は、 shortestTree() メソッドで作成したグラフを使用して、最短経路木を表示する非常に簡単なコードです（ レイヤ パネルでラインストリングレイヤを選択し、座標をあなたのものに置き換えてください）。

警告

このコードは例としてだけ使ってください。これは 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. 最短経路を見つける

2点間の最適経路を求めるには、次のようなアプローチが用いられます。両方の点（始点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コンソールで shortestTree() メソッドを使うサンプルコードです（TOCにあるラインストリングレイヤを読み込んで選択し、コード内の座標を自分のものに置き換える必要があるかもしれません）

 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分で到達できるのは市内のどの部分ですか？」。これらの質問への回答が消防署の利用可能領域です。

利用可能領域を見つけるには、 QgsGraphAnalyzer クラスの dijkstra() メソッドを使用します。cost配列の要素をあらかじめ定義された値と比較するだけです。もしcost[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))