19. Bibliotheek Netwerkanalyse

Hint

De codesnippers op deze pagina hebben de volgende import nodig als u buiten de console van PyQGIS bent:

from qgis.core import (
  QgsVectorLayer,
  QgsPointXY,
)

De bibliotheek Netwerkanalyse kan worden gebruikt voor:

• rekenkundige grafieken maken uit geografische gegevens (polylijn vectorlagen)

• basismethoden implementeren vanuit grafiektheorie (momenteel alleen Dijkstra’s algoritme)

De bibliotheek Network analysis werd gemaakt door het exporteren van basisfuncties vanuit de bronplug-in RoadGraph en nu kunt u de methoden daarvan in plug-ins gebruiken of direct vanuit de console voor Python.

19.1. Algemene informatie

In het kort kan een typisch gebruik worden omschreven als:

1. maakt rekenkundige grafiek uit geo-gegevens (gewoonlijk polylijn vectorlaag)

2. voert grafiekanalyse uit

3. gebruikt resultaten van analyse (door ze, bijvoorbeeld, te visualiseren)

19.2. Een grafiek bouwen

Het eerste dat u moet doen — is om invoergegevens voor te bereiden, dat is een vectorlaag converteren naar een grafiek. Alle verdere acties zullen deze grafiek gebruiken, niet de laag.

Als een bron kunnen we elke polylijn vectorlaag gebruiken. Knopen van de polylijnen worden punten in de grafiek, en segmenten van de polylijnen worden randen van de grafiek. Indien verscheidene knopen dezelfde coördinaten hebben dan zijn zij dezelfde knop in de grafiek. Dus twee lijnen die een gemeenschappelijk knoop hebben worden aan elkaar verbonden.

Aanvullend, gedurende het maken van de grafiek is het mogelijk om de een willekeurige aantal aanvullende punten “vast te zetten” (“te verbinden”) aan de invoer vectorlaag. Voor elk aanvullend punt zal een overeenkomst worden gevonden — het dichtstbijzijnde punt in de grafiek of de dichtstbijzijnde gelegen rand. In het laatste geval zal de rand worden gesplitst en een nieuw punt worden toegevoegd.

Attributen van de vectorlaag en de lengte van een rand kunnen worden gebruikt als de eigenschappen van een rand.

Converteren van een vectorlaag naar de grafiek wordt gedaan met behulp van het Builder programmeringspatroon. Een grafiek wordt geconstrueerd met behulp van een zogenaamde Director. Er is nu nog slechts één Director: QgsVectorLayerDirector. De director stelt de basisinstellingen in die zullen worden gebruikt om een grafiek uit een lijn-vectorlaag te maken, gebruikt door de builder om de grafiek te maken. Momenteel, net zoals in het geval van de Director, bestaat er slechts één builder: QgsGraphBuilder, die objecten QgsGraph maakt. U wilt misschien uw eigen builders implementeren die een grafiek bouwen die compatibel is met bibliotheken zoals BGL of NetworkX.

Voor het berekenen van de eigenschappen van de rand wordt het programmeringspatroon strategy gebruikt. Momenteel zijn alleen beschikbaar QgsNetworkDistanceStrategy strategie (die rekening houdt met de lengte van de route) en QgsNetworkSpeedStrategy (die ook rekening houdt met de snelheid). U kunt uw eigen strategie implementeren die alle noodzakelijke parameters zal gebruiken. De plug-in RoadGraph gebruikt bijvoorbeeld een strategie die de reistijd berekent met behulp van de lengte van de rand en een waarde voor de snelheid uit attributen.

Het is tijd om het proces in te duiken.

Als eerste, om deze bibliotheek te kunnen gebruiken, zouden we de module analysis moeten importeren

from qgis.analysis import *

Dan enkele voorbeelden voor het maken van een 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)

We zouden, om een director te construeren, een vectorlaag door moeten geven, die zal worden gebruikt als de bron voor de structuur van de grafiek en informatie over toegestane bewegingen over elke segment van de weg (één richting of beide, directe of tegengestelde richting). De aanroep ziet er uit zoals deze

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

En hier is een volledige lijst van wat deze parameters betekenen:

• vectorLayer — vectorlaag gebruikt om de grafiek te bouwen

• directionFieldId — index van het attribuut tabelveld, waar informatie over de richting van de wegen is opgeslagen. Indien -1, gebruik deze informatie dan helemaal niet. Een integer.

• directDirectionValue — veldwaarde voor wegen met een directe richting (verplaatsen vanaf het eerste punt op de lijn tot het laatste). Een string.

• reverseDirectionValue — veldwaarde voor wegen met een tegengestelde richting (verplaatsen vanaf het laatste punt op de lijn tot het eerste). Een string.

• bothDirectionValue — veldwaarde voor wegen in beide richtingen (voor dergelijke wegen kunnen we verplaatsen van het eerste punt naar het laatste en van laatste naar eerste). Een string.

• defaultDirection — standaard richting van de weg. Deze waarde zal worden gebruikt voor die wegen waar het veld directionFieldId niet is ingesteld of een andere waarde heeft dan een van de hierboven gespecificeerde drie waarden. Mogelijke waarden zijn:

  • QgsVectorLayerDirector.DirectionForward — Éénrichting vooruit

  • QgsVectorLayerDirector.DirectionBackward — Éénrichting tegengesteld

  • QgsVectorLayerDirector.DirectionBoth — Twee richtingen

Het is dan nodig om een strategie te maken voor het berekenen van de eigenschappen van de rand

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

En de director vertellen over deze strategie

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

Nu kunnen we de builder gebruiken, wat de grafiek zal maken. De klasseconstructor QgsGraphBuilder kan verschillende argumenten aannemen:

• crs — te gebruiken coördinaten referentiesysteem. Verplicht argument.

• otfEnabled — opnieuw projecteren met “Gelijktijdige CRS-transformatie gebruiken” gebruiken of niet. Standaard True (OTF gebruiken).

• topologyTolerance — topologische tolerantie. Standaardwaarde is 0.

• ellipsoidID — te gebruiken ellipsoïde. Standaard “WGS84”.

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

Ook kunnen we verscheidene punten definiëren, die zullen worden gebruikt in de analyse. Bijvoorbeeld

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

Nu is alles op zijn plaats dus kunnen we de grafiek bouwen en deze punten daaraan “verbinden”

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

Bouwen van de grafiek kan enige tijd vergen (wat afhankelijk is van het aantal objecten in een laag en de grootte van de laag). tiedPoints is een lijst met coördinaten van de “verbonden” punten. Als de bewerking van het bouwen is voltooid kunnen we de grafiek nemen en die gebruiken voor de analyse

graph = builder.graph()

Met de volgende code kunnen we de vertex-indexen verkrijgen van onze punten

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

19.3. Grafiekanalyse

Netwerkanalyse wordt gebruikt om antwoord te vinden op twee vragen: welke punten zijn verbonden en hoe het kortste pad te vinden. De bibliotheek Network analysis verschaft Dijkstra’s algoritme om deze problemen op te lossen.

Dijkstra’s algoritme zoekt de kortste route van één van de punten van de grafiek naar alle andere en de waarden van de parameters voor optimalisatie. De resultaten kunnen worden weergegeven als een kortste pad-boom.

De kortste pad-boom is een gedirigeerde gewogen grafiek (of meer precies een boom) met de volgende eigenschappen:

• slechts één punt heeft geen inkomende randen — de wortel van de boom

• alle andere punten hebben slechts één inkomende rand

• als punt B bereikbaar is vanuit punt A, dan is het pad van A naar B het enige beschikbare pad en is het optimaal (kortste) op deze grafiek

Gebruik, om de boom van het kortste pad te verkrijgen, de methoden shortestTree() en dijkstra() van de klasse QgsGraphAnalyzer. Aanbevolen wordt om de methode dijkstra() te gebruiken omdat die sneller werkt en het geheugen meer efficiënt gebruikt.

De methode shortestTree() is handig wanneer u over de boom van het kortste pad wilt wandelen. Het maakt altijd een nieuw grafiekobject (QgsGraph) en accepteert drie variabelen:

• source — grafiek voor invoer

• startVertexIdx — index van het punt op de boom (de wortel van de boom)

• criterionNum — nummer van te gebruiken eigenschap van de rand (beginnend vanaf 0).

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

De methode dijkstra() heeft dezelfde argumenten, maar geeft twee arrays terug. In het eerste array bevat element n de index van de inkomende rand of -1 als er geen inkomende randen zijn. In de tweede array bevat element n de afstand van de wortel van de boom tot het punt n of DOUBLE_MAX als het punt n onbereikbaar is vanuit de wortel.

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

Hier is enige eenvoudige code om de boom van het kortste pad weer te geven met de grafiek die is gemaakt met de methode shortestTree() (selecteer de lijnenlaag in het paneel Lagen en vervang de coördinaten door die van uzelf).

Waarschuwing

Gebruik deze code alleen als voorbeeld, Het maakt heel veel objecten QgsRubberBand en zou zeer traag kunnen zijn voor grote gegevenssets.

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

Hetzelfde, maar met de methode 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. Kortste pad zoeken

De volgende benadering wordt gebruikt om het optimale pad tussen twee punten te zoeken. Beide punten (begin A en einde B) zijn “verbonden” met de grafiek wanneer die wordt gebouwd. Dan bouwen we met de methode shortestTree() of dijkstra() de boom voor het kortste pad met de wortel in beginpunt A. In dezelfde boom zoeken we ook naar eindpunt B en beginnen te lopen door de boom vanaf punt B naar punt A. Het gehele algoritme kan worden geschreven als:

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

Op dit punt hebben we het pad, in de vorm van de geïnverteerde lijst van punten (punten zijn vermeld in de omgekeerde volgorde van eindpunt naar beginpunt) die zullen worden bezocht gedurende het lopen over dit pad.

Hier is de voorbeeldcode voor de console van Python in QGIS (u zou misschien een lijnenlaag willen laden en selecteren in de inhoudsopgave en de coördinaten in de code te vervangen door die van uzelf) dat de methode shortestTree() gebruikt

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)

En hier is hetzelfde voorbeeld, maar voor de methode 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. Beschikbare gebieden

Het beschikbare gebied voor punt A is de subset van punten op de grafiek die toegankelijk zijn vanuit punt A en de kosten van de paden van A naar deze punten zijn niet groter dan een bepaalde waarde.

Dit kan duidelijker worden weergegeven met behulp van het volgende voorbeeld: “Er is een brandweergarage. Welke delen van de stad kan een brandweerauto bereiken in 5 minuten? 10 minuten? 15 minuten?”. De antwoorden op deze vragen zijn de beschikbare gebieden voor deze brandweergarage.

We kunnen de methode dijkstra() van de klasse QgsGraphAnalyzer gebruiken om de beschikbare gebieden te zoeken. Het is voldoende om de elementen van de array met kosten te vergelijken met een vooraf gedefinieerde waarde. Als de kosten[i] minder zijn dan of gelijk zijn aan een vooraf gedefinieerde waarde, dan ligt punt i binnen het beschikbare gebied, anders ligt het er buiten.

Een wat moeilijker probleem is om de grenzen van de beschikbare gebieden te verkrijgen. De ondergrens is de set punten die nog steeds toegankelijk zijn, en de bovengrens is de set punten die niet toegankelijk zijn. In feite is dit eenvoudig: het is de grens van beschikbaarheid, gebaseerd op de randen van de boom van het kortste pad waarvoor het bronpunt van de rand toegankelijk is en het doelpunt van de rand is dat niet.

Hier is een voorbeeld

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