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)

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.

19.1. Información general 

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

Crear gráfica de geodatos (normalmente de capa vectorial de polilíneas)
ejecutar análisis gráfico
utilizar resultados de análisis (por ejemplo, visualizarlos)

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

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.

Para calcular las propiedades del borde, se usó el patrón de programación estrategia. Por ahora solo estrategia QgsNetworkDistanceStrategy (que tiene en cuenta la longitud de la ruta) y QgsNetworkSpeedStrategy (que también considera la velocidad) están disponibles. Puede implementar su propia estrategia que utilizará todos los parámetros necesarios. Por ejemplo, el complemento RoadGraph utiliza una estrategia que calcula el tiempo de viaje utilizando la longitud del borde y el valor de velocidad de los atributos.

Es tiempo de sumergirse en el proceso.

En primer lugar, para utilizar esta biblioteca debemos importar el módulo de análisis

from qgis.analysis import *

Después algunos ejemplos para crear un director

# 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 una capa vectorial que se utilizará como fuente para la estructura del gráfico y la información sobre el movimiento permitido en cada segmento de la carretera (movimiento unidireccional o bidireccional, dirección directa o inversa). La llamada se ve así

director = QgsVectorLayerDirector(vectorLayer,
                                  directionFieldId,
                                  directDirectionValue,
                                  reverseDirectionValue,
                                  bothDirectionValue,
                                  defaultDirection)

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

vectorLayer — capa vectorial usada para construir la gráfica
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 — dirección de la carretera predeterminada. Este valor se usará para aquellas carreteras donde el campo directionFieldId no está configurado o tiene algún valor diferente de cualquiera de los tres valores especificados anteriormente. Los posibles valores son:
- QgsVectorLayerDirector.DirectionForward — dirección de un sentido
- QgsVectorLayerDirector.DirectionBackward — Sentido único reversible
- QgsVectorLayerDirector.DirectionBoth — Bi-direccional

Es necesario entonces crear una estrategia para calcular propiedades de borde

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

Ahora podemos usar el constructor, que creará el gráfico. La clase QgsGraphBuilder constructor toma varios argumentos:

crs — sistema de referencia de coordenadas a utilizar. Argumento obligatorio.
otfEnabled — usar reproyección «al vuelo» o no. Por defecto True (usa OTF).
topologyTolerance — tolerancia topológica. Por defecto es 0.
ellipsoidID — elipsoide a usar. Por defecto «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])

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

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

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.

El método shortestTree() es útil cuando desea caminar alrededor del árbol del camino más corto. Siempre crea un nuevo objeto gráfico (QgsGraph) y acepta tres 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 dos matrices. En el primer elemento de matriz, n contiene el índice del borde entrante o -1 si no hay bordes entrantes. En el segundo elemento de matriz, n contiene la distancia desde la raíz del árbol al vértice` n` o DOUBLE_MAX si el vértice n es inalcanzable desde la raíz.

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

Aquí hay un código muy simple para mostrar el árbol de ruta más corto usando el gráfico creado con el método shortestTree() (seleccione la capa de cadena de líneas en el panel Capas y reemplace las coordenadas por las suyas).

Advertencia

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

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

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

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

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.

Aquí está el código de muestra para Consola de Pyhton QGIS (es posible que deba cargar y seleccionar una capa de cadena de líneas en TOC y reemplazar las coordenadas en el código por las suyas) que usa el método:meth:shortestTree() <qgis.analysis.QgsGraphAnalyzer.shortestTree>

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)

Y aquí está la misma muestra pero usando el método dijkstra()

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)

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

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.

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

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

19. Biblioteca de análisis de redes

19.1. Información general

19.2. Contruir un gráfico

19.3. Análisis gráfico

19.3.1. Encontrar la ruta más corta

19.3.2. Áreas de disponibilidad