Wichtig

Übersetzen ist eine Gemeinschaftsleistung Sie können mitmachen. Diese Seite ist aktuell zu 35.71% übersetzt.

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)

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

As a source we can use any polyline vector layer. Nodes of the polylines become graph vertices, and segments of the polylines are graph edges. If several nodes have the same coordinates then they are the same graph vertex. So two lines that have a common node become connected to each other.

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 available. You can implement your own strategy that will use all necessary parameters.

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

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

    from qgis.analysis import *
    
  2. Then some examples for creating a director:

    1# Don't use information about road direction from layer attributes,
    2# all roads are treated as two-way
    3director = QgsVectorLayerDirector(
    4    vectorLayer, -1, "", "", "", QgsVectorLayerDirector.DirectionBoth
    5)
    
    1# Use field with index 5 as source of information about road direction.
    2# one-way roads with direct direction have attribute value "yes",
    3# one-way roads with reverse direction have the value "1", and accordingly
    4# bidirectional roads have "no". By default roads are treated as two-way.
    5# This scheme can be used with OpenStreetMap data
    6director = QgsVectorLayerDirector(
    7    vectorLayer, 5, "yes", "1", "no", QgsVectorLayerDirector.DirectionBoth
    8)
    

    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 (find more details on the parameters at qgis.analysis.QgsVectorLayerDirector):

    1director = QgsVectorLayerDirector(
    2    vectorLayer,
    3    directionFieldId,
    4    directDirectionValue,
    5    reverseDirectionValue,
    6    bothDirectionValue,
    7    defaultDirection,
    8)
    
  3. Nun muss eine Strategie für die Ermittlung der Kanteneigenschaften festgelegt werden

    1# The index of the field that contains information about the edge speed
    2attributeId = 1
    3# Default speed value
    4defaultValue = 50
    5# Conversion from speed to metric units ('1' means no conversion)
    6toMetricFactor = 1
    7strategy = QgsNetworkSpeedStrategy(attributeId, defaultValue, toMetricFactor)
    
  4. und der Director über diese Strategie informiert werden

    director = QgsVectorLayerDirector(vectorLayer, -1, "", "", "", 3)
    director.addStrategy(strategy)
    
  5. Now we can use the builder, which will create the graph, using the QgsGraphBuilder class constructor.

    # only CRS is set, all other values are defaults
    builder = QgsGraphBuilder(vectorLayer.crs())
    
  6. Also we can define several points, which will be used in the analysis. For example:

    startPoint = QgsPointXY(1179720.1871, 5419067.3507)
    endPoint = QgsPointXY(1180616.0205, 5419745.7839)
    
  7. Now all is in place so we can build the graph and „tie“ these points to it:

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

    Building the graph can take some time (which depends on the number of features in a layer and layer size). tiedPoints is a list with coordinates of „tied“ points.

  8. When the build operation is finished we can get the graph and use it for the analysis:

    graph = builder.graph()
    
  9. With the next code we can get the vertex indexes of our points:

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

19.2. Graphenanalyse

Networks analysis is used to find answers to two questions: which vertices are connected and how to find a shortest path. To solve these problems the network analysis library provides Dijkstra’s algorithm.

Dijkstra’s algorithm finds the shortest route from one of the vertices of the graph to all the others and the values of the optimization parameters. The results can be represented as a shortest path tree.

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

  • all other vertices have only one incoming edge

  • 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 a tuple of 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() or the dijkstra() 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.

 1from qgis.core import *
 2from qgis.gui import *
 3from qgis.analysis import *
 4from qgis.PyQt.QtCore import *
 5from qgis.PyQt.QtGui import *
 6
 7vectorLayer = QgsVectorLayer(
 8    "testdata/network.gpkg|layername=network_lines", "lines"
 9)
10director = QgsVectorLayerDirector(
11    vectorLayer, -1, "", "", "", QgsVectorLayerDirector.DirectionBoth
12)
13strategy = QgsNetworkDistanceStrategy()
14director.addStrategy(strategy)
15builder = QgsGraphBuilder(vectorLayer.crs())
16
17pStart = QgsPointXY(1179661.925139, 5419188.074362)
18tiedPoint = director.makeGraph(builder, [pStart])
19pStart = tiedPoint[0]
20
21graph = builder.graph()
22
23idStart = graph.findVertex(pStart)
24
25tree = QgsGraphAnalyzer.shortestTree(graph, idStart, 0)
26
27i = 0
28while i < tree.edgeCount():
29    rb = QgsRubberBand(iface.mapCanvas())
30    rb.setColor(Qt.red)
31    rb.addPoint(tree.vertex(tree.edge(i).fromVertex()).point())
32    rb.addPoint(tree.vertex(tree.edge(i).toVertex()).point())
33    i = i + 1

19.2.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:

1assign T = B
2while T != B
3    add point T to path
4    get incoming edge for point T
5    look for point TT, that is start point of this edge
6    assign T = TT
7add point A to path

At this point we have the path, in the form of the inverted list of vertices (vertices are listed in reversed order from end point to start point) that will be visited during traveling by this path.

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() or dijkstra() method:

 1from qgis.core import *
 2from qgis.gui import *
 3from qgis.analysis import *
 4
 5from qgis.PyQt.QtCore import *
 6from qgis.PyQt.QtGui import *
 7
 8vectorLayer = QgsVectorLayer(
 9    "testdata/network.gpkg|layername=network_lines", "lines"
10)
11director = QgsVectorLayerDirector(
12    vectorLayer, -1, "", "", "", QgsVectorLayerDirector.DirectionBoth
13)
14strategy = QgsNetworkDistanceStrategy()
15director.addStrategy(strategy)
16
17builder = QgsGraphBuilder(vectorLayer.sourceCrs())
18
19startPoint = QgsPointXY(1179661.925139, 5419188.074362)
20endPoint = QgsPointXY(1180942.970617, 5420040.097560)
21
22tiedPoints = director.makeGraph(builder, [startPoint, endPoint])
23tStart, tStop = tiedPoints
24
25graph = builder.graph()
26idxStart = graph.findVertex(tStart)
27
28tree = QgsGraphAnalyzer.shortestTree(graph, idxStart, 0)
29
30idxStart = tree.findVertex(tStart)
31idxEnd = tree.findVertex(tStop)
32
33if idxEnd == -1:
34    raise Exception("No route!")
35
36# Add last point
37route = [tree.vertex(idxEnd).point()]
38
39# Iterate the graph
40while idxEnd != idxStart:
41    edgeIds = tree.vertex(idxEnd).incomingEdges()
42    if len(edgeIds) == 0:
43        break
44    edge = tree.edge(edgeIds[0])
45    route.insert(0, tree.vertex(edge.fromVertex()).point())
46    idxEnd = edge.fromVertex()
47
48# Display
49rb = QgsRubberBand(iface.mapCanvas())
50rb.setColor(Qt.green)
51
52# This may require coordinate transformation if project's CRS
53# is different from layer's CRS
54for p in route:
55    rb.addPoint(p)

19.2.2. Erreichbarkeitsgebiete

The area of availability for vertex A is the subset of graph vertices that are accessible from vertex A and the cost of the paths from A to these vertices are not greater that some value.

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.

A more difficult problem is to get the borders of the area of availability. The bottom border is the set of vertices that are still accessible, and the top border is the set of vertices that are not accessible. In fact this is simple: it is the availability border based on the edges of the shortest path tree for which the source vertex of the edge is accessible and the target vertex of the edge is not.

Here is an example:

 1director = QgsVectorLayerDirector(
 2  vectorLayer, -1, "", "", "", QgsVectorLayerDirector.DirectionBoth
 3)
 4strategy = QgsNetworkDistanceStrategy()
 5director.addStrategy(strategy)
 6builder = QgsGraphBuilder(vectorLayer.crs())
 7
 8
 9pStart = QgsPointXY(1179661.925139, 5419188.074362)
10delta = iface.mapCanvas().getCoordinateTransform().mapUnitsPerPixel() * 1
11
12rb = QgsRubberBand(iface.mapCanvas())
13rb.setColor(Qt.green)
14rb.addPoint(QgsPointXY(pStart.x() - delta, pStart.y() - delta))
15rb.addPoint(QgsPointXY(pStart.x() + delta, pStart.y() - delta))
16rb.addPoint(QgsPointXY(pStart.x() + delta, pStart.y() + delta))
17rb.addPoint(QgsPointXY(pStart.x() - delta, pStart.y() + delta))
18
19tiedPoints = director.makeGraph(builder, [pStart])
20graph = builder.graph()
21tStart = tiedPoints[0]
22
23idStart = graph.findVertex(tStart)
24
25(tree, cost) = QgsGraphAnalyzer.dijkstra(graph, idStart, 0)
26
27upperBound = []
28r = 1500.0
29i = 0
30tree.reverse()
31
32while i < len(cost):
33    if cost[i] > r and tree[i] != -1:
34        outVertexId = graph.edge(tree [i]).toVertex()
35        if cost[outVertexId] < r:
36            upperBound.append(i)
37    i = i + 1
38
39for i in upperBound:
40    centerPoint = graph.vertex(i).point()
41    rb = QgsRubberBand(iface.mapCanvas())
42    rb.setColor(Qt.red)
43    rb.addPoint(QgsPointXY(centerPoint.x() - delta, centerPoint.y() - delta))
44    rb.addPoint(QgsPointXY(centerPoint.x() + delta, centerPoint.y() - delta))
45    rb.addPoint(QgsPointXY(centerPoint.x() + delta, centerPoint.y() + delta))
46    rb.addPoint(QgsPointXY(centerPoint.x() - delta, centerPoint.y() + delta))