Issues using arcpy for network analyst tasks

[GregKing1]

Hi,

I am writing a script with ArcPy to look at modelling closed
roads and the length of possible diversions.  This is being done with the network analyst
tools.

The basic process I am using is to:

• Create a new route analysis
• Add route end points and a polygon barriers to
  block the road in question
• Solve the problem to create the result of an alternative
  route to the closed road.
• Save this result to a geodatabase feature class.

I have created a script that works well for a single
polyline feature.  However, when I put
this code into a ‘for loop’ to run over multiple polylines, I get errors.  The resulting feature class has route results
for all of the polylines in my input feature class.  I don’t get a separate feature class for each
time the loop iterates through a separate polyline in my input feature class.

I think it may have something to do with removing the end
points and barrier polygons from my route layer after each iteration of the
loop.

Can someone have a quick look at this code and let me know
what the issue may be?

Thanks

Greg

import arcpy
from arcpy import na
from arcpy import env

env.overwriteOutput = True
arcpy.CheckOutExtension('Network')

#___set global variables__

# road segments layer
roadSegs = r'Z:\BNE\Projects\246000\246777-00 Prioritising road segments\Work\Internal\GIS\Data\ALL_SEGMENTS_PROCESSING\Diversions_Assessment.gdb\State_Road_Network_Results\Sample_Segment_Five_Digit_IDs'
# network dataset (local road network, state road network, freight road network).  Needs to have a 'time' attribute for each link
networkDS = r'Z:\BNE\Projects\246000\246777-00 Prioritising road segments\Work\Internal\GIS\Data\ALL_SEGMENTS_PROCESSING\Diversions_Assessment.gdb\State_Road_Network\Qld_State_Road_Network_CL_ND'
# temp geodatabase for working layers
tempFGDB = r'Z:\BNE\Projects\246000\246777-00 Prioritising road segments\Work\Internal\GIS\Data\ALL_SEGMENTS_PROCESSING\Diversions_TEMP.gdb'
# output route line object
routeFC = r'Z:\BNE\Projects\246000\246777-00 Prioritising road segments\Work\Internal\GIS\Data\ALL_SEGMENTS_PROCESSING\Diversions_Assessment.gdb'

#Loop through each segment in the feature class to calculate the diversion
for segment in roadSegs:
   
    # Get the segment ID value from the segment line layer
    segIDCursor = arcpy.UpdateCursor(roadSegs)

    for segRow in segIDCursor:
        segID = segRow.getValue('Segment_ID')

    #___create the working layers for the assessment___

    # Local loop variables
    segsEnds = tempFGDB + '\RoadSegsEndPts_'+ str(segID)[:5]
    endsBuffer = segsEnds + '_Buffer100m_'+ str(segID)[:5]
    segsErase = tempFGDB + '\RoadSegsErase_'+ str(segID)[:5]
    segsEraseBuffer = tempFGDB + '\SegmentClosurePolygon_'+ str(segID)[:5]

    # Convert road segment to points (end)
    arcpy.FeatureVerticesToPoints_management(roadSegs,segsEnds,'BOTH_ENDS')

    # Create the blockage geometry:
    # Buffer the end points by 100m
    arcpy.Buffer_analysis(segsEnds,endsBuffer,'100 Meters')
    # Erase road segment line by 100m point buffers
    arcpy.Erase_analysis(roadSegs,endsBuffer,segsErase)
    # Buffer the erased road segment line by 50m - POLYGON BARRIER
    arcpy.Buffer_analysis(segsErase,segsEraseBuffer,'50 Meters')

    #___setup and solve the route scenario layer___

    # MakeRouteLayer tool
    routeLayer = arcpy.MakeRouteLayer_na(networkDS,'ReRoute','Time')
    routeOutput = routeLayer.getOutput(0)
    subLayerNames = arcpy.na.GetNAClassNames(routeOutput)
    routesLayerName = subLayerNames["Routes"]
    # AddLocations - End Points
    arcpy.AddLocations_na(routeLayer,'Stops',segsEnds,'','100 Meters','','SHAPE','','CLEAR','SNAP','','','')
    # AddLocations - barrier polygon
    arcpy.AddLocations_na(routeLayer,'Polygon Barriers',segsEraseBuffer,'', '#','','SHAPE','','','','','INCLUDE','')
    # Solve the route scenario
    arcpy.Solve_na(routeLayer,'SKIP','CONTINUE','')
    # Save the result as a Feature Layer (set the file name as the 'segment ID' + 'network type')
    RoutesSubLayer = arcpy.mapping.ListLayers(routeOutput, routesLayerName)[0]
    arcpy.management.CopyFeatures(RoutesSubLayer,routeFC + '\State_Road_Network_Results\Deviation_Result_Segment_' + str(segID)[:5])

    # variable for the resulting deviation layer.  SEE IF THIS CAN BE INCLUDED IN THE LOCAL VARIABLES ABOVE
    routeDeviation = routeFC + '\State_Road_Network_Results\Deviation_Result_Segment_' + str(segID)[:5]
    #routeDeviation.strip('.0')

    print routeDeviation

    # Add a field to the deviation layer and add in the segment ID from the segment line layer
    arcpy.AddField_management(routeDeviation,'Seg_ID','TEXT','','',50)
    deviateCursor = arcpy.UpdateCursor(routeDeviation)
    for devRow in deviateCursor:
        devRow.setValue('Seg_ID',segID)
        deviateCursor.updateRow(devRow)

    # Delete the temporary working files and objects
    arcpy.Delete_management('im_memory')
   
    del segIDCursor
    del segRow

    del segsEnds
    del endsBuffer
    del segsErase
    del segsEraseBuffer

    del deviateCursor
    del devRow

# ___Close Loop___

arcpy.CheckInExtension('Network')

[NeilAyres]

Firstly,

use the advanced editor and post the code with python syntax highlighting.

That might make it a bit easier to read.

Those paths are very complex and contain spaces and other nasty bits. Have you tried running this locally instead?

What is happening here:

#Loop through each segment in the feature class to calculate the diversion
for segment in roadSegs:

    # Get the segment ID value from the segment line layer
    segIDCursor = arcpy.UpdateCursor(roadSegs)
    for segRow in segIDCursor:
        segID = segRow.getValue('Segment_ID')

What is the nature of roadSegs? This does not seem to be a list, it points to something on disk (I think).

Then you open an Update cursor, but actually just reads the data.

Then it reads everything, so segID will always be the last one in segIDCursor..

[GregKing1]

Hi Neil,

Thanks for your reply.

When I was initially testing this code I was using a single polyline feature.  The code works well for a single feature and I get the result I am looking for.  So the file paths I have set here work fine in this situation.

The issue is when I put the main body of code into the 'for loop' I'm not getting the results I thought I would be.  My input feature class has multiple polylines representing road centrelines.  What I would like to get is a separate feature class with the result from the network analyst operation for each feature in the input feature class .

'roadSegs' is the feature class with the road polyline features.  Within this feature class is an attribute called 'Segment_ID' with the ID value for each road polyline.

Following your comment below, how can I setup the cursor so it only reads the ID value for the polyline being processed in each particular iteration of the loop?

Also, when I run a network analyst operation using python, should I delete the network dataset at the end of each iteration of the loop?  Then create a new network dataset at the start of the next iteration?

thanks for your assistance.

Greg

[NeilAyres]

Havn't been through your code, it is rather difficult to follow.

Consider using arcpy.da style cursor to access data.

If you need the separate features of roadSeg and to process them individually you could first read them into a dictionary.

Like

roadSegDict = {}
with arcpy.da.SearchCursor(fc, ["Segment_ID", "SHAPE@"]) as Cur:
    for row in Cur:
        roadSegDict[row[0]] = row[1]  # row[0] is the Segment_ID, row[1] is the geometry.

In this scenario, Segment_ID must be unique.

Then you can loop through the contents of roadSegDict.

for ID, geom in roadSegDict.iteritems():

When I was last messing with NA datasets, I couldn't find a "Create Network Dataset" tool to use in python.

It maybe there somewhere, but I couldn't find it.

Ended up having a pre created one (using the wizard), then deleting content, appending content, then solving.