네트워크 분석 라이브러리

리비전 ee19294562 (QGIS 1.8 버전 이상)부터 네트워크 분석 라이브러리가 QGIS 코어의 분석 라이브러리에 추가되었습니다. 이 라이브러리의 기능은 다음과 같습니다.

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

일반 정보

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

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

  2. 그래프 분석 실행

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

그래프 만들기

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

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

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

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

빌더(builder) 프로그래밍 패턴을 사용해서 벡터 레이어를 그래프로 변환시키는 작업을 수행합니다. director라고 불리는 클래스를 사용해서 그래프를 만듭니다. 현재 director는 QgsLineVectorLayerDirector 하나밖에 없습니다. director는 builder가 라인 벡터 레이어에서 그래프를 생성할 때 쓰일 기본 설정들을 지정합니다. 현재 director와 함께 작업할 수 있는 builder는 QgsGraph 오브젝트를 생성하는 QgsGraphBuilder 하나뿐입니다. 여러분이 직접 BGL 또는 NetworkX 같은 라이브러리와 호환되는 그래프를 만들 수 있는 자신만의 builder를 구현하고 싶을 수도 있겠군요.

엣지 속성을 계산하려면 스트래티지(strategy) 프로그래밍 패턴을 사용합니다. 현재 사용할 수 있는 것은 QgsDistanceArcProperter 스트래티지뿐이며, 이 클래스는 경로의 길이를 받습니다. 모든 필요한 파라미터를 사용할 수 있도록 자신만의 스트래티지를 구현할 수도 있습니다. 예를 들면, RoadGraph 플러그인은 경계선 길이와 속성에 있는 속력 값을 이용해서 여행 시간을 계산하는 스트래티지를 사용합니다.

이제 실제로 해 봅시다.

제일 먼저, 라이브러리를 이용하기 위해 networkanalysis 모듈을 임포트해야 합니다.

from qgis.networkanalysis import *

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

# 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를 생성하려면 벡터 레이어를 넘겨줘야 하는데, 이 레이어는 그래프 구조 및 각 도로 엣지에서 허용되는 움직임(단방향 또는 양방향, 순방향 또는 역방향)에 대한 정보를 위한 자료로 쓰이게 됩니다. 다음과 같이 함수를 호출합니다.

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

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

  • vl — 그래프 생성에 사용되는 벡터 레이어

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

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

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

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

  • defaultDirection — 기본 도로 방향. directionFieldId 항목이 설정되지 않거나, 앞의 항목들에서 설정한 3가지 값이 아닐 경우에 이 값을 쓰게 됩니다. 정수형입니다. 1 은 순방향, 2 는 역방향, 3 은 양방향을 의미합니다.

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

properter = QgsDistanceArcProperter()

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

director.addProperter(properter)

이제 그래프를 생성할 builder를 사용할 수 있습니다. QgsGraphBuilder 클래스 생성자는 다음 몇 가지 인자를 받습니다.

  • crs — 사용할 좌표계. 필수적인 인자입니다.

  • otfEnabled — 실시간(on the fly) 재투영 사용 여부 결정. 기본값은 True.

  • topologyTolerance — 위상오차 허용치. 기본값은 0.

  • ellipsoidID — 사용할 타원체. 기본값은 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)

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

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

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

graph = builder.graph()

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

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

그래프 분석

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

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

최단 경로 트리는 다음과 같은 속성을 가진 방향성과 가중치가 적용된 그래프(더 정확히는 트리)입니다.

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

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

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

최단 경로 트리를 얻으려면 QgsGraphAnalyzer 클래스의 shortestTree() 또는 dijkstra() 메소드를 사용하십시오. 작업 속도가 빠르고 메모리를 좀 더 효율적으로 이용하기에 dijkstra() 메소드를 이용할 것을 추천합니다.

사용자가 최단 경로 트리를 다각적으로 검토해보고 싶다면 shortestTree() 메소드가 유용합니다. 이 메소드는 항상 (QgsGraph 클래스의) 새 그래프 오브젝트를 생성하며, 다음 3개의 변수를 받습니다.

  • source — 입력 그래프

  • startVertexIdx — 트리에 있는 포인트의 인덱스 (트리의 루트)

  • criterionNum — 사용할 경계선 속성의 개수 (0부터 시작)

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

dijkstra() 메소드의 경우 인자는 같지만 2개의 배열(array)을 반환합니다. 첫 번째 배열의 i 요소는 들어오는 엣지의 인덱스를 담고 있고, 만일 들어오는 엣지가 없을 경우 -1 값을 담게 됩니다. 두 번째 배열의 i 요소는 트리의 루트에서 버텍스 i까지의 거리를 담고 있고, 만일 루트에서 버텍스 i까지 도달할 수 없는 경우에는 DOUBLE_MAX 값을 담게 됩니다.

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

다음은 shortestTree() 메소드로 생성된 그래프를 사용해 최단 경로 트리를 표출하는 매우 간단한 코드입니다. (QGIS의 범례창에서 라인스트링 레이어를 선택하고 좌표를 사용자 데이터에 맞게 바꿔 주세요.) 경고: 이 코드는 예제로만 사용하십시오. QgsRubberBand 클래스 오브젝트를 많이 생성하기 때문에 대용량 데이터셋 상에서는 느려질 수 있습니다.

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를 찾아 포인트 B에서 포인트 A까지 트리를 따라갑니다. 이 전체 알고리즘을 다음과 같이 작성할 수 있습니다.

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

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

다음은 QGIS 파이썬 콘솔에서 shortestTree() 메소드를 사용여 알고리즘을 구현한 예시 코드입니다. (범례창에서 라인스트링 레이어를 선택하고 코드 내부의 좌표를 사용자 데이터에 맞게 바꿔야 할 것입니다.)

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)

도달 가능 범위

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

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

도달 가능 범위를 찾는 데 QgsGraphAnalyzer 클래스의 dijkstra() 메소드를 사용할 수 있습니다. 비용 배열의 요소들을 지정된 값과 비교하기만 하면 됩니다. 비용[i]가 지정된 값보다 작거나 같을 경우, vertex i는 도달 가능 범위 안에 있는 겁니다. 초과할 경우는 바깥입니다.

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

다음은 그 예시입니다.

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