19. 네트워크 분석 라이브러리

힌트

여러분이 PyQGIS 콘솔을 사용하지 않는 경우 이 페이지에 있는 코드 조각들을 다음과 같이 가져와야 합니다:

from qgis.core import (
  QgsVectorLayer,
  QgsPointXY,
)

The network analysis library can be used to:

• create mathematical graph from geographical data (polyline vector layers)

• implement basic methods from graph theory (currently only Dijkstra’s algorithm)

네트워크 분석 라이브러리는 코어 플러그인 RoadGraph 에서 기본 함수들을 가져와 생성됐습니다. 이제 플러그인에서나 파이썬 콘솔에서 직접 이 라이브러리의 메소드를 사용할 수 있습니다.

19.1. 일반 정보

이 라이브러리의 전형적인 용도를 간단히 설명하면 다음과 같습니다.

1. 공간 데이터(일반적으로 폴리라인 벡터 레이어)로부터 그래프 생성

2. 그래프 분석 실행

3. 분석 결과 이용(예를 들어 분석 결과의 시각화 등)

19.2. 그래프 만들기

가장 먼저 해야 할 일은 입력 데이터를 준비하는 것인데, 벡터 레이어를 그래프(수학적으로 간략화된 연결관계)로 변환하는 것입니다. 이후의 모든 작업은 레이어가 아니라 이 그래프를 사용합니다.

어떤 폴리라인 벡터 레이어라도 소스로 사용할 수 있습니다. 폴리라인의 노드(node)는 그래프의 버텍스(vertex)가 되고, 폴리라인의 선분(segment)은 그래프의 엣지(edge)가 됩니다. 노드 몇 개가 동일한 좌표에 있을 경우 그 노드들은 동일한 그래프 버텍스가 됩니다. 따라서 공통 노드를 가진 2개의 선분은 서로 연결됩니다.

또, 그래프 생성 중에 입력 벡터 레이어에 추가적인 포인트를 몇 개라도 “고정”(다른 용어로는 “결속”) 시킬 수 있습니다. 각 추가 포인트에 대응하는, 가장 가까운 그래프 버텍스 또는 가장 가까운 그래프 엣지를 찾을 것입니다. 후자의 경우 엣지가 나뉘어 새 버텍스가 추가됩니다.

벡터 레이어의 속성과 엣지 길이를 그래프 엣지의 속성으로 쓸 수 있습니다.

Converting from a vector layer to the graph is done using the Builder programming pattern. A graph is constructed using a so-called Director. There is only one Director for now: QgsVectorLayerDirector. The director sets the basic settings that will be used to construct a graph from a line vector layer, used by the builder to create the graph. Currently, as in the case with the director, only one builder exists: QgsGraphBuilder, that creates QgsGraph objects. You may want to implement your own builders that will build a graph compatible with such libraries as BGL or 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 availabile. You can implement your own strategy that will use all necessary parameters. For example, RoadGraph plugin uses a strategy that computes travel time using edge length and speed value from attributes.

이제 실제로 해 봅시다.

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

from qgis.analysis import *

다음은 director를 생성하는 몇 가지 방법의 예시입니다.

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

# 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 = QgsVectorLayerDirector(vectorLayer, 5, 'yes', '1', 'no', QgsVectorLayerDirector.DirectionBoth)

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

director = QgsVectorLayerDirector(vectorLayer,
                                  directionFieldId,
                                  directDirectionValue,
                                  reverseDirectionValue,
                                  bothDirectionValue,
                                  defaultDirection)

각각의 파라미터가 의미하는 바는 다음과 같습니다.

• vectorLayer — vector layer used to build the graph

• directionFieldId — 도로 방향에 관한 정보가 저장된 속성 테이블 필드의 인덱스. 값이 -1 일 경우 방향 정보를 전혀 사용하지 않습니다. 정수형입니다.

• directDirectionValue — 순방향(첫 번째 라인 포인트에서 마지막 라인 포인트로 이동)인 도로의 필드값. 문자열입니다.

• reverseDirectionValue — 역방향(마지막 라인 포인트에서 첫 번째 포인트로 이동)인 도로의 필드값. 문자열입니다.

• bothDirectionValue — 양방향(첫 번째 포인트에서 마지막으로도 마지막에서 첫 번째로도 이동 가능)인 도로의 필드값. 문자열입니다.

• defaultDirection — default road direction. This value will be used for those roads where field directionFieldId is not set or has some value different from any of the three values specified above. Possible values are:

  • QgsVectorLayerDirector.DirectionForward — One-way direct

  • QgsVectorLayerDirector.DirectionBackward — One-way reverse

  • QgsVectorLayerDirector.DirectionBoth — Two-way

그 다음 edge 속성을 계산하기 위해 strategy를 생성해야 합니다.

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

그리고 drirector에게 이 strategy에 대해 알려줍니다.

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

Now we can use the builder, which will create the graph. The QgsGraphBuilder class constructor takes several arguments:

• crs — coordinate reference system to use. Mandatory argument.

• otfEnabled — use “on the fly” reprojection or no. By default True (use OTF).

• topologyTolerance — topological tolerance. Default value is 0.

• ellipsoidID — ellipsoid to use. By default “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)

이제 모든 준비가 끝났으므로 그래프를 만들고 이 포인트들을 그래프에 “결속(tie)”시킬 수 있습니다.

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

그래프를 만드는 데 시간이 좀 걸릴 수도 있습니다. (레이어에 있는 피처의 개수 및 레이어 크기에 따라 다릅니다.) tiedPoints 는 “결속”된 포인트들의 좌표 목록입니다. builder의 작업이 완료되면 분석에 이용할 수 있는 그래프를 얻게 됩니다.

graph = builder.graph()

다음 코드를 이용하면 포인트들의 vertex 인덱스들을 얻을 수 있습니다.

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

19.3. 그래프 분석

네트워크 분석은 다음 두 가지 질문에 대한 답을 찾는 데 사용됩니다. 어떤 vertex들이 연결되어 있는가? 그리고 어떻게 최단 경로를 찾을 것인가? 네트워크 분석 라이브러리는 이 문제를 해결하기 위해 데이크스트라 알고리즘(Dijkstra’s algorithm)을 제공합니다.

데이크스트라 알고리즘은 그래프의 한 vertex에서 다른 모든 vertex로 가는 최단 경로와 최적화 파라미터의 값을 찾습니다. 그 결과는 최단 경로 트리로 나타낼 수 있습니다.

The shortest path tree is a directed weighted graph (or more precisely a tree) with the following properties:

• 들어오는 edge가 없는 vertex는 단 하나, 트리의 루트뿐입니다.

• 다른 모든 vertex는 들어오는 edge를 딱 하나 가지고 있습니다.

• vertex A에서 vertex B에 도달할 수 있다면, A에서 B로의 경로는 사용할 수 있는 단 하나의 경로이며 이 그래프에서 최적(최단) 경로입니다.

To get the shortest path tree use the methods shortestTree() and dijkstra() of the QgsGraphAnalyzer class. It is recommended to use the dijkstra() method because it works faster and uses memory more efficiently.

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 — input graph

• startVertexIdx — index of the point on the tree (the root of the tree)

• criterionNum — number of edge property to use (started from 0).

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

The dijkstra() method has the same arguments, but returns two arrays. In the first array element n contains index of the incoming edge or -1 if there are no incoming edges. In the second array element n contains the distance from the root of the tree to vertex n or DOUBLE_MAX if vertex n is unreachable from the root.

(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() method (select linestring layer in Layers panel and replace coordinates with your own).

경고

Use this code only as an example, it creates a lot of QgsRubberBand objects and may be slow on large datasets.

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

Same thing but using the dijkstra() method

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. 최단 경로 탐색

To find the optimal path between two points the following approach is used. Both points (start A and end B) are “tied” to the graph when it is built. Then using the shortestTree() or dijkstra() method we build the shortest path tree with root in the start point A. In the same tree we also find the end point B and start to walk through the tree from point B to point A. The whole algorithm can be written as:

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

이 시점에서 이 경로를 지나가는 동안 거치게 될 vertex의 역순 목록의 형태로 경로를 얻게 됩니다. (vertex들이 종료점에서 시작점의 순서로 역순으로 나열됩니다.)

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() method

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)

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!')

# Add last point
route = [tree.vertex(idxEnd).point()]

# Iterate the graph
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()

# Display
rb = QgsRubberBand(iface.mapCanvas())
rb.setColor(Qt.green)

# This may require coordinate transformation if project's CRS
# is different than layer's CRS
for p in route:
    rb.addPoint(p)

And here is the same sample but using the dijkstra() method

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]

# Add last point
route = [graph.vertex(idxEnd).point()]

# Iterate the graph
while idxEnd != idxStart:
    idxEnd = graph.edge(tree[idxEnd]).fromVertex()
    route.insert(0, graph.vertex(idxEnd).point())

# Display
rb = QgsRubberBand(iface.mapCanvas())
rb.setColor(Qt.red)

# This may require coordinate transformation if project's CRS
# is different than layer's CRS
for p in route:
    rb.addPoint(p)

19.3.2. 도달 가능 범위

vertex A의 도달 가능 범위(area of availability)란 vertex A에서 접근할 수 있고, vertex A에서 이 vertex들까지의 경로 비용이 지정된 값을 초과하지 않는, 그래프 vertex들의 부분집합을 말합니다.

다음 질문을 통해 이를 더 명확히 알 수 있습니다. “소방서가 있다. 소방차가 5분/10분/15분 안에 도착할 수 있는 도시의 구역은 어디인가?” 이 질문에 대한 답이 바로 소방서의 도달 가능 범위입니다.

To find the areas of availability we can use the dijkstra() method of the QgsGraphAnalyzer class. It is enough to compare the elements of the cost array with a predefined value. If cost[i] is less than or equal to a predefined value, then vertex i is inside the area of availability, otherwise it is outside.

도달 가능 범위의 경계를 구하는 일은 좀 더 어려운 문제입니다. 하단 경계는 도달 가능한 vertex들의 집합이고, 상단 경계는 도달 불가능한 vertex들의 집합입니다. 사실 단순합니다. 도달 가능 범위의 경계는 edge의 윈본 vertex가 접근 가능한 vertex이고, edge의 대상 vertex가 접근 불가능한 vertex인 최단경로 트리의 edge들에 기반한 경계선입니다.

다음은 그 예시입니다.

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))