Importante

La traducción es un esfuerzo comunitario al que puedes unirte. Esta página está actualmente traducida en un 62.50%.

19. Biblioteca de análisis de redes

Consejo

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

La biblioteca de análisis de red se puede utilizar para:

  • crear un gráfico matemático a partir de datos geográficos (capas vectoriales de polilínea)

  • implementar métodos básicos de la teoría de grafos (actualmente solo el algoritmo de Dijkstra)

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)

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

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.

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.

La conversión de una capa vectorial al gráfico se realiza mediante el patrón de programación Constructor. Un gráfico se construye utilizando el llamado Director. Solo hay un Director por ahora: QgsVectorLayerDirector. El director establece la configuración básica que se utilizará para construir un gráfico a partir de una capa vectorial de línea, utilizada por el constructor para crear el gráfico. Actualmente, como en el caso del director, solo existe un constructor: QgsGraphBuilder, que crea objetos QgsGraph. Es posible que desee implementar sus propios constructores que crearán un gráfico compatible con bibliotecas como BGL o 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 tiempo de sumergirse en el proceso.

  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. Es necesario entonces crear una estrategia para calcular propiedades de borde

    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. Y decirle al director sobre esta estrategia

    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. Análisis gráfico

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.

El árbol de ruta más corto es un gráfico ponderado dirigido (o más precisamente un árbol) con las siguientes propiedades:

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

  • all other vertices have only one incoming edge

  • 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

Para obtener la ruta del árbol más corta, use los métodos shortestTree() y dijkstra() de la clase QgsGraphAnalyzer . Es recomendable usar el método dijkstra() porque funciona más rápido y usa la memoria de manera más eficiente.

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 — gráfica entrante

  • startVertexIdx — índice del punto en el árbol (la raíz del árbol)

  • criterionNum — número de propiedad de borde a usar (comenzando desde 0).

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

El método dijkstra() tiene los mismos argumentos, pero devuelve una tupla de matrices:

  • En la primera matriz, el elemento n contiene el índice del borde entrante o -1 si no hay bordes entrantes.

  • En la segunda matriz, el elemento n contiene la distancia desde la raíz del árbol hasta el vértice n o DOUBLE_MAX si el vértice n es inaccesible desde la raíz.

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

Advertencia

Use este código solo como ejemplo, crea muchos objetos QgsRubberBand y puede ser lento en conjuntos de datos extensos.

 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. Encontrar la ruta más corta

Para encontrar la ruta óptima entre dos puntos, se utiliza el siguiente enfoque. Ambos puntos (inicio A y final B) están «vinculados» al gráfico cuando se construye. Luego, usando el método shortestTree () o dijkstra() construimos el árbol de ruta más corto con la raiz en el punto inicial A. En el mismo árbol también encontramos el punto final B y comenzamos a caminar a través del árbol desde el punto B al punto A. Todo el algoritmo se puede escribir como:

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. Áreas de disponibilidad

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.

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.

Para encontrar las áreas de disponibilidad podemos usar el método dijkstra() de la clase QgsGraphAnalyzer. Es suficiente comparar los elementos de la matriz de costos con un valor predefinido. Si el costo [i] es menor o igual a un valor predefinido, entonces el vértice i está dentro del área de disponibilidad, de lo contrario, está fuera.

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