Los fragmentos de código en esta página necesitan las siguientes adiciones si está fuera de la consola de pyqgis:

from qgis.core import (
  QgsVectorLayer,
  QgsPointXY,
)

18. Biblioteca de análisis de redes

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)

La librería de análisis de redes fue creada por funciones básicas de exportación del complemento núcleo RoadGraph y ahora se puede utilizar los metodos en complementos o directamente de la consola Python.

18.1. Información general

Brevemente, un caso de uso típico se puede describir como:

  1. Crear gráfica de geodatos (normalmente de capa vectorial de polilíneas)

  2. ejecutar análisis gráfico

  3. utilizar resultados de análisis (por ejemplo, visualizarlos)

18.2. Contruir un gráfico

Lo primero que hay que hacer — es preparar la entrada de datos, que es convertir una capa vectorial en un gráfico. Todas las acciones adicionales utilizarán esta gráfica, no la capa.

Como fuente podemos utilizar una capa vectorial de polilínea. Los nodos de las polilíneas se convierten en vértices del gráfico, y los segmentos de la polilínea son bordes de gráfico. Si varios nodos tienen la misma coordenada entonces ellos tienen el mimso vértice gráfico. Por lo que dos líneas que tienen un nodo en común se conectaran entre si.

Además durante la creación del gráfico se puede «arreglar» («atar») a la capa vectorial de entrada cualquier número de puntos adicionales. Para cada punto adicional se encontrará una coincidencia — el vértice gráfico más cercano o el borde gráfico más cercano. En el último caso el borde será dividido y un nuevo vértice se añadirá.

Los atributos de la capa vectorial y la longitud de un borde se puede utilizar como las propiedades de un borde.

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 graphs 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 tiempo de sumergirse en el proceso.

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

from qgis.analysis import *

Después algunos ejemplos para crear un director

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
# 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)

Para construir un director debemos pasar a una capa vectorial, que se utilizará como fuente para la estructura gráfica y la información sobre el movimiento permitido en cada segmento de carretera (movimiento unidireccional o bidireccional, dirección directa o inversa). La llamada se parece a esto

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

Y aquí esta la lista completa de lo que significan estos parámetros:

  • vectorLayer — vector layer used to build the graph

  • directionFieldId — índice de la tabla de atributos de campo, donde se almacena información acerca de dirección de carreteras. Si -1, entonces no utilice esta información en absoluto. Un entero.

  • directDirectionValue — el valor del campo de carreteras con dirección directa (mover desde la primer punto de línea a la última). Un texto.

  • reverseDirectionValue — valor del campo de carreteras con dirección inversa (mover del último punto de línea al primero). Un texto.

  • bothDirectionValue — valor de campo para carreteras bidireccionales (para cada carretera podemos mover del primer punto al último y del último al primero). Un texto.

  • 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

Es necesario entonces crear una estrategia para calcular propiedades de borde

1
2
3
4
5
6
7
# 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)

Y decirle al director sobre esta estrategia

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

También podemos definir varios puntos, que se utilizarán en el análisis. Por ejemplo

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

Ahora todo está en su lugar para que podamos construir el gráfico y «atar» a estos puntos

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

Construir el grafo puede tomar tiempo (que depende del número de elementos y tamaño de una capa). tiedPoints es una lista con coordenadas de puntos «tied». Cuando la operación de construcción se finalizo podemos obtener la gráfica y utilizarlo para el análisis

graph = builder.graph()

Con el siguiente código podemos obtener el índice del vértice de nuestros puntos

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

18.3. Análisis gráfico

El análisis de redes es utilizado para encontrar respuestas a dos preguntas: que vértices estan conectados y cómo encontrar la ruta más corta. Para resolver estos problemas la librería de análisis de redes proporciona el algoritmo Dijkstra.

El algoritmo Dijkstra encuentra la ruta más corta de uno de los vértices del grafo a todos los otros y los valores de los parámetros de optimización, El resultado puede ser representado como un árbol de la ruta más corta.

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

  • sólo un vértice no tiene bordes entrantes — la raíz del árbol

  • todos los otros vértices sólo tienen un borde entrante

  • Si el vértice B es accesible desde el vértice A, entonces el camino de A a B es la única ruta disponible y es optima (más corta) en este grafo

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

Advertencia

Use this code only as an example, it creates a lot of QgsRubberBand objects and may be slow on large datasets.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
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

Lo mismo pero usando el método: meth: dijkstra <qgis.analysis.QgsGraphAnalyzer.dijkstra>

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
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())

18.3.1. Encontrar la ruta más corta

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:

1
2
3
4
5
6
7
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

En este punto tenemos la ruta, en el formulario de la lista invertida de vértices (los vértices están listados en orden invertida del punto final al punto inicial) que serán visitados durante el viaje por este camino.

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
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)

18.3.2. Áreas de disponibilidad

El área de la disponibilidad para el vértice A es el subconjunto de vértices del grafo que son accesibles desde el vértice A y el costo de los caminos de la A a estos vértices son no es mayor que cierto valor.

Más claramente esto se puede demostrar con el siguiente ejemplo: «Hay una estación de bomberos ¿Qué partes de la ciudad puede un camión de bomberos alcanzar en 5 minutos? 10 minutos? 15 minutos?». Las respuestas a estas preguntas son las zonas de la estación de bomberos de la disponibilidad.

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 problema más difícil es conseguir los límites de la zona de disponibilidad. El borde inferior es el conjunto de vértices que son todavía accesibles, y el borde superior es el conjunto de vértices que no son accesibles. De hecho esto es simple: es la frontera disponibilidad basado en los bordes del árbol de ruta más corta para los que el vértice origen del contorno es más accesible y el vértice destino del borde no lo es.

Aquí tiene un ejemplo

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
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(), True)
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(), True)
    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))