5

I'm a new user of QGIS (and have experience using ArcGIS but zero programming skills) and 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?

Hoping not to have to do each step manually, one at a time...

Black squares - addresses, white circles - subway stops, gray lines - roads

3

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 networkanalysis libray. (That is the same library used by the Road graph plugin.)

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.