重要

翻訳は あなたが参加できる コミュニティの取り組みです。このページは現在 75.00% 翻訳されています。

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

ヒント

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

from qgis.core import (
  QgsVectorLayer,
  QgsPointXY,
)

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

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

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

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

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

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

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

19.1. グラフを構築する

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

ソースとしてどんなポリラインベクタレイヤも使用できます。ポリラインの頂点は、グラフの頂点となり、ポリラインのセグメントは、グラフの辺です。いくつかのノードが同じ座標を持っている場合、それらは同じグラフの頂点です。だから共通のノードを持つ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 などのライブラリと互換性のあるグラフを構築する独自のビルダを実装することもできます。

To calculate edge properties the programming pattern strategy is used. For now only QgsNetworkDistanceStrategy strategy (that takes into account the length of the route) and QgsNetworkSpeedStrategy (that also considers the speed) are available. You can implement your own strategy that will use all necessary parameters.

では手順に行きましょう。

  1. First of all, to use this library we should import the analysis module:

    from qgis.analysis import *

  2. Then some examples for creating a director:

    1# Don't use information about road direction from layer attributes,
2# all roads are treated as two-way
3director = QgsVectorLayerDirector(
4    vectorLayer, -1, "", "", "", QgsVectorLayerDirector.DirectionBoth
5)

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

    To construct a director, we should pass a vector layer that will be used as the source for the graph structure and information about allowed movement on each road segment (one-way or bidirectional movement, direct or reverse direction). The call looks like this (find more details on the parameters at qgis.analysis.QgsVectorLayerDirector):

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

  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)

  4. そして、ディレクタにこの戦略について教えます

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

  5. Now we can use the builder, which will create the graph, using the QgsGraphBuilder class constructor.

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

  6. Also we can define several points, which will be used in the analysis. For example:

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

  7. Now all is in place so we can build the graph and "tie" these points to it:

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

    Building the graph can take some time (which depends on the number of features in a layer and layer size). tiedPoints is a list with coordinates of "tied" points.

  8. When the build operation is finished we can get the graph and use it for the analysis:

    graph = builder.graph()

  9. With the next code we can get the vertex indexes of our points:

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

19.2. グラフ分析

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

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

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

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

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

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

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

The shortestTree() method is useful when you want to walk around the shortest path tree. It always creates a new graph object (QgsGraph) and accepts three variables:

  • source --- 入力グラフ

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

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

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

qgis.analysis.QgsGraphAnalyzer.dijkstra() メソッドは同じ引数を受け取りますが、配列のタプルを返します：

  • 最初の配列では、要素 n は入力エッジのインデックスを含み、入力エッジが存在しない場合は -1 を含む。

  • 2番目の配列は、要素 n に木の根から頂点 n までの距離、または頂点 n が根から到達不能である場合はDOUBLE_MAXを格納します。

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

Here is some very simple code to display the shortest path tree using the graph created with the shortestTree() or the dijkstra() method (select linestring layer in Layers panel and replace coordinates with your own).

警告

このコードは例としてだけ使ってください。これは 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(
 8    "testdata/network.gpkg|layername=network_lines", "lines"
 9)
10director = QgsVectorLayerDirector(
11    vectorLayer, -1, "", "", "", QgsVectorLayerDirector.DirectionBoth
12)
13strategy = QgsNetworkDistanceStrategy()
14director.addStrategy(strategy)
15builder = QgsGraphBuilder(vectorLayer.crs())
16
17pStart = QgsPointXY(1179661.925139, 5419188.074362)
18tiedPoint = director.makeGraph(builder, [pStart])
19pStart = tiedPoint[0]
20
21graph = builder.graph()
22
23idStart = graph.findVertex(pStart)
24
25tree = QgsGraphAnalyzer.shortestTree(graph, idStart, 0)
26
27i = 0
28while i < tree.edgeCount():
29    rb = QgsRubberBand(iface.mapCanvas())
30    rb.setColor(Qt.red)
31    rb.addPoint(tree.vertex(tree.edge(i).fromVertex()).point())
32    rb.addPoint(tree.vertex(tree.edge(i).toVertex()).point())
33    i = i + 1

19.2.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

この時点において、この経路で走行中に訪問される頂点の反転リストの形（頂点は逆順で終点から始点へと列挙されている）で、経路が得られます。

Here is the sample code for QGIS Python Console (you may need to load and select a linestring layer in TOC and replace coordinates in the code with yours) that uses the shortestTree() or dijkstra() method:

 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(
 9    "testdata/network.gpkg|layername=network_lines", "lines"
10)
11director = QgsVectorLayerDirector(
12    vectorLayer, -1, "", "", "", QgsVectorLayerDirector.DirectionBoth
13)
14strategy = QgsNetworkDistanceStrategy()
15director.addStrategy(strategy)
16
17builder = QgsGraphBuilder(vectorLayer.sourceCrs())
18
19startPoint = QgsPointXY(1179661.925139, 5419188.074362)
20endPoint = QgsPointXY(1180942.970617, 5420040.097560)
21
22tiedPoints = director.makeGraph(builder, [startPoint, endPoint])
23tStart, tStop = tiedPoints
24
25graph = builder.graph()
26idxStart = graph.findVertex(tStart)
27
28tree = QgsGraphAnalyzer.shortestTree(graph, idxStart, 0)
29
30idxStart = tree.findVertex(tStart)
31idxEnd = tree.findVertex(tStop)
32
33if idxEnd == -1:
34    raise Exception("No route!")
35
36# Add last point
37route = [tree.vertex(idxEnd).point()]
38
39# Iterate the graph
40while idxEnd != idxStart:
41    edgeIds = tree.vertex(idxEnd).incomingEdges()
42    if len(edgeIds) == 0:
43        break
44    edge = tree.edge(edgeIds[0])
45    route.insert(0, tree.vertex(edge.fromVertex()).point())
46    idxEnd = edge.fromVertex()
47
48# Display
49rb = QgsRubberBand(iface.mapCanvas())
50rb.setColor(Qt.green)
51
52# This may require coordinate transformation if project's CRS
53# is different from layer's CRS
54for p in route:
55    rb.addPoint(p)

19.2.2. 利用可能領域

頂点Aに対する利用可能領域とは、頂点Aから到達可能であり、Aからこれらの頂点までの経路のコストがある値以下になるような、グラフの頂点の部分集合です。

より明確に、これは次の例で示すことができます。「消防署があります。消防車が5分、10分、15分で到達できるのは市内のどの部分ですか？」。これらの質問への回答が消防署の利用可能領域です。

利用可能領域を見つけるには、 QgsGraphAnalyzer クラスの dijkstra() メソッドを使用します。cost配列の要素をあらかじめ定義された値と比較するだけです。もしcost[i]があらかじめ定義された値以下であれば、頂点iは利用可能領域内であり、そうでなければ領域外です。

より難しい問題は、利用可能領域の境界を取得することです。下限はまだ到達できる頂点の集合であり、上限は到達できない頂点の集合です。実際にはこれは簡単です：それは、辺の元頂点が到達可能であり辺の先頂点が到達可能ではないような、最短経路木の辺に基づく利用可能境界です。

Here is an example:

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