19. Netzwerkanalysebibliothek

Hinweis

Die Codefragmente auf dieser Seite müssen wie folgt importiert werden, wenn Sie sich außerhalb der pyqgis-Konsole befinden:

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)

Die Netzwerkanalysebibliothek wurde erstellt, indem grundlegende Funktionen aus dem RoadGraph-Kern-Plugin exportiert wurden. Jetzt können Sie die Methoden in Plugins oder direkt von der Python-Konsole aus verwenden.

19.1. Allgemeine Informationen

Ein typischer Anwendungsfall sieht so aus:

1. Graphen aus Geodaten erstellen (Typischerweise Polylinien Vektorlayer)

2. Graphenanalyse durchführen

3. Analyseergebnisse verwenden (z.B. für Visualisierungszwecke)

19.2. Einen Graphen erstellen

Als erstes müssen die Eingabedaten vorbereitet werden, d.h. ein Vektorlayer muss in einen Graphen konvertiert werden. Alle weiteren Aktionen verwenden diesen Graphen und nicht den Vektorlayer.

Als Quelle können wir jeden Polylinienvektorlayer verwenden. Vertices der Polylinien werden zu Knoten des Graphen und Segmente der Polylinien werden zu Kanten des Graphen. Wenn mehrere Vertices die gleichen Koordinaten haben, dann repräsentieren sie den selben Knoten. So werden zwei Kanten, die einen gemeinsamen Knoten haben, miteinander verbunden.

Darüber hinaus ist es während der Graphenerstellung möglich, eine beliebige Anzahl zusätzlicher Punkte auf den Eingabevektorlayer zu „fixieren“ (engl. „tie“). Für jeden zusätzlichen Punkt wird eine Übereinstimmung gefunden, entweder der nächstgelegene Knoten des Graphen oder die nächstgelegene Kante. Im letzteren Fall wird die Kante geteilt und ein neuer Knoten hinzugefügt.

Vektorlayerattribute und die Länge eines Segements können als Eigenschaften einer Kante verwendet werden.

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.

Es ist Zeit, den Prozess näher zu beleuchten.

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

from qgis.analysis import *

Ein paar Beispiele einen director zu erzeugen

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

Folgende Liste fasst die Bedeutung der Parameter zusammen:

• vectorLayer — vector layer used to build the graph

• directionFieldId — Index des Feldes der Attributtabelle, in dem Informationen zur Straßenrichtung gespeichert sind. Falls `` -1`` wird diese Informationen nicht benutzt. Ganzzahl (integer).

• directDirectionValue — Feldwert für Straßen mit direkter Richtung (vom ersten zum letzten Punkt). Zeichenkette (string).

• reverseDirectionValue — Feldwert für Straßen mit umgekehrter Richtung (vom letzten zum ersten Punkt). Zeichenkette (string).

• bothDirectionValue — Feldwert für bidirektionale Straßen (in beide Richtungen befahrbar, vom ersten Punkt zum letzten und vom letzten zum ersten). Zeichenkette (string).

• 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

Nun muss eine Strategie für die Ermittlung der Kanteneigenschaften festgelegt werden

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

und der Director über diese Strategie informiert werden

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

Zusätzlich könne verschiedene Punkte definiert werden, die in der Analyse verwendet werden sollen, z.B.

startPoint = QgsPointXY(1179720.1871, 5419067.3507)
endPoint = QgsPointXY(1180616.0205, 5419745.7839)

Nun ist alles bereit um den Graphen erstellen zu können und diese Punkte an den Graphen zu „binden“

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

Das Erstellen des Graphen kann einige Zeit dauern (abhängig von der Anzahl der Features der Eingabelayers und der Layergröße). `` tiedPoints`` ist eine Liste mit Koordinaten von „angebundenen“ Punkten. Wenn der Build-Vorgang abgeschlossen ist, kann der Graph abgerufen und für die Analyse verwendet werden.

graph = builder.graph()

MIt diesem Code erhalten wir die Indizes unserer Knoten

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

19.3. Graphenanalyse

Die Netzwerkanalyse wird verwendet, um Antworten auf zwei Fragen zu finden: Welche Knoten sind verbunden und wie findet man einen kürzesten Pfad? Um diese Probleme zu lösen, stellt die Netzwerkanalyse-Bibliothek den Dijkstra-Algorithmus zur Verfügung.

Dijkstras Algorithmus findet den kürzesten Pfad von einem der Knoten des Graphen zu allen anderen und die Werte der Optimierungsparameter. Die Ergebnisse können als ein Baum für die kürzesten Pfade dargestellt werden.

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

• Genau ein Knoten besitzt keine ankommende Kante — Die Wurzel des Baumes

• Alle anderen Knoten haben genau eine ankommende Kante

• Wenn Knoten B von Knoten A aus erreichbar ist, dann ist der Pfad von A nach B der einzige verfügbare und optimale (kürzeste) auf diesem Graph

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

Warnung

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. Kürzeste Pfade ermitteln

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

Nun verfügen wir über den Pfad in Form der umgekehrten Liste von Konten (Knoten sind in umgekehrter Reihenfolge vom Endknoten zum Startknoten aufgelistet), die beim Durchlaufen des Pfades besucht werden.

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. Erreichbarkeitsgebiete

Das Erreichbarkeitsgebiet für den Knoten A ist die Teilmenge der Graphknoten, die von dem Knoten A aus erreichbar sind, und die Kosten der Pfade von A zu diesen Knoten einen bestimmten Wert nicht überschreiten.

Deutlicher wird das am folgenden Beispiel: „Es gibt eine Feuerwache. Welche Teile der Stadt kann ein Feuerwehrauto in 5 Minuten erreichen? 10 Minuten? 15 Minuten?“. Die Antworten auf diese Fragen sind die Erreichbarkeitsgebiete der Feuerwache.

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.

Ein schwierigeres Problem besteht darin, die Begrenzung des Erreichbarkeitsgebietes zu ermitteln. Die untere Begrenzung ist die Gruppe von Knoten, die noch erreichbar sind, und die obere Begrenzung die Gruppe von Knoten, welche nicht mehr erreichbar sind. Tatsächlich ist dies einfach: Es ist die Grenze der Erreichbarkeit basierend auf den Kanten des kürzesten Pfadbaums, für den der Quellknoten der Kante erreichbar ist, und der Zielknoten der Kante nicht.

Hier ein Beispiel:

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