Bibliothèque d’analyse de réseau

Avertissement

Despite our constant efforts, information beyond this line may not be updated for QGIS 3. Refer to https://qgis.org/pyqgis/master for the python API documentation or, give a hand to update the chapters you know about. Thanks.

Depuis la révision ee19294562 (QGIS >= 1.8), la nouvelle bibliothèque d’analyse de réseau a été ajoutée à la bibliothèque principale d’analyse de QGIS. La bibliothèque :

  • créé un graphe mathématique à partir de données géographiques (couches vecteurs de polylignes)

  • implémente des méthodes simples de la théorie des graphes (pour l’instant, uniquement avec l’algorithme Dijkstra).

La bibliothèque d’analyse de réseau a été créée en exportant les fonctions de l’extension principale RoadGraph. Vous pouvez en utiliser les méthodes dans des extensions ou directement dans la console Python.

Information générale

Voici un résumé d’un cas d’utilisation typique:

  1. créer un graphe depuis les données géographiques (en utilisant une couche vecteur de polylignes)

  2. lancer une analyse de graphe

  3. utiliser les résultats d’analyse (pour les visualiser par exemple)

Construire un graphe

La première chose à faire est de préparer les données d’entrée, c’est à dire de convertir une couche vecteur en graphe. Les actions suivantes utiliseront ce graphe et non la couche.

Comme source de données, on peut utiliser n’importe quelle couche vecteur de polylignes. Les nœuds des polylignes deviendront les sommets du graphe et les segments des polylignes seront les arcs du graphes. Si plusieurs nœuds ont les mêmes coordonnées alors ils composent le même sommet de graphe. Ainsi, deux lignes qui ont en commun un même nœud sont connectées ensemble.

Pendant la création d’un graphe, il est possible de « forcer » (« lier ») l’ajout d’un ou de plusieurs points additionnels à la couche vecteur d’entrée. Pour chaque point additionnel, un lien sera créé: le sommet du graphe le plus proche ou l’arc de graphe le plus proche. Dans le cas final, l’arc sera séparé en deux et un nouveau sommet sera ajouté.

Les attributs de la couche vecteur et la longueur d’un segment peuvent être utilisés comme propriétés du segment.

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: QgsLineVectorLayerDirector. 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 graphs compatible with such libraries as BGL or NetworkX.

To calculate edge properties the programming pattern strategy is used. For now only QgsDistanceArcProperter strategy is available, that takes into account the length of the route. 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.

Il est temps de plonger dans le processus.

D’abord, nous devrions importer le module networkanalysis pour utiliser la bibliothèque

from qgis.networkanalysis import *

Ensuite, quelques exemples pour créer un directeur

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

Pour construire un directeur, il faut lui fournir une couche vecteur qui sera utilisée comme source pour la structure du graphe ainsi que des informations sur les mouvements permis sur chaque segment de route (sens unique ou déplacement bidirectionnel, direct ou inversé). L’appel au directeur se fait de la manière suivante

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

Voici la liste complète de la signification de ces paramètres:

  • vl — couche vecteur utilisée pour construire le graphe

  • directionFieldId — index du champ de la table d’attribut où est stockée l’information sur la direction de la route. Si -1 est utilisé, cette information n’est pas utilisée. Un entier.

  • directDirectionValue — valeur du champ utilisé pour les routes avec une direction directe (déplacement du premier point de la ligne au dernier). Une chaîne de caractères.

  • reverseDirectionValue — valeur du champ utilisé pour les routes avec une direction inverse (déplacement du dernier point de la ligne au premier). Une chaîne de caractères.

  • bothDirectionValue — valeur du champ utilisé pour les routes bidirectionelles (pour ces routes, on peut se déplacer du premier point au dernier et du dernier au premier). Une chaîne de caractères.

  • defaultDirection — direction par défaut de la route. Cette valeur sera utilisée pour les routes où le champ directionFieldId` n'est pas paramétré ou qui a une valeur différente des trois valeurs précédentes. Un entier ``1 indique une direction directe, 2 indique une direction inverse et 3 indique les deux directions.

Il est ensuite impératif de créer une stratégie de calcul des propriétés des arcs:

properter = QgsDistanceArcProperter()

Et d’informer le directeur à propos de cette stratégie:

director.addProperter(properter)

Now we can use the builder, which will create the graph. The QgsGraphBuilder class constructor takes several arguments:

  • crs — système de coordonnées de référence à utiliser. Argument obligatoire.

  • otfEnabled — utiliser ou non la projection « à la volée ». La valeur par défaut est const:True (oui, utiliser OTF).

  • topologyTolerance — la tolérance topologique. La valeur par défaut est 0.

  • ellipsoidID — ellipsoïde à utiliser. Par défaut « WGS84 ».

# only CRS is set, all other values are defaults
builder = QgsGraphBuilder(myCRS)

Nous pouvons également définir plusieurs points qui seront utilisés dans l’analyse, par exemple:

startPoint = QgsPoint(82.7112, 55.1672)
endPoint = QgsPoint(83.1879, 54.7079)

Maintenant que tout est en place, nous pouvons construire le graphe et lier ces points dessus:

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

La construction du graphe peut prendre du temps (qui dépend du nombre d’entités dans la couche et de la taille de la couche). tiedPoints est une liste qui contient les coordonnées des points liés. Lorsque l’opération de construction est terminée, nous pouvons récupérer le graphe et l’utiliser pour l’analyse:

graph = builder.graph()

Avec le code qui suit, nous pouvons récupérer les index des arcs de nos points:

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

Analyse de graphe

L’analyse de graphe est utilisée pour trouver des réponses aux deux questions: quels arcs sont connectés et comment trouver le plus court chemin ? Pour résoudre ces problèmes la bibliothèque d’analyse de graphe fournit l’algorithme de Dijkstra.

L’algorithme de Dijkstra trouve le plus court chemin entre un des arcs du graphe par rapport à tous les autres en tenant compte des paramètres d’optimisation. Ces résultats peuvent être représentés comme un arbre du chemin le plus court.

L’arbre du plus court chemin est un graphe pondéré de direction (plus précisément un arbre) qui dispose des propriétés suivantes:

  • Seul un arc n’a pas d’arcs entrants: la racine de l’arbre.

  • Tous les autres arcs n’ont qu’un seul arc entrant.

  • Si un arc B est atteignable depuis l’arc A alors le chemin de A vers B est le seul chemin disponible et il est le chemin optimal (le plus court) sur ce graphe.

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 — graphe en entrée

  • startVertexIdx — index du point sur l’arbre (la racine de l’arbre)

  • criterionNum — nombre de propriétés d’arc à utiliser (en partant de 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)

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

Avertissement

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.networkanalysis import *
from qgis.PyQt.QtCore import *
from qgis.PyQt.QtGui 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 the dijkstra method

from qgis.core import *
from qgis.gui import *
from qgis.networkanalysis import *
from qgis.PyQt.QtCore import *
from qgis.PyQt.QtGui 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())

Trouver les chemins les plus courts

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

A ce niveau, nous avons le chemin, sous la forme d’une liste inversée d’arcs (les arcs sont listés dans un ordre inversé, depuis le point de la fin vers le point de démarrage) qui seront traversés lors de l’évolution sur le chemin.

Here is the sample code for QGIS Python Console (you will need to select 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.networkanalysis import *
from qgis.PyQt.QtCore import *
from qgis.PyQt.QtGui 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 the dijkstra method

from qgis.core import *
from qgis.gui import *
from qgis.networkanalysis import *
from qgis.PyQt.QtCore import *
from qgis.PyQt.QtGui 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)

Surfaces de disponibilité

La surface de disponibilité d’un arc A est le sous-ensemble des arcs du graphe qui sont accessibles à partir de l’arc A et où le coût des chemins à partir de A vers ces arcs ne dépasse pas une certaine valeur.

Plus clairement, cela peut être illustré par l’exemple suivant: « Il y a une caserne de pompiers. Quelles parties de la ville peuvent être atteintes par un camion de pompier en 5 minutes ? 10 minutes ? 15 minutes ? » La réponse à ces questions correspond aux surface de disponibilité de la caserne de pompiers.

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.

Un problème plus difficile à régler est d’obtenir les frontières de la surface de disponibilité. La frontière inférieure est constituée par l’ensemble des arcs qui sont toujours accessibles et la frontière supérieure est composée des arcs qui ne sont pas accessibles. En fait, c’est très simple: c’est la limite de disponibilité des arcs de l’arbre du plus court chemin pour lesquels l’arc source de l’arc est accessible et l’arc cible ne l’est pas.

Voici un exemple:

from qgis.core import *
from qgis.gui import *
from qgis.networkanalysis import *
from qgis.PyQt.QtCore import *
from qgis.PyQt.QtGui 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))