Finding the shortest route using Dijkstra's algorithm

2864
5
12-23-2014 04:09 PM
XanderBakker
Esri Esteemed Contributor
5 5 2,864

Last month a user mentioned the Dijkstra algorithm to determine the shortest path between two nodes in a network. My first reaction was to use Network Analyst instead, but to honor the essence of a previous useless post: Japanese visual multiplication with lines using arcpy and a dash of patriotism (Dijkstra was Dutch) made me look into this a little more.

I came across this post; Dijkstra's algorithm for shortest paths « Python recipes « ActiveState Code  and from this post I actually used the code from this one; Rebrained! » Blog Archive » Shortest Path – Python . Remember, basic skills for developing code are search, copy and paste

It is basically recursive code that will ultimately get you the route between two nodes on the network. What you need is a "graph". This is a nested dictionary that lists all the connections between the from and to nodes and the corresponding distances. To be able to use the code I need to convert my "network" into this graph format.

What you need is a network (a simple line featureclass), in which all intersections of lines consist of nodes. You will also need a point feature class containing the start and end point. These points should correspond with nodes of the network.This is the result:

DijkstraRoute.png

I came up with the following code:

#-------------------------------------------------------------------------------
# Name:        dijkstra.py
# Purpose:     determine route
#
# Author:      xbakker
#
# Created:     02/12/2014
#-------------------------------------------------------------------------------

import arcpy
import sys

def main():
    # input featureclasses
    fc_red = r"D:\Xander\GeoNet\Dijkstra\fgdb\test.gdb\red"
    fc_stops = r"D:\Xander\GeoNet\Dijkstra\fgdb\test.gdb\stops2"

    # create line dictionary
    dct_lines = createLineDict(fc_red)

    # create node and inverted node dictionary
    dct_nodes = createNodes(dct_lines)

    # update line dictionary with nodes
    dct_lines_nodes = combineLinesAndNodes(dct_lines, dct_nodes)

    # create graph
    dct_graph = createGraph(dct_nodes, dct_lines_nodes)

    # determine the node of start and end of route
    node_start, node_end = getStartEndNodes(fc_stops, dct_nodes)

    # determine route
    result = shortestpath(dct_graph, node_start, node_end)

    # create line with the result
    dist = result[0]
    route = result[1]
    print "Resulting route    : {0}".format(route)
    print "Length of the route: {0}\n".format(dist)
    cnt = 0
    for node_id in route:
        cnt += 1
        print "{0}\t{1}\t{2}".format(cnt, dct_nodes[node_id][0], dct_nodes[node_id][1])


def createLineDict(fc):
    """Crear dict de lineas con datos relevantes"""
    #OBJECTID (line): length, (Xstart, Ystart), (Xend, Yend)
    dct_lines = {}
    flds = ("OID@", "SHAPE@", "SHAPE@LENGTH")
    dct_lines = {r[0]: [r[2], (r[1].firstPoint.X, r[1].firstPoint.Y),
                (r[1].lastPoint.X, r[1].lastPoint.Y)]
                for r in arcpy.da.SearchCursor(fc, flds)}
    return dct_lines

def createNodes(dct):
    """Crear dict de nodos"""
    node_format = "node_{0}"
    cnt = 0
    dct_nodes = {}
    dct_nodes_inv = {}

    for line in dct.values():
        start = line[1]
        if not start in dct_nodes.values():
            cnt += 1
            node_id = node_format.format(cnt)
            dct_nodes[node_id] = start

        end = line[2]
        if not end in dct_nodes.values():
            cnt += 1
            node_id = node_format.format(cnt)
            dct_nodes[node_id] = end
    return dct_nodes

def combineLinesAndNodes(dct_lines, dct_nodes):
    """Crear dict con lineas y codigo de nodos"""
    dct_lines_nodes = {}
    for oid, line in dct_lines.items():
        start = getNodeID(line[1], dct_nodes)
        end = getNodeID(line[2], dct_nodes)
        dct_lines_nodes[oid] = [line[0], start, end]
    return dct_lines_nodes

def getNodeID(pnt, dct_nodes):
    for node_id, coords in dct_nodes.items():
        if coords == pnt:
            return node_id
            break

def createGraph(dct_nodes, dct_lines_nodes):
    res = {}
    for node_id in dct_nodes.keys():
        links = {}
        nodes_start = [[line[2], line[0]] for line in dct_lines_nodes.values() if line[1] == node_id]
        for node in nodes_start:
            links[node[0]] = node[1]
        nodes_end = [[line[1], line[0]] for line in dct_lines_nodes.values() if line[2] == node_id]
        for node in nodes_end:
            links[node[0]] = node[1]
        res[node_id] = links
    return res

def getStartEndNodes(fc, dct_nodes):
    cnt = 0
    with arcpy.da.SearchCursor(fc, ("SHAPE@")) as curs:
        for row in curs:
            cnt += 1
            if cnt == 1:
                pnt = row[0].firstPoint
                node_start = (pnt.X, pnt.Y)
                node_start = getNodeID(node_start, dct_nodes)
            elif cnt == 2:
                pnt = row[0].lastPoint
                node_end = (pnt.X, pnt.Y)
                node_end = getNodeID(node_end, dct_nodes)
            else:
                break
    return node_start, node_end

def shortestpath(graph,start,end,visited=[],distances={},predecessors={}):
    """Find the shortest path between start and end nodes in a graph"""
    # we've found our end node, now find the path to it, and return
    if start==end:
        path=[]
        while end != None:
            path.append(end)
            end=predecessors.get(end,None)
        return distances[start], path[::-1]
    # detect if it's the first time through, set current distance to zero
    if not visited:
        distances[start]=0
    # process neighbors as per algorithm, keep track of predecessors
    for neighbor in graph[start]:
        if neighbor not in visited:
            neighbordist = distances.get(neighbor,sys.maxint)
            tentativedist = distances[start] + graph[start][neighbor]
            if tentativedist < neighbordist:
                distances[neighbor] = tentativedist
                predecessors[neighbor]=start
    # neighbors processed, now mark the current node as visited
    visited.append(start)
    # finds the closest unvisited node to the start
    unvisiteds = dict((k, distances.get(k,sys.maxint)) for k in graph if k not in visited)
    closestnode = min(unvisiteds, key=unvisiteds.get)
    # now we can take the closest node and recurse, making it current
    return shortestpath(graph,closestnode,end,visited,distances,predecessors)

if __name__ == '__main__':
    main()

Have fun!

5 Comments
About the Author
Solution Engineer for the Utilities Sector @ Esri Colombia - Ecuador - Panamá sr GIS Advisor / Python - Arcpy developer / GIS analyst / technical project leader / lecturer and GeoNet moderator, focusing on innovations in the field of GIS. Specialties: ArcGIS, Python, ArcGIS Enterprise, ArcGIS Online, Arcade, Configurable Apps, WAB, Mobile Apps, Insights, Spatial Analysis, LiDAR / 3D Laser Scanning / Point Clouds. UNME http://nl.linkedin.com/in/xanderbakker/ http://www.slideshare.net/XanderBakker http://www.scribd.com/xbakker http://twitter.com/#!/XanderBakker
Labels