Бібліотека аналізу мереж

Starting from revision ee19294562 (QGIS >= 1.8) the new network analysis library was added to the QGIS core analysis library. The library:

  • дозволяє створювати математичний граф з географічних даних (лінійних векторних шарів)

  • реалізує базові методи теорії графів (на сьогодні лише алгоритм Дейкстри)

Бібліотека аналізу мереж з’явилась як результат експорту основних функцій плаґіна RoadGraph і тепер ви можете використовувати його можливості у своїх модулях або з консолі Python.

Загальна інформація

Типовий алгоритм використання бібліотеки складається з наступних кроків:

  1. створити граф з просторових даних (зазвичай лінійних векторних шарів)

  2. проаналізувати граф

  3. використати результати аналізу (наприклад. візуалізувати їх)

Побудова графу

Перш за все необхідно підготувати вхідні дані, тобто конвертувати векторний шар у граф. Всі подальші операції будуть виконуватися з цим графом, а не з шаром.

В якості джерела даних може виступати будь-який векторний шар. Вузли ліній стануть вузлами графу, а сегменти ліній — ребрами. Якщо декілька вершин мають однакові координати, то вони будуть відповідати одній вершині графу. Тобто дві лінії, що мають спільний вузол, будуть зв’язані між собою.

Крім того, під час створення графу можна «прив’язати» до векторного шару будь-яку кількість додаткових точок. Для кожної такої додаткової точки буде знайдено відповідність — найближчий вузол чи ребро графу. В останньому випадку ребро буде розбито на дві частини, і з’явиться новий вузол.

В якості характеристик ребер можуть використовуватися атрибути векторного шару або довжина ребра.

Конвертор з векторного шару у граф реалізовано з використанням шаблону програмування будівельник. А за побудову графу відповідає так званий «директор». На сьогодні доступний лише один директор QgsLineVectorLayerDirector. Директор задає основні налаштування, які будуть використовуватися під час конвертації шару в граф, та за допомогою будівельника будує граф. Як і у випадку з директором, на сьогодні доступний лише один будівельник QgsGraphBuilder, який генерує об’єкти QgsGraph. Ви можете реалізувати своїх власних будівельників, які будуть генерувати графи, сумісні з такими бібліотеками як BGL та NetworkX.

Для обчислення характеристик ребер використовується шаблон проектування стратегія. На сьогодні доступна лише стратегія QgsDistanceArcProperter, яка враховує довжину маршруту. За необхідності ви можете реалізувати власну стратегію, яка буде використовувати все необхідні параметри. Наприклад, плаґін RoadGraph використовує стратегію, що обчислює час подорожі на основі довжини ребер та швидкості з атрибутів шару.

Розглянемо процес створення графу докладніше.

Спочатку необхідно імпортувати модуль аналізу мереж

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 — значення поля, яке відповідає прямому напрямку руху (рух від першої вершини до останньої)

  • directDirectionValue — значення поля, яке відповідає прямому напрямку руху (рух від першої вершини до останньої)

  • bothDirectionValue —значення поля, яке відповідає двосторонньому руху (тобто допускається як рух від першої вершини до останньої, так і рух в зворотньому напрямі, від останньої вершини до першої)

  • defaultDirection — напрям руху за замовчанням. Ця величина використовується для тих ділянок доріг, для яких значення поля directionFieldId не задане або не співпадає з жодним з вищенаведених.

Далі необхідно створити стратегію обчислення характеристик ребер графу

properter = QgsDistanceArcProperter()

Та повідомити директору про неї. Один директор може використовувати декілька стратегій

director.addProperter(properter)

Нарешті створюємо будівельника, який буде будувати граф. Конструктор QgsGraphBuilder приймає наступні параметри:

  • crs — система координат, що буде використовуватися. Обов’язковий параметр.

  • otfEnabled — вказує на необхідність використання перепроектування «на льоту». За замовчанням True (використовувати перепроектування).

  • topologyTolerance — топологічний допуск. За замовчанням 0.

  • ellipsoidID — еліпсоїд, який буде використовуватися. За замовчанням WGS 84.

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

Аналіз графу

В основі аналізу мереж лежить задача зв’язності графу та задача пошуку найкоротшого маршруту. Для вирішення цих задач у бібліотеці аналізу мереж реалізовано алгоритм Дейкстри.

Алгоритм Дейсктри знаходить найкоротший маршрут від однієї вершини графу до всіх інших а також значення параметру, що оптимізується. Наочно результат можна представити як дерево найкоротших маршрутів.

Дерево найкоротших маршрутів — це орієнтований зважений граф (або більш точно — дерево) з наступними властивостями:

  • тільки одна вершина не має вхідних ребер — корінь дерева

  • всі інші вершини мають лише одне вхідне ребро

  • якщо до вершини B можна дістатися з вершини A, то шлях, який їх з’єднує, єдиний і він же найкоротший (оптимальний) на вихідному графі

Дерево найкоротших шляхів можна отримати за допомогою методів shortestTree() та dijkstra() класу QgsGraphAnalyzer. Рекомендується використовувати метод dijkstra(), оскільки він працює швидше та ефективніше використовує пам’ять.

Метод shortestTree() може бути корисний, коли необхідно обійти дерево найкоротших маршрутів. Він створює новий об’єкт (завжди QgsGraph) та приймає три параметри:

  • source — вихідний граф

  • startVertexIdx — індекс точки на графі (корінь дерева)

  • criterionNum — порядковий номер характеристики ребра (відлік починається з 0).

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

Метод dijkstra() має такі ж параметри, але повертає два масиви. У першому масиві i-тий елемент містить індекс ребра, що входить в i-ту вершину, в протилежному випадку — -1. У другому масиві i-тий елемент містить відстань від i-ї вершини, якщо до вершини можна дістатися з кореня, або максимально велике значення, яке може вміститися у тип C++ double, якщо вершина недосяжна.

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

Ось дуже простий спосіб відобразити дерево найкоротших маршрутів з використанням графу, отриманого в результаті роботи метода shortestTree() (не забудьте замінити координати початкової точки на свої, а також виберіть шар доріг у легенді). Увага: використовуйте цей код тільки як приклад і лише для невеликих шарів, так як він генерує велику кількість об’єктів 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

Same thing but using dijkstra() method

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

Пошук найкоротшого маршруту

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 methods shortestTree() or dijkstra() 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 Т = 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

На цьому пошук маршруту завершено. Ми отримали інвертований список вершин (тобто вершини йдуть у зворотньому порядку, від кінцевої точки до початкової), які будуть відвідані під час руху по цьому маршруту.

Ось приклад пошуку найкоротшого маршруту для консолі Python 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)

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

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 хвилин? Відповіддю на ці питання і буде область досяжності пожежного депо.

Пошук областей досяжності легко організувати за допомогою методу dijksta() класу QgsGraphAnalyzer. Достатньо порівняти елементи масиву вартості з необхідною величиною. Якщо cost[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))