As I currently only need all shortest paths between two set of points, I ended up by calculating the shortest paths on my own with the Python binding of the igraph library. What I have is an input data set with points and a dataset with the street network as set of polylines with start node, end node and length. What I do is first find for each data point the nearest point on the road network. Then I create the graph from the road network. For each data point, I add two new edges to the graph. This is actually the biggest part of the code. Maybe there are more elegant and shorter versions. I know for each data point the id of the polyline nearest to it. For each such polyline I iterate over the segments and check if the nearest point of the data point is on the segment or the points between.
'''
Created on 13.05.2013
Create graph from road network dataset and find shortest path
'''
import arcpy
import igraph
from collections import defaultdict
poi_layer = 1 # Layer with points of interests.
street_layer = 2 # Layer with streets as polylines. Start and end vertex ids in fields
# Find for every poi closest point on road entwork
arcpy.analysis.Near(poi_layer, street_layer, location="LOCATION")
# Set of ids of all polylines containing a nearest point
intersect_points = defaultdict(set)
with arcpy.da.SearchCursor(poi_layer, ("SHAPE@", "NEAR_FID", "NEAR_DIST", "NEAR_X", "NEAR_Y", "OBJECTID")) as cursor:
for row in cursor:
near_fid = row[1]
near_dist = row[2]
near_x = row[3]
near_y = row[4]
objectId = row[5]
id = str(objectId)
intersect_points[near_fid].add(
(near_x,
near_y,
near_dist,
id,
objectId))
del row
del cursor
# Create Road Network
vertices = set()
edges, weights = [], []
fields = ["FNODE_", "TNODE_", "SHAPE@", "LENGTH", "OBJECTID_1"]
with arcpy.da.SearchCursor(street_layer, fields) as cursor:
for row in cursor:
s = str(int(row[0])) # Start node id
t = str(int(row[1])) # End node id
line = row[2] # Polyline object of road
length = row[3] # Road length
lineid = row[4] # Road id
vertices.add(s)
vertices.add(t)
# Check if road near to some data points
if lineid in intersect_points:
# For each intersection point
for ip in intersect_points[lineid]:
ip_near_distance = ip[2]
ip_source_id = ip[3]
intersect_point = arcpy.Point(ip[0], ip[1])
last_point = None
found = False
length_from = 0.0
length_to = 0.0
# Check for each line segment if it intersects
for part in line:
for current_point in part:
current_point_pg = arcpy.PointGeometry(current_point)
# Check if intersection point is on line point
if current_point_pg.contains(intersect_point):
found = True
length_from += ip_near_distance
length_to = ip_near_distance
if last_point is not None:
# Create line segment
array = arcpy.Array()
array.add(last_point)
array.add(current_point)
last_point_pg = arcpy.PointGeometry(last_point)
segment = arcpy.Polyline(array)
# Check if intersection point is on line segment
if not found and segment.contains(intersect_point):
found = True
length_from += ip_near_distance + \
last_point_pg.distanceTo(intersect_point)
length_to = ip_near_distance + \
current_point_pg.distanceTo(
intersect_point)
last_point = current_point
continue
if found:
length_to += segment.length
else:
length_from += segment.length
last_point = current_point
if not found:
print "Int point not found!"
# If interection point found road, add new node to
# network
else:
vertices.add(ip_source_id)
# Add edge from start node to data point
edges.append((s, ip_source_id))
weights.append(length_from)
# Add edge from data point to end node
edges.append((ip_source_id, t))
weights.append(length_to)
# Add road to network
e = (s, t)
w = length
edges.append(e)
weights.append(w)
del row
del cursor
# Create road network graph
g = igraph.Graph()
for v in vertices:
g.add_vertex(v)
for e in edges:
g.add_edge(e[0], e[1])
g.es["weight"] = weights
# Find location points
sources = # List of vertices ids
targets = # List of vertices ids
# Find shortest paths
shortest_paths = g.shortest_paths(
sources,
targets,
weights="weight",
mode="ALL")
