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

改訂 `ee19294562 <https://github.com/qgis/QGIS/commit/ee19294562b00c6ce957945f14c1727210cffdf7> `_ (QGIS> = 1.8)から始め、新しいネットワーク解析ライブラリがQGISコア解析ライブラリに追加されました。ライブラリ:

  • 地理データから数学的なグラフ(ポリラインベクタレイヤー)を作成します。

  • グラフ理論からの基本的なメソッドを実装します(現在はダイクストラ法のみ)

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

一般情報

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

  1. 地理データから地理学的なグラフ(たいていはポリラインベクタレイヤー)を作成します。

  2. グラフ分析の実行

  3. 分析結果の利用(例えば、これらの可視化)

グラフの構築

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

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

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

ベクターレイヤ属性と辺の長さは、辺の性質として使用できます。

ベクターレイヤーからグラフへの変換は `ビルダー<http://en.wikipedia.org/wiki/Builder_pattern>`_ プログラミングパターンを使用して行われます。グラフは、いわゆるディレクターを用いて構成されています。今のところ唯一のディレクターがあります: `QgsLineVectorLayerDirector<http://qgis.org/api/classQgsLineVectorLayerDirector.html>`_ 。ディレクターは、グラフを作成するためにビルダーによって使用される線ベクトルレイヤーからグラフを構築するために使用される基本的な設定を設定します。現在、ディレクターとの場合のように、ビルダーは一つだけ存在します: `QgsGraphBuilder<http://qgis.org/api/classQgsGraphBuilder.html>`_ 、それは`QgsGraph <http://qgis.org/api/classQgsGraph.html>`_ オブジェクトを作成します。自作のビルダーは`BGL<http://www.boost.org/doc/libs/1_48_0/libs/graph/doc/index.html>`_ か、NetworkX などのライブラリと互換性のあるグラフを構築するように実装できます。

辺性質を計算するにはプログラミングパターン `戦略<http://en.wikipedia.org/wiki/Strategy_pattern>`_ が使用されます。今のところ利用できるのは、経路の長さを考慮に入れる `QgsDistanceArcProperter<http://qgis.org/api/classQgsDistanceArcProperter.html>`_ 戦略だけです。すべての必要なパラメータを使用する独自の戦略を実装できます。例えば、RoadGraphプラグインは、属性から辺長と速度値を使用して旅行時間を計算する戦略を使用します。

プロセスに飛び込む時間です。

まず第一に、このライブラリを使用するために、networkanalysisモジュールをインポートする必要があります

from qgis.networkanalysis import *

ディレクターを作成するためのその後いくつかの例

# don't use information about road direction from layer attributes,
# all roads are treated as two-way
director = QgsLineVectorLayerDirector(vLayer, -1, '', '', '', 3)

# use field with index 5 as source of information about road direction.
# one-way roads with direct direction have attribute value "yes",
# one-way roads with reverse direction have the value "1", and accordingly
# bidirectional roads have "no". By default roads are treated as two-way.
# This scheme can be used with OpenStreetMap data
director = QgsLineVectorLayerDirector(vLayer, 5, 'yes', '1', 'no', 3)

ディレクタを構築するために、ベクターレイヤーを渡さなければなりません。これはグラフ構造および各道路セグメント上で許される移動(一方向または双方向の動き、順または逆の方向)についての情報のソースとして使用されるでしょう。呼び出しは次のようになります

director = QgsLineVectorLayerDirector(vl, directionFieldId,
                                      directDirectionValue,
                                      reverseDirectionValue,
                                      bothDirectionValue,
                                      defaultDirection)

そして、ここでこれらのパラメータは何を意味するかの完全なリストは以下のとおりです。

  • vl —グラフを構築するために使用される ベクターレイヤー

  • directionFieldId —道路の方向に関する情報が格納されている属性テーブルのフィールドのインデックス。`` -1``、その後、まったくこの情報を使用しない場合。整数。

  • directDirectionValue —順方向(線の最初の点から最後の点へ移動)の道に対するフィールド値。文字列。

  • reverseDirectionValue —逆方向(線の最後の点から最初の点へ移動)の道に対するフィールド値。文字列。

  • bothDirectionValue —双方向道路に対するフィールド値(例えば最初の点から最後まで、また最後から最初まで移動できる道に対する)。文字列。

  • defaultDirection —デフォルトの道路の方向。この値は、フィールド directionFieldId```が設定されていない、または上で指定された3つの値のいずれとも異なる値を有する道路に使用されるであろう。整数。``1 は順方向、 2 は逆方向、 3 は両方向を示します。

それから、辺性質を計算するための戦略を作成することが必要です

properter = QgsDistanceArcProperter()

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

director.addProperter(properter)

これで、グラフを作成するビルダーを使用できます。QgsGraphBuilderクラスのコンストラクタには、いくつかの引数を取ります。

  • crs —使用する空間参照系。必須の引数です。

  • otfEnabled — “その場で” 再投影を使用するかどうか。デフォルトでは True (OTFを使用).

  • トポロジ許容値—トポロジ的な許容値です。デフォルト値は0です。

  • 楕円体ID — 使用する楕円体です。デフォルトは “WGS84”です。

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

分析に使用されるいくつかのポイントを、定義することもできます。例えば

startPoint = QgsPoint(82.7112, 55.1672)
endPoint = QgsPoint(83.1879, 54.7079)

これですべてが整いましたので、グラフを構築し、それにこれらの点を「結びつける」ことができます。

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

グラフを構築するには(レイヤー中の図形数とレイヤーのサイズのに応じて)いくらか時間がかかることがあります。 tiedPoints は「結ばれた」点の座標とのリストです。ビルド操作が完了すると、グラフは取得して分析のために使用できます

graph = builder.graph()

次のコードで、点の頂点インデックスを取得できます

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

グラフ分析

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

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

最短経路木は、次の性質を有する有向重み付きグラフ(より正確には、木)です。

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

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

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

最短経路木を得るには、 `QgsGraphAnalyzer<http://qgis.org/api/classQgsGraphAnalyzer.html>`_ クラスの shortestTree()dijkstra() メソッドを使用してください。より速く動作し、より多くのメモリを効率的に使用するため、 メソッド dijkstra() を使用することをお勧めします。

最短経路木の周りを歩くしたい場合 shortestTree() メソッドが便利です。それは、常に新しいグラフオブジェクト(QgsGraph)を作成し、三つの変数を受けとります。

  • source — 入力グラフ

  • startVertexIdx — ツリー上のポイントのインデックス(ツリーのルート)

  • criterionNum — 使用するエッジプロパティの数(0から始まる)。

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

dijkstra() メソッドは引数は同じですが2つの配列を返します。流入辺が存在しない場合は、最初の配列要素iに流入辺のインデックスまたは-1を含有します。二番目の配列要素iには木の根から頂点iまでの距離が入ります。頂点iが根から到達不能である場合はDOUBLE_MAXが入ります。

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

shortestTree() メソッド(TOCでラインストリングのレイヤーを選択し、自分で座標を置き換える)ここで作成したグラフを使用して最短パスツリーを表示するには、いくつかの非常に単純なコードです。 警告 :このコードはほんの一例として使用してください、それは多くの `QgsRubberBand<http://qgis.org/api/classQgsRubberBand.html>`_ オブジェクトを作成し、大規模なデータセットでは低速になることがあります。

from PyQt4.QtCore import *
from PyQt4.QtGui import *

from qgis.core import *
from qgis.gui import *
from qgis.networkanalysis import *

vl = qgis.utils.iface.mapCanvas().currentLayer()
director = QgsLineVectorLayerDirector(vl, -1, '', '', '', 3)
properter = QgsDistanceArcProperter()
director.addProperter(properter)
crs = qgis.utils.iface.mapCanvas().mapRenderer().destinationCrs()
builder = QgsGraphBuilder(crs)

pStart = QgsPoint(-0.743804, 0.22954)
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.arcCount()):
  rb = QgsRubberBand(qgis.utils.iface.mapCanvas())
  rb.setColor (Qt.red)
  rb.addPoint (tree.vertex(tree.arc(i).inVertex()).point())
  rb.addPoint (tree.vertex(tree.arc(i).outVertex()).point())
  i = i + 1

同じですが dijkstra() メソッドを使用して

from PyQt4.QtCore import *
from PyQt4.QtGui import *

from qgis.core import *
from qgis.gui import *
from qgis.networkanalysis import *

vl = qgis.utils.iface.mapCanvas().currentLayer()
director = QgsLineVectorLayerDirector(vl, -1, '', '', '', 3)
properter = QgsDistanceArcProperter()
director.addProperter(properter)
crs = qgis.utils.iface.mapCanvas().mapRenderer().destinationCrs()
builder = QgsGraphBuilder(crs)

pStart = QgsPoint(-1.37144, 0.543836)
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(qgis.utils.iface.mapCanvas())
  rb.setColor (Qt.red)
  rb.addPoint (graph.vertex(graph.arc(edgeId).inVertex()).point())
  rb.addPoint (graph.vertex(graph.arc(edgeId).outVertex()).point())

最短経路を見つける

2点間の最適経路を見つけるために下のアプローチが使用されます。両方の点(始点Aおよび終点B)は構築されたときにグラフに「結ばれます」。それから、メソッド shortestTree() または dijkstra() を使用して、開始点Aに根を持つ最短経路木を構築します。また、同じ木の中に終点Bを見つけ、木を点Aから点Bまで歩き始めます。全体のアルゴリズムは次のように書けます

assign Т = B
while Т != A
    add point Т to path
    get incoming edge for point Т
    look for point ТТ, that is start point of this edge
    assign Т = ТТ
add point А to path

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

これはメソッド shortestTree() を使用するQGIS Pythonコンソールのためのサンプルコードです(TOC中のラインストリングレイヤーを選択して、コードの座標を自分のものに置き換える必要があります)

from PyQt4.QtCore import *
from PyQt4.QtGui import *

from qgis.core import *
from qgis.gui import *
from qgis.networkanalysis import *

vl = qgis.utils.iface.mapCanvas().currentLayer()
director = QgsLineVectorLayerDirector(vl, -1, '', '', '', 3)
properter = QgsDistanceArcProperter()
director.addProperter(properter)
crs = qgis.utils.iface.mapCanvas().mapRenderer().destinationCrs()
builder = QgsGraphBuilder(crs)

pStart = QgsPoint(-0.835953, 0.15679)
pStop = QgsPoint(-1.1027, 0.699986)

tiedPoints = director.makeGraph(builder, [pStart, pStop])
graph = builder.graph()

tStart = tiedPoints[0]
tStop = tiedPoints[1]

idStart = graph.findVertex(tStart)
tree = QgsGraphAnalyzer.shortestTree(graph, idStart, 0)

idStart = tree.findVertex(tStart)
idStop = tree.findVertex(tStop)

if idStop == -1:
  print "Path not found"
else:
  p = []
  while (idStart != idStop):
    l = tree.vertex(idStop).inArc()
    if len(l) == 0:
      break
    e = tree.arc(l[0])
    p.insert(0, tree.vertex(e.inVertex()).point())
    idStop = e.outVertex()

  p.insert(0, tStart)
  rb = QgsRubberBand(qgis.utils.iface.mapCanvas())
  rb.setColor(Qt.red)

  for pnt in p:
    rb.addPoint(pnt)

そしてこれは同じサンプルですが dijkstra() メソッドを使用しています

from PyQt4.QtCore import *
from PyQt4.QtGui import *

from qgis.core import *
from qgis.gui import *
from qgis.networkanalysis import *

vl = qgis.utils.iface.mapCanvas().currentLayer()
director = QgsLineVectorLayerDirector(vl, -1, '', '', '', 3)
properter = QgsDistanceArcProperter()
director.addProperter(properter)
crs = qgis.utils.iface.mapCanvas().mapRenderer().destinationCrs()
builder = QgsGraphBuilder(crs)

pStart = QgsPoint(-0.835953, 0.15679)
pStop = QgsPoint(-1.1027, 0.699986)

tiedPoints = director.makeGraph(builder, [pStart, pStop])
graph = builder.graph()

tStart = tiedPoints[0]
tStop = tiedPoints[1]

idStart = graph.findVertex(tStart)
idStop = graph.findVertex(tStop)

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

if tree[idStop] == -1:
  print "Path not found"
else:
  p = []
  curPos = idStop
  while curPos != idStart:
    p.append(graph.vertex(graph.arc(tree[curPos]).inVertex()).point())
    curPos = graph.arc(tree[curPos]).outVertex();

  p.append(tStart)

  rb = QgsRubberBand(qgis.utils.iface.mapCanvas())
  rb.setColor(Qt.red)

  for pnt in p:
    rb.addPoint(pnt)

利用可能領域

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

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

利用可能領域を見つけるためには、 QgsGraphAnalyzer クラスのメソッド dijkstra() を使用できます。事前に定義された値とコスト配列の要素を比較すれば十分です。コスト[i]が所定値以下の場合は、頂点iは利用可能領域内に、そうでない場合は利用可能領域外にあります。

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

例があります

from PyQt4.QtCore import *
from PyQt4.QtGui import *

from qgis.core import *
from qgis.gui import *
from qgis.networkanalysis import *

vl = qgis.utils.iface.mapCanvas().currentLayer()
director = QgsLineVectorLayerDirector(vl, -1, '', '', '', 3)
properter = QgsDistanceArcProperter()
director.addProperter(properter)
crs = qgis.utils.iface.mapCanvas().mapRenderer().destinationCrs()
builder = QgsGraphBuilder(crs)

pStart = QgsPoint(65.5462, 57.1509)
delta = qgis.utils.iface.mapCanvas().getCoordinateTransform().mapUnitsPerPixel() * 1

rb = QgsRubberBand(qgis.utils.iface.mapCanvas(), True)
rb.setColor(Qt.green)
rb.addPoint(QgsPoint(pStart.x() - delta, pStart.y() - delta))
rb.addPoint(QgsPoint(pStart.x() + delta, pStart.y() - delta))
rb.addPoint(QgsPoint(pStart.x() + delta, pStart.y() + delta))
rb.addPoint(QgsPoint(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 = 2000.0
i = 0
while i < len(cost):
  if cost[i] > r and tree[i] != -1:
    outVertexId = graph.arc(tree [i]).outVertex()
    if cost[outVertexId] < r:
      upperBound.append(i)
  i = i + 1

for i in upperBound:
  centerPoint = graph.vertex(i).point()
  rb = QgsRubberBand(qgis.utils.iface.mapCanvas(), True)
  rb.setColor(Qt.red)
  rb.addPoint(QgsPoint(centerPoint.x() - delta, centerPoint.y() - delta))
  rb.addPoint(QgsPoint(centerPoint.x() + delta, centerPoint.y() - delta))
  rb.addPoint(QgsPoint(centerPoint.x() + delta, centerPoint.y() + delta))
  rb.addPoint(QgsPoint(centerPoint.x() - delta, centerPoint.y() + delta))