Shortest path from points in one layer to points in another layer via roads in QGIS

Question

I'm looking for a way to calculate the shortest path from addresses to subway stations via a network of roads.

I have the Road Graph plugin, however, the addresses and subway stations are point layers that do not overlap with vertices in the roads layer.

In QGIS, is there a way to calculate the shortest distance from the address/subway station points to the roads, create vertices at those locations, then calculate the shortest distance via the roads network?

• black squares - addresses
• white circles - subway stops
• gray lines - roads

Answer 1

You'll need some Python scripting because there is currently no GUI solution for this issue.

To get you started, have a look at my point layer to route script for Processing which is based on the QGIS network analysis library. (That is the same library used by the Road graph plugin.)

<pre>
##Routing tools=group
##points=vector point
##network=vector line
##route=output vector line

from PyQt4.QtCore import *
from PyQt4.QtGui import *

from qgis.core import *
from qgis.gui import *
from qgis.networkanalysis import *

from processing.tools.vector import VectorWriter

point_layer = processing.getObject(points)
network_layer = processing.getObject(network)
writer = VectorWriter(route, None, [QgsField("order", QVariant.Int)], network_layer.dataProvider().geometryType(), network_layer.crs() )

# prepare graph
vl = network_layer
director = QgsLineVectorLayerDirector( vl, -1, '', '', '', 3 )
properter = QgsDistanceArcProperter()
director.addProperter( properter )
crs = vl.crs()
builder = QgsGraphBuilder( crs )

# prepare points
features = processing.features(point_layer)
point_count = point_layer.featureCount()

points = []

for f in features:
  points.append(f.geometry().asPoint())

tiedPoints = director.makeGraph( builder, points )
graph = builder.graph()

route_vertices = []

for i in range(0,point_count-1):
    progress.setPercentage(int(100 * i/ point_count))

    from_point = tiedPoints[i]
    to_point = tiedPoints[i+1]

    from_id = graph.findVertex(from_point)
    to_id = graph.findVertex(to_point)

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

    if tree[to_id] == -1:
        continue # ignore this point pair
    else:
        # collect all the vertices between the points
        route_points = []
        curPos = to_id 
        while (curPos != from_id):
            route_points.append( graph.vertex( graph.arc( tree[ curPos ] ).inVertex() ).point() )
            curPos = graph.arc( tree[ curPos ] ).outVertex()

        route_points.append(from_point)

    # add a feature
    fet = QgsFeature()
    fet.setGeometry(QgsGeometry.fromPolyline(route_points))
    fet.setAttributes([i])
    writer.addFeature(fet)

del writer
</pre>

Answer 2

I had a very similar problem, finding the amenities (restaurants etc) within a certain distance by road from access points along a bike path. I solved it with the toolbox script shown below, which is a modified version of @underdark's script, that I found through a similar question here.

The script asks for two-point vector layers, a network line vector layer and a maximum acceptable route distance, and finds the shortest path between each pair of access/amenity points. Then, for each path shorter than the given distance, the script writes the path, id's for the associated access and amenity points, and the route distance, to the output line vector. (If the given distance is zero or less, all paths are output regardless of length.)

You'll need to modify your point vector's attributes, for both access (subway entrance) and amenity layers (addresses), by adding a column called "acam_id" containing a unique identifying number for each point.

#Definition of inputs and outputs
#==================================
##Routing tools=group
##access=vector point
##amenities=vector point
##network=vector line
##distance=number 1600
##Routes=output vector line

#Algorithm body
#==================================
from PyQt4.QtCore import *
from PyQt4.QtGui import *

from qgis.core import *
from qgis.gui import *
from qgis.networkanalysis import *

from processing.tools.vector import VectorWriter
import processing


progress.setText("Building Network...")

access_layer = processing.getObject(access)
amenities_layer = processing.getObject(amenities)
network_layer = processing.getObject(network)
writer = VectorWriter(Routes, None, [
    QgsField("acc_id", QVariant.Int), QgsField("amen_id", QVariant.Int),   QgsField("length", QVariant.Double)], \
    network_layer.dataProvider().geometryType(), network_layer.crs() )

# prepare graph
vl = network_layer
director = QgsLineVectorLayerDirector( vl, -1, '', '', '', 3 )
properter = QgsDistanceArcProperter()
director.addProperter( properter )
crs = vl.crs()
builder = QgsGraphBuilder( crs )

# prepare points
access_features = processing.features(access_layer)
access_count = access_layer.featureCount()
amenities_features = processing.features(amenities_layer)
amenities_count = amenities_layer.featureCount()
point_count = access_count + amenities_count

points = []
ids = []

for f in access_features:
    points.append(f.geometry().asPoint())
    ids.append(f['acam_id'])
for f in amenities_features:
    points.append(f.geometry().asPoint())
    ids.append(f['acam_id'])


tiedPoints = director.makeGraph( builder, points )
graph = builder.graph()

route_vertices = []
progress.setText("Processing Routes...")
for i in range(0,access_count):
    progress.setPercentage(int(100 * i/ access_count))
    from_point = tiedPoints[i]
    from_id = graph.findVertex(from_point)
    (tree,cost) = QgsGraphAnalyzer.dijkstra(graph,from_id,0)

    for j in range(access_count,point_count):
        to_point = tiedPoints[j]
        to_id = graph.findVertex(to_point)

        if tree[to_id] != -1 and (cost[to_id] <= distance or distance <= 0):
            # collect all the vertices between the points
            route_points = []
            curPos = to_id 
            while (curPos != from_id):
                route_points.append( graph.vertex( graph.arc( tree[ curPos ] ).inVertex() ).point() )
                curPos = graph.arc( tree[ curPos ] ).outVertex()

            route_points.append(from_point)

            # add a feature
            fet = QgsFeature()
            fet.setGeometry(QgsGeometry.fromPolyline(route_points))
            fet.setAttributes([ids[i], ids[j], cost[to_id]])
            writer.addFeature(fet)

del writer