Începând cu revizia ee19294562 (QGIS >= 1.8) noua bibliotecă de analiză de rețea a fost adaugată la biblioteca de analize de bază din QGIS. Biblioteca:
creează graful matematic din datele geografice (straturi vectoriale de tip polilinie)
implementează metode de bază din teoria grafurilor (în prezent, doar algoritmul lui Dijkstra)
Biblioteca analizelor de rețea a fost creată prin exportarea funcțiilor de bază ale plugin-ului RoadGraph, iar acum aveți posibilitatea să-i utilizați metodele în plugin-uri sau direct în consola Python.
Pe scurt, un caz tipic de utilizare poate fi descris astfel:
crearea grafului din geodate (de obicei un strat vectorial de tip polilinie)
rularea analizei grafului
folosirea rezultatelor analizei (de exemplu, vizualizarea lor)
Primul lucru pe care trebuie să-l faceți — este de a pregăti datele de intrare, ceea ce înseamnă conversia stratului vectorial într-un graf. Toate acțiunile viitoare vor folosi acest graf, și nu stratul.
Ca și sursă putem folosi orice strat vectorial de tip polilinie. Nodurile poliliniilor devin noduri ale grafului, segmentele poliliniilor reprezentând marginile grafului. În cazul în care mai multe noduri au aceleași coordonate, atunci ele sunt în același nod al grafului. Astfel, două linii care au un nod comun devin conectate între ele.
În plus, în timpul creării grafului este posibilă “fixarea” (‘legarea”) de stratul vectorial de intrare a oricărui număr de puncte suplimentare. Pentru fiecare punct suplimentar va fi găsită o potrivire — cel mai apropiat nod sau cea mai apropiată muchie a grafului. În ultimul caz muchia va fi divizată iar noul nod va fi adăugat.
Atributele stratului vectorial și lungimea unei muchii pot fi folosite ca proprietăți ale marginii.
Convertorul din strat vectorial în graf este dezvoltat folosind modelul de programare al Constructorului. De construirea grafului răspunde așa-numitul Director. Există doar un singur Director pentru moment: QgsLineVectorLayerDirector. Directorul stabilește setările de bază ale Constructorului, care vor fi folosite pentru a construi un graf dintr-un strat vectorial de tip linie. În prezent, ca și în cazul Directorului, există doar un singur Constructor: QgsGraphBuilder, care creează obiecte <http://qgis.org/api/classQgsGraph.html>`_. Este posibil să doriți implementarea propriilor constructori care să construiască grafuri compatibile cu bibliotecile, cum ar fi BGL sau NetworkX.
Pentru a calcula proprietățile marginii este utilizată strategia modelului de programare. Pentru moment doar strategia QgsDistanceArcProperter este disponibilă, care ia în considerare lungimea traseului. Puteți implementa propria strategie, care va folosi toți parametrii necesari. De exemplu, plugin-ul RoadGraph folosește strategia care calculează timpul de călătorie, folosind lungimea muchiei și valoarea vitezei din atribute.
Este timpul de a aprofunda acest proces.
Înainte de toate, pentru a utiliza această bibliotecă ar trebui să importăm modulul networkanalysis
from qgis.networkanalysis import *
Apoi, câteva exemple pentru crearea unui 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)
Pentru a construi un director ar trebui să transmitem stratul vectorial, care va fi folosit ca sursă pentru structura grafului și informațiile despre mișcările permise pe fiecare segment de drum (circulație unilaterală sau bilaterală, sens direct sau invers). Acest apel arată în felul următor
director = QgsLineVectorLayerDirector(vl, directionFieldId,
directDirectionValue,
reverseDirectionValue,
bothDirectionValue,
defaultDirection)
Iată lista completă a ceea ce înseamnă acești parametri:
vl — stratul vectorial utilizat pentru a construi graf
directionFieldId — indexul câmpului din tabelul de atribute, în care sunt stocate informații despre direcțiile drumurilor. Dacă este -1, atunci aceste informații nu se folosesc deloc. Număr întreg.
directDirectionValue — valoarea câmpului pentru drumurile cu sens direct (trecere de la primul punct de linie la ultimul). Șir de caractere.
reverseDirectionValue — valoarea câmpului pentru drumurile cu sens invers (în mișcare de la ultimul punct al liniei până la primul). Șir de caractere.
bothDirectionValue — valoarea câmpului pentru drumurile bilaterale (pentru astfel de drumuri putem trece de la primul la ultimul punct și de la ultimul la primul). Șir de caractere.
defaultDirection — direcția implicită a drumului. Această valoare va fi folosită pentru acele drumuri în care câmpul directionFieldId nu este setat sau are o valoare diferită de oricare din cele trei valori specificate mai sus. Număr întreg. 1 indică sensul direct, 2 indică sensul inversă, iar 3 indică ambele sensuri.
Este necesară, apoi, crearea unei strategii pentru calcularea proprietăților marginii
properter = QgsDistanceArcProperter()
Apoi spuneți directorului despre această strategie
director.addProperter(properter)
Acum putem crea constructorul, care va crea graful. Constructorul clasei QgsGraphBuilder ia mai multe argumente:
crs — sistemul de coordonate de referință de utilizat. Argument obligatoriu.
otfEnabled — utilizați sau nu reproiectarea “din zbor”. În mod implicit const:True (folosiți OTF).
topologyTolerance — toleranța topologică. Valoarea implicită este 0.
ellipsoidID — elipsoidul de utilizat. În mod implicit “WGS84”.
# only CRS is set, all other values are defaults
builder = QgsGraphBuilder(myCRS)
De asemenea, putem defini mai multe puncte, care vor fi utilizate în analiză. De exemplu
startPoint = QgsPoint(82.7112, 55.1672)
endPoint = QgsPoint(83.1879, 54.7079)
Acum că totul este la locul lui, putem să construim graful și să “legăm” aceste puncte la el
tiedPoints = director.makeGraph(builder, [startPoint, endPoint])
Construirea unui graf poate dura ceva timp (depinzând de numărul de entități dintr-un strat și de dimensiunea stratului). tiedPoints reprezintă o listă cu coordonatele punctelor “asociate”. Când s-a terminat operațiunea de construire putem obține graful și să-l utilizăm pentru analiză
graph = builder.graph()
Cu următorul cod putem obține indecșii punctelor noastre
startId = graph.findVertex(tiedPoints[0])
endId = graph.findVertex(tiedPoints[1])
Analiza de rețea este utilizată pentru a găsi răspunsuri la două întrebări: care noduri sunt conectate și identificarea celei mai scurte căi. Pentru a rezolva această problemă, biblioteca de analiză de rețea oferă algoritmul lui Dijkstra.
Algoritmul lui Dijkstra găsește cea mai bună cale între unul dintre vârfurile grafului și toate celelalte, precum și valorile parametrilor de optimizare. Rezultatele pot fi reprezentate ca cel mai scurt arbore.
Arborele drumurilor cele mai scurte reprezintă un graf (sau mai precis — arbore) orientat, ponderat, cu următoarele proprietăți:
doar un singur nod nu are muchii de intrare — rădăcina arborelui
toate celelalte noduri au numai o margine de intrare
dacă nodul B este accesibil din nodul A, apoi calea de la A la B este singura disponibilă și este optimă (cea mai scurtă) în acest graf
Pentru a obține cel mai scurt arbore folosiți metodele shortestTree() și dijkstra() ale clasei QgsGraphAnalyzer. Se recomandă utilizarea metodei dijkstra(), deoarece lucrează mai rapid și utilizează memoria mult mai eficient.
Metoda shortestTree() este utilă atunci când doriți să vă plimbați de-a lungul celei mai scurte căi. Aceasta creează mereu un nou obiect de tip graf (QgsGraph) care acceptă trei variabile:
source — graf de intrare
startVertexIdx — Indexul punctului de pe arbore (rădăcina arborelui)
criterionNum — numărul de proprietății marginii de folosit (începând de la 0).
tree = QgsGraphAnalyzer.shortestTree(graph, startId, 0)
The dijkstra() method has the same arguments, but returns two arrays. In the first array element i contains index of the incoming edge or -1 if there are no incoming edges. In the second array element i contains distance from the root of the tree to vertex i or DOUBLE_MAX if vertex i is unreachable from the root.
(tree, cost) = QgsGraphAnalyzer.dijkstra(graph, startId, 0)
Iată un cod foarte simplu pentru a afișa arborele celei mai scurte căi, folosind graful creat cu metoda shortestTree() (selectați stratul linie în TOC și înlocuiți coordonatele cu ale dvs). Atenție: folosiți acest cod doar ca exemplu, deoarece el va crea o mulțime de obiecte QgsRubberBand , putând fi lent pe seturi de date de mari dimensiuni.
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
Același lucru, dar cu ajutorul metodei 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())
Pentru a găsi calea optimă între două puncte este utilizată următoarea abordare. Ambele puncte (se începe cu A și se termină cu B) sunt “legate” de un graf, atunci când se construiește. Apoi folosind metodele shortestTree() sau dijkstra() vom construi cel mai scurt arbore cu rădăcina în punctul de pornire A. În același arbore am găsit, de asemenea, punctul de final B și începem parcurgerea arborelui de la punctul B la punctul A. Întregul algoritm poate fi scris ca
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
În acest moment avem calea, sub formă de listă inversată de noduri (nodurile sunt listate în ordine inversă, de la punctul de final către cel de start), ele fiind vizitate în timpul parcurgerii căii.
Aici este codul de test pentru consola Python a QGIS (va trebui să selectați stratul linie în TOC și să înlocuiți coordonate din cod cu ale dvs.), care folosește metoda 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)
Iar aici este același exemplu, dar folosind metoda 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)
Aria de disponibilitate a nodului A este un subset de noduri ale graf-ului, care sunt accesibile din nodul A iar costurile căii de la A la aceste noduri nu sunt mai mari decât o anumită valoare.
Mai clar, acest lucru poate fi dovedit cu următorul exemplu: “Există o echipă de intervenție în caz de incendiu. Ce zone ale orașului acoperă această echipă în 5 minute? Dar în 10 minute? Dar în 15 minute?”. Răspunsul la aceste întrebări îl reprezintă zonele de disponibilitate ale echipei de intervenție.
Pentru a găsi zonele de disponibilitate putem folosi metoda dijksta() a clasei QgsGraphAnalyzer. Este suficientă compararea elementelor matricei de costuri cu valoarea predefinită. În cazul în care costul[i] este mai mic sau egal decât o valoare predefinită, atunci nodul i se află în zona de disponibilitate, în caz contrar este în afară.
Mai dificilă este obținerea granițelor zonei de disponibilitate. Marginea de jos reprezintă un set de noduri care încă sunt accesibile, iar marginea de sus un set de noduri inaccesibile. De fapt, acest lucru este simplu: marginea disponibilă a atins aceste margini parcurgând arborele cel mai scurt, pentru care nodul de start este accesibil, spre deosebire de celelalt capăt, care nu este accesibil.
Iată un exemplu
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))