19. Libreria analisi di rete

Suggerimento

I frammenti di codice di questa pagina necessitano delle seguenti importazioni se sei è al di fuori della console pyqgis:

from qgis.core import (
  QgsVectorLayer,
  QgsPointXY,
)

La libreria di analisi di rete può essere utilizzata per:

• creare un grafico matematico a partire da dati geografici (layer vettoriali polilineari)

• implementare i metodi di base della teoria dei grafi (attualmente solo l’algoritmo di Dijkstra)

La libreria di analisi di rete è stata creata esportando le funzioni di base dal plugin principale di RoadGraph e ora puoi utilizzare i suoi metodi nei plugin o direttamente dalla console Python.

19.1. Informazioni generali

In breve, un caso d’uso tipico può essere descritto come segue:

1. creare un grafo da geodati (solitamente layer vettoriali polilineari)

2. eseguire l’analisi del grafo

3. utilizzare i risultati dell’analisi (ad esempio, visualizzarli)

19.2. Costruire un grafo

La prima cosa da fare è preparare i dati di input, cioè convertire un layer vettoriale in un grafo. Tutte le azioni successive utilizzeranno questo grafo, non il layer.

Come sorgente possiamo utilizzare qualsiasi layer vettoriale di polilinee. I nodi delle polilinee diventano vertici del grafo e i segmenti delle polilinee sono bordi del grafo. Se più nodi hanno le stesse coordinate, allora sono lo stesso vertice del grafo. Quindi due linee che hanno un nodo in comune diventano collegate tra loro.

Inoltre, durante la creazione del grafo è possibile «fissare» («legare») al layer del vettore di input un numero qualsiasi di punti aggiuntivi. Per ogni punto aggiuntivo verrà trovata una corrispondenza — con il vertice del grafo più vicino o con il bordo del grafo più vicino. In quest’ultimo caso, il bordo verrà diviso e verrà aggiunto un nuovo vertice.

Gli attributi del layer vettoriale e la lunghezza di un bordo possono essere utilizzati come proprietà di un bordo.

La conversione da un layer vettoriale a grafo viene effettuata utilizzando il modello di programmazione Builder. Un grafo viene costruito usando un cosiddetto Director. Per ora esiste un solo Director: QgsVectorLayerDirector. Il Director imposta le impostazioni di base che verranno utilizzate per costruire un grafico da un layer vettoriale di linee, utilizzato dal costruttore per creare il grafico. Attualmente, come nel caso del director, esiste un solo costruttore: QgsGraphBuilder, che crea oggetti QgsGraph. Puoi implementare i propri costruttori per creare un grafo compatibile con librerie come BGL o NetworkX.

Per calcolare le proprietà dei bordi si utilizza lo schema di programmazione strategy. Per ora sono disponibili solo le strategie QgsNetworkDistanceStrategy (che tiene conto della lunghezza del percorso) e QgsNetworkSpeedStrategy (che considera anche la velocità). Puoi implementare la tua strategia personale, che utilizzerà tutti i parametri necessari. Ad esempio, il plugin RoadGraph utilizza una strategia che calcola il tempo di percorrenza utilizzando la lunghezza dei bordi e il valore della velocità dagli attributi.

È il momento di addentrarsi nel processo.

Prima di tutto, per utilizzare questa libreria dobbiamo importare il modulo di analisi

from qgis.analysis import *

Poi alcuni esempi di creazione di 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)

Per costruire un director, occorre passare un layer vettoriale che sarà usato come fonte per la struttura del grafo e le informazioni sul movimento consentito su ogni segmento stradale (movimento unidirezionale o bidirezionale, direzione diretta o inversa). La call si presenta come segue

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

Ecco l’elenco completo del significato di questi parametri:

• vectorLayer — layer vettoriale usato per costruire il grafo

• directionFieldId — indice del campo della tabella degli attributi, dove sono memorizzate le informazioni sulla direzione delle strade. Se -1, non utilizzare affatto questa informazione. Un numero intero.

• directDirectionValue — Valore del campo per le strade con direzione diretta (che si muovono dal primo punto della linea all’ultimo). Una stringa.

• reverseDirectionValue — valore del campo per le strade con direzione inversa (che si muovono dall’ultimo punto della linea al primo). Una stringa.

• bothDirectionValue — valore del campo per le strade bidirezionali (per tali strade ci si può muovere dal primo punto all’ultimo e dall’ultimo al primo). Una stringa.

• defaultDirection — direzione predefinita della strada. Questo valore sarà utilizzato per le strade in cui il campo directionFieldId non è impostato o ha un valore diverso da uno dei tre valori specificati sopra. I valori possibili sono:

  • QgsVectorLayerDirector.DirectionForward — Diretto a senso unico

  • QgsVectorLayerDirector.DirectionBackward — Inversione unidirezionale

  • QgsVectorLayerDirector.DirectionBoth — bidirezionale

È quindi necessario creare una strategia per il calcolo delle proprietà dei bordi

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

E informare il director di questa strategia

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

Ora possiamo usare il costruttore, che creerà il grafo. Il costruttore della classe QgsGraphBuilder utilizza diversi parametri:

• crs — sistema di riferimento delle coordinate da usare. Parametro obbligatorio.

• otfEnabled — utilizzare la riproiezione «al volo» o no. Per impostazione predefinita True (usa OTF).

• topologyTolerance — tolleranza topologica. Il valore predefinito è 0.

• ellipsoid — ellissoide da usare. Per impostazione predefinita «WGS84».

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

Inoltre, possiamo definire diversi punti che verranno utilizzati nell’analisi. Ad esempio

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

Ora tutto è a posto e possiamo costruire il grafo e «legare» questi punti ad esso

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

La costruzione del grafo può richiedere un certo tempo (che dipende dal numero di elementi in un layer e dalla dimensione del layer). tiedPoints è un elenco con le coordinate dei punti «legati». Quando l’operazione di costruzione è terminata, si può ottenere il grafico e utilizzarlo per l’analisi

graph = builder.graph()

Con il codice che segue possiamo ottenere gli indici dei vertici dei nostri punti

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

19.3. Analisi grafo

L’analisi di rete viene utilizzata per trovare le risposte a due domande: quali vertici sono connessi e come trovare il percorso più breve. Per risolvere questi problemi, la libreria di analisi di rete fornisce l’algoritmo di Dijkstra.

L’algoritmo di Dijkstra trova il percorso più breve da un vertice del grafo a tutti gli altri e i valori dei parametri di ottimizzazione. I risultati possono essere rappresentati come un albero del percorso più breve.

L’albero del percorso più breve è un grafo diretto pesato (o più precisamente un albero) con le seguenti proprietà:

• solo un vertice non ha bordi in entrata: la radice dell’albero.

• tutti gli altri vertici hanno un solo bordo in entrata

• se il vertice B è raggiungibile dal vertice A, allora il percorso da A a B è il solo percorso disponibile ed è ottimale (il più breve) su questo grafo

Per ottenere l’albero del percorso più breve, utilizzare i metodi shortestTree() e dijkstra() della classe QgsGraphAnalyzer. Si consiglia di utilizzare il metodo dijkstra() perché funziona più velocemente e utilizza la memoria in modo più efficiente.

Il metodo shortestTree() è utile quando si vuole percorrere l’albero del percorso più breve. Crea sempre un nuovo oggetto grafo (QgsGraph) e accetta tre variabili:

• source — grafo in ingresso

• startVertexIdx — indice del punto sull” albero (la radice dell’albero)

• criterionNum — numero di proprietà del bordo da usare (a partire da 0).

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

Il metodo dijkstra() ha gli stessi argomenti, ma restituisce due array. Il primo elemento n contiene l’indice del bordo entrante o -1 se non ci sono bordi entranti. Nel secondo elemento n contiene la distanza dalla radice dell’albero al vertice n o DOUBLE_MAX se il vertice n non è raggiungibile dalla radice.

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

Ecco un codice molto semplice per visualizzare l’albero del percorso più breve utilizzando il grafo creato con il metodo shortestTree() (seleziona il layer linestring nel pannello Layer e sostituisci le coordinate con le tue).

Avvertimento

Usa questo codice solo come esempio; crea molti oggetti QgsRubberBand e può essere lento su grandi insiemi di dati.

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

Stessa cosa, ma utilizzando il metodo 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. Trovare i percorsi più brevi

Per trovare il percorso ottimale tra due punti si utilizza il seguente approccio. Entrambi i punti (inizio A e fine B) sono «legati» al grafo quando viene costruito. Quindi, utilizzando il metodo shortestTree() o dijkstra() si costruisce l’albero del percorso più breve con radice nel punto iniziale A. Nello stesso albero si trova anche il punto finale B e si inizia a percorrere l’albero dal punto B al punto A. L’intero algoritmo può essere scritto come:

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 questo punto abbiamo il percorso, sotto forma di elenco invertito di vertici (i vertici sono elencati in ordine inverso, dal punto finale al punto iniziale) che saranno visitati durante il percorso.

Ecco il codice di esempio per QGIS Python Console (potrebbe essere necessario caricare e selezionare un layer linestring in legenda e sostituire le coordinate nel codice con le tue) che utilizza il metodo 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)

Ed ecco lo stesso esempio, ma utilizzando il metodo 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. Zone di disponibilità

La zona di disponibilità per il vertice A è il sottoinsieme dei vertici del grafo che sono accessibili dal vertice A e il costo dei percorsi da A a questi vertici non è superiore a un certo valore.

Questo può essere mostrato più chiaramente con il seguente esempio: «C’è una stazione dei vigili del fuoco. Quali zone della città può raggiungere un camion dei pompieri in 5 minuti? 10 minuti? 15 minuti?». Le risposte a queste domande sono le zone di disponibilità della caserma dei pompieri.

Per trovare le zone di disponibilità possiamo usare il metodo dijkstra() della classe QgsGraphAnalyzer. È sufficiente confrontare gli elementi dell’array cost con un valore predefinito. Se cost[i] è minore o uguale a un valore predefinito, allora il vertice i è all’interno della zona di disponibilità, altrimenti è fuori.

Un problema più difficile è quello di ottenere i confini della zona di disponibilità. Il confine inferiore è l’insieme dei vertici ancora accessibili, mentre il confine superiore è l’insieme dei vertici non accessibili. In realtà è semplice: si tratta del confine di disponibilità basato sui bordi dell’albero del cammino più breve per i quali il vertice di origine del bordo è accessibile e il vertice di destinazione del bordo no.

Ecco un esempio

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