19. Bibliothèque d’analyse de réseau

Indication

Les extraits de code sur cette page nécessitent les importations suivantes si vous êtes en dehors de la console pyqgis :

from qgis.core import (
  QgsVectorLayer,
  QgsPointXY,
)

La bibliothèque d’analyse de réseau peut être utilisée pour :

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

  • mettre en œuvre les méthodes de base de la théorie des graphes (actuellement, seul l’algorithme de Dijkstra est utilisé)

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.

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

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

La conversion d’une couche vectorielle vers le graphe se fait en utilisant le modèle de programmation Builder. Un graphe est construit en utilisant ce qu’on appelle un Director. Il n’y a qu’un seul Director pour le moment : QgsVectorLayerDirector. Le Director définit les paramètres de base qui seront utilisés pour construire un graphique à partir d’une couche de vecteurs linéaires, utilisée par le builder pour créer le graphique. Actuellement, comme dans le cas du director, un seul constructeur existe : QgsGraphBuilder, qui crée des objets QgsGraph. Vous pouvez vouloir implémenter vos propres constructeurs qui construiront un graphe compatible avec des bibliothèques telles que BGL ou NetworkX.

Pour calculer les propriétés des bords, on utilise le modèle de programmation strategy. Pour l’instant, seules sont disponibles QgsNetworkDistanceStrategy stratégie (qui prend en compte la longueur du trajet) et QgsNetworkSpeedStrategy (qui prend également en compte la vitesse). Vous pouvez mettre en œuvre votre propre stratégie qui utilisera tous les paramètres nécessaires. Par exemple, le plugin RoadGraph utilise une stratégie qui calcule le temps de parcours en utilisant la longueur des bords et la valeur de la vitesse à partir des attributs.

Il est temps de plonger dans le processus.

Tout d’abord, pour utiliser cette bibliothèque, nous devons importer le module d’analyse

from qgis.analysis import *

Ensuite, quelques exemples pour créer un directeur

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

Pour construire un director, nous devons passer une couche vecteur qui servira de source pour la structure du graphique et les informations sur les mouvements autorisés sur chaque segment de route (mouvement unidirectionnel ou bidirectionnel, direction directe ou inverse). L’appel ressemble à ceci

1director = QgsVectorLayerDirector(vectorLayer,
2                                  directionFieldId,
3                                  directDirectionValue,
4                                  reverseDirectionValue,
5                                  bothDirectionValue,
6                                  defaultDirection)

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

  • vectorLayer — 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 de la route par défaut. Cette valeur sera utilisée pour les routes où le champ directionFieldId n’est pas défini ou a une valeur différente de l’une des trois valeurs spécifiées ci-dessus. Les valeurs possibles sont :

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

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)

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

director = QgsVectorLayerDirector(vectorLayer, -1, '', '', '', 3)
director.addStrategy(strategy)

Nous pouvons maintenant utiliser le constructeur, qui va créer le graphique. Le constructeur de la classe QgsGraphBuilder prend plusieurs arguments :

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

  • otfEnabled — utiliser la reprojection « à la volée » ou non. Par défaut True (utilisez OTF).

  • TopologyTolerance — 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(vectorLayer.crs())

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

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

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

19.3. 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 chemin le plus court est un graphique pondéré dirigé (ou plus précisément un arbre) ayant les 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.

Pour obtenir l’arbre du chemin le plus court, utilisez les méthodes shortestTree() et dijkstra() de la classe QgsGraphAnalyzer. Il est recommandé d’utiliser la méthode dijkstra() car elle fonctionne plus rapidement et utilise la mémoire de manière plus efficace.

La méthode shortestTree() est utile lorsque vous voulez vous promener dans l’arbre qui a le chemin le plus court. Elle crée toujours un nouvel objet graphique (QgsGraph) et accepte trois variables :

  • source — Graphique d’entrée

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

  • criterionNum — nombre de propriétés de bord à utiliser (à partir de 0).

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

La méthode dijkstra() a les mêmes arguments, mais renvoie deux tableaux. Dans le premier tableau, l’élément n contient l’index du front entrant ou -1 s’il n’y a pas de bord entrant. Dans le second tableau, l’élément n contient la distance entre la racine de l’arbre et le sommet n ou DOUBLE_MAX si le sommet n est inaccessible depuis la racine.

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

Voici un code très simple pour afficher l’arbre du chemin le plus court en utilisant le graphe créé avec la méthode shortestTree() (sélectionner la couche de ligne dans le panneau Couches et remplacer les coordonnées par les vôtres).

Avertissement

Utilisez ce code uniquement à titre d’exemple, il crée beaucoup d’objets QgsRubberBand et peut être lent sur de grands ensembles de données.

 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('testdata/network.gpkg|layername=network_lines', 'lines')
 8director = QgsVectorLayerDirector(vectorLayer, -1, '', '', '', QgsVectorLayerDirector.DirectionBoth)
 9strategy = QgsNetworkDistanceStrategy()
10director.addStrategy(strategy)
11builder = QgsGraphBuilder(vectorLayer.crs())
12
13pStart = QgsPointXY(1179661.925139,5419188.074362)
14tiedPoint = director.makeGraph(builder, [pStart])
15pStart = tiedPoint[0]
16
17graph = builder.graph()
18
19idStart = graph.findVertex(pStart)
20
21tree = QgsGraphAnalyzer.shortestTree(graph, idStart, 0)
22
23i = 0
24while (i < tree.edgeCount()):
25  rb = QgsRubberBand(iface.mapCanvas())
26  rb.setColor (Qt.red)
27  rb.addPoint (tree.vertex(tree.edge(i).fromVertex()).point())
28  rb.addPoint (tree.vertex(tree.edge(i).toVertex()).point())
29  i = i + 1

Même chose mais en utilisant la méthode dijkstra().

 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('testdata/network.gpkg|layername=network_lines', 'lines')
 8
 9director = QgsVectorLayerDirector(vectorLayer, -1, '', '', '', QgsVectorLayerDirector.DirectionBoth)
10strategy = QgsNetworkDistanceStrategy()
11director.addStrategy(strategy)
12builder = QgsGraphBuilder(vectorLayer.crs())
13
14pStart = QgsPointXY(1179661.925139,5419188.074362)
15tiedPoint = director.makeGraph(builder, [pStart])
16pStart = tiedPoint[0]
17
18graph = builder.graph()
19
20idStart = graph.findVertex(pStart)
21
22(tree, costs) = QgsGraphAnalyzer.dijkstra(graph, idStart, 0)
23
24for edgeId in tree:
25  if edgeId == -1:
26    continue
27  rb = QgsRubberBand(iface.mapCanvas())
28  rb.setColor (Qt.red)
29  rb.addPoint (graph.vertex(graph.edge(edgeId).fromVertex()).point())
30  rb.addPoint (graph.vertex(graph.edge(edgeId).toVertex()).point())

19.3.1. Trouver les chemins les plus courts

Pour trouver le chemin optimal entre deux points, l’approche suivante est utilisée. Les deux points (début A et fin B) sont « liés » au graphique lors de sa construction. Ensuite, en utilisant la méthode shortestTree() ou dijkstra() nous construisons l’arbre du chemin le plus court avec une racine au point de départ A. Dans le même arbre, nous trouvons également le point de fin B et commençons à parcourir l’arbre du point B au point A. L’algorithme complet peut être écrit sous la forme :

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

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.

Voici un exemple de code pour la console Python QGIS (vous devrez peut-être charger et sélectionner une couche de lignes dans la table des matières et remplacer les coordonnées dans le code par les vôtres) qui utilise la méthode shortestTree().

 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('testdata/network.gpkg|layername=network_lines', 'lines')
 9builder = QgsGraphBuilder(vectorLayer.sourceCrs())
10director = QgsVectorLayerDirector(vectorLayer, -1, '', '', '', QgsVectorLayerDirector.DirectionBoth)
11
12startPoint = QgsPointXY(1179661.925139,5419188.074362)
13endPoint = QgsPointXY(1180942.970617,5420040.097560)
14
15tiedPoints = director.makeGraph(builder, [startPoint, endPoint])
16tStart, tStop = tiedPoints
17
18graph = builder.graph()
19idxStart = graph.findVertex(tStart)
20
21tree = QgsGraphAnalyzer.shortestTree(graph, idxStart, 0)
22
23idxStart = tree.findVertex(tStart)
24idxEnd = tree.findVertex(tStop)
25
26if idxEnd == -1:
27    raise Exception('No route!')
28
29# Add last point
30route = [tree.vertex(idxEnd).point()]
31
32# Iterate the graph
33while idxEnd != idxStart:
34    edgeIds = tree.vertex(idxEnd).incomingEdges()
35    if len(edgeIds) == 0:
36        break
37    edge = tree.edge(edgeIds[0])
38    route.insert(0, tree.vertex(edge.fromVertex()).point())
39    idxEnd = edge.fromVertex()
40
41# Display
42rb = QgsRubberBand(iface.mapCanvas())
43rb.setColor(Qt.green)
44
45# This may require coordinate transformation if project's CRS
46# is different than layer's CRS
47for p in route:
48    rb.addPoint(p)

Et voici le même échantillon mais en utilisant la méthode dijkstra().

 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('testdata/network.gpkg|layername=network_lines', 'lines')
 9director = QgsVectorLayerDirector(vectorLayer, -1, '', '', '', QgsVectorLayerDirector.DirectionBoth)
10strategy = QgsNetworkDistanceStrategy()
11director.addStrategy(strategy)
12
13builder = QgsGraphBuilder(vectorLayer.sourceCrs())
14
15startPoint = QgsPointXY(1179661.925139,5419188.074362)
16endPoint = QgsPointXY(1180942.970617,5420040.097560)
17
18tiedPoints = director.makeGraph(builder, [startPoint, endPoint])
19tStart, tStop = tiedPoints
20
21graph = builder.graph()
22idxStart = graph.findVertex(tStart)
23idxEnd = graph.findVertex(tStop)
24
25(tree, costs) = QgsGraphAnalyzer.dijkstra(graph, idxStart, 0)
26
27if tree[idxEnd] == -1:
28    raise Exception('No route!')
29
30# Total cost
31cost = costs[idxEnd]
32
33# Add last point
34route = [graph.vertex(idxEnd).point()]
35
36# Iterate the graph
37while idxEnd != idxStart:
38    idxEnd = graph.edge(tree[idxEnd]).fromVertex()
39    route.insert(0, graph.vertex(idxEnd).point())
40
41# Display
42rb = QgsRubberBand(iface.mapCanvas())
43rb.setColor(Qt.red)
44
45# This may require coordinate transformation if project's CRS
46# is different than layer's CRS
47for p in route:
48    rb.addPoint(p)

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

Pour trouver les zones de disponibilité, nous pouvons utiliser la méthode dijkstra() de la classe QgsGraphAnalyzer. Il suffit de comparer les éléments du tableau des coûts avec une valeur prédéfinie. Si le coût [i] est inférieur ou égal à une valeur prédéfinie, alors le sommet i est à l’intérieur de la zone de disponibilité, sinon il est à l’extérieur.

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:

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