How to iterate through vertices in line segments, and all possible branches in arcpy?

Question

I have an interesting problem. I am using arcpy with arcGIS desktop 10.2.

I am given a starting point and distance to a fault on electrical lines. My task is to automate the process of plotting these points into a feature class and displaying them on a webmap for other non-GIS users to view. An example: I receive tabular info that has these two relevant pieces of information:

Starting point: Breaker #113 (Which I know the location of, and touches an electric line segment)
Fault Location: 18.3 Miles

It is important to note that the line features that I must be working on have a commonly shared attribute, called 'lineGroup', but there are many touching line features or segments that may be in that lineGroup. After I select the first line feature that touches my starting point, I iterate recursively through each vertice and measure the distance between them. I keep a running total of the distance, and the code plots a point when it reaches the fault location.

I can get my code to work fine when there is only a single line segment to iterate over, but my dataset contains many different line segments that I may have to iterate over to reach the fault distance. Also, there may be branches off of a line that head in a different direction. So I can potentially have multiple points plotted from one run of this script.

Here is the code I currently have, which works well for a single line segment. What I need help with is how to continue onward with the function when I hit the end of the first line segment, as well as how to branch off and find a possible fault location down a branch.

import arcpy, sys, os, math, time
from math import radians, sin, cos
print("Running against: {}".format(sys.version))
def main_funct(breaker_number, faultMiles):
    try:
        arcpy.env.overwriteOutput = True     
        faultDistance = faultMiles * 1609.34 # Mile to meter conversion rate       
        # Set feature class variables
        ohLine='Database Connections\T_OverheadTransmissionLine'
        breakers= r'C:\GIS\Data\temp_data\temp.gdb\Breakers' # Local copy for testing    
        fault_Locations= r'C:\GIS\Data\temp_data\temp.gdb\Fault_Locations'     
        # Make Feature layers for selections
        breaker_FL = 'in_memory\\breaker_FL'
        arcpy.MakeFeatureLayer_management(breakers, breaker_FL)   
        line_FL = 'in_memory\\line_FL'
        arcpy.MakeFeatureLayer_management(ohLine, line_FL)        
        # Select Starting Breaker and overhead line segments then get breaker XY 
        selectBreakerAndLine(breaker_FL, line_FL, breaker_number)        
        # Measure and total up distance between all vertices in the line segment
        totalDistance = 0
        startingGeometry = getPointXY(breaker_FL, arcpy.Describe(line_FL).spatialReference)
        (startPoint, endPoint, distance2EndPoint) = measureBetweenVertices(faultDistance, line_FL, startingGeometry)
        distance = faultDistance - distance2EndPoint
        angle = findAngleBetweenPoints(startPoint, endPoint)
        faultCoordinates = findPointCoordinates(startPoint, angle, distance)
        updateFaultLocations(fault_Locations, faultCoordinates, breaker_number, faultMiles, line_FL)                 
    except arcpy.ExecuteError:
        print arcpy.GetMessages(2)
    except Exception as e:
        print e.args[0]
    arcpy.AddMessage('\nCompleted')
def selectBreakerAndLine(breaker_FL, line_FL, breaker_number):
    expression = "EquipmentID = '{}'".format(breaker_number)
    arcpy.SelectLayerByAttribute_management(breaker_FL, 'NEW_SELECTION', expression)
    if int(arcpy.GetCount_management(breaker_FL).getOutput(0)) == 1: #Alter this later? For breakers on both ends of line with same ID
        arcpy.SelectLayerByLocation_management(line_FL, 'INTERSECT', breaker_FL)
        #print 'Count of selected Lines: {}'.format(int(arcpy.GetCount_management(line_FL).getOutput(0)))
    else:
        print('Error in selecting breakers')
def getPointXY(startPoint, spatial_ref):
    # Returns the X and Y of the starting point's XY
    shape = arcpy.da.SearchCursor(startPoint, ("SHAPE@",), where_clause=None, spatial_reference= spatial_ref ).next()[0]  
    point = arcpy.Point(shape.centroid.X, shape.centroid.Y)
    geometry = arcpy.PointGeometry(point)
    return geometry    
def measureBetweenVertices(faultDistance, line_FL, startingGeometry, skipList = [0] , totalDistance = 0, previousVertice = None):
    if totalDistance > faultDistance:
        print('Reached the fault location distance\nStart Point: ({}, {})\nEnd Point: ({}, {})\nTotal Distance: {}'.format(previousVertice.centroid.X, previousVertice.centroid.Y,\
            startingGeometry.centroid.X, startingGeometry.centroid.Y, totalDistance))    
        return (startingGeometry, previousVertice, totalDistance)                
    else:    
        skipList.append((startingGeometry.centroid.X, startingGeometry.centroid.Y))
        for row in arcpy.da.SearchCursor(line_FL, ['SHAPE@',], where_clause=None, spatial_reference= arcpy.Describe(line_FL).spatialReference):
            for part in row[0]:
                verticeDict = {}
                for pnt in part:
                    point = arcpy.Point(pnt.X, pnt.Y)
                    verticePoint = arcpy.PointGeometry(point)
                    measuredDistance = verticePoint.distanceTo(startingGeometry)
                    #print( 'Measured distance from start to end point = {}'.format(measuredDistance))
                    if measuredDistance == 0: 
                        pass
                    elif (verticePoint.centroid.X, verticePoint.centroid.Y) in skipList:
                        pass
                    else: 
                        verticeDict[verticePoint] = measuredDistance 
                totalDistance += (min(verticeDict.values()))
                closestVertice = (min(verticeDict.iterkeys(), key=(lambda key: verticeDict[key])))
                skipList.append((closestVertice.centroid.X, closestVertice.centroid.Y))   
        #print(skipList)
        print('Total Distance = {} meters, Required distance = {} meters | Original vertice location = ({}, {}) | Closest vertice location = ({}, {})\n{}'.format(totalDistance, faultDistance,\                                                                                                                                                       startingGeometry.centroid.X, startingGeometry.centroid.Y,\                                                                                                                                                       closestVertice.centroid.X, closestVertice.centroid.Y, '-'*200))
        # Recursive call
        return measureBetweenVertices(faultDistance, line_FL, closestVertice, skipList, totalDistance, startingGeometry)
def findIntersectingLines(startingGeometry, closestVertice, line_FL):
    #Make blank geometry points and an array to store them in for line vertices
    point = arcpy.Point()
    array = arcpy.Array()    
    #Add 1st point to array
    point.X = startingGeometry.centroid.X
    point.Y = startingGeometry.centroid.Y
    array.add(point)    
    #Add 2nd point to array
    point.X = closestVertice.centroid.X
    point.Y = closestVertice.centroid.Y
    array.add(point)    
    #Create a polyline object from array
    polyline = arcpy.Polyline(array)   
    arcpy.SelectLayerByLocation_management(line_FL, 'INTERSECT', polyline)
    print int(arcpy.GetCount_management(line_FL).getOutput(0))
def findAngleBetweenPoints(point1, point2):
    x1, x2 = point1.centroid.X, point2.centroid.X
    y1, y2 = point1.centroid.Y, point2.centroid.Y
    yDelta = y2 - y1
    xDelta = x2 - x1
    radians = math.atan2(yDelta, xDelta)
    degrees = math.degrees(radians)
    adjusted_degrees = (degrees + 360) % 360
    return adjusted_degrees - 90
def findPointCoordinates(startPoint, angle, distance):
    #Find point from an angle and a distance
    origin_x = startPoint.centroid.X
    origin_y = startPoint.centroid.Y    
    # calculate offsets with light trig
    (disp_x, disp_y) = (distance * sin(radians(angle)),\
                        distance * cos(radians(angle)))
    (end_x, end_y) = (origin_x + disp_x, origin_y + disp_y)
    return (end_x, end_y)
def updateFaultLocations(fault_Locations, faultCoordinates, breakerNum, distanceToFault, line_FL):
# Create an edit session and use an updateCursor to insert a row with point geometry for new breaker
    edit = arcpy.da.Editor(os.path.dirname(fault_Locations))
    edit.startEditing(False, True)
    edit.startOperation()                
    updateCursor = arcpy.da.InsertCursor(fault_Locations, ["SHAPE@XY", 'X', 'Y', 'BreakerNumber', 'DistanceFromBreaker',],)
    updateCursor.insertRow([faultCoordinates, faultCoordinates[0], faultCoordinates[1], breakerNum, distanceToFault])
    edit.stopOperation()
    edit.stopEditing(True)
    arcpy.Snap_edit(fault_Locations, [[line_FL, "EDGE", "5 Meters"]])
    return   
# This test allows the script to be used from the operating
# system command prompt (stand-alone), in a Python IDE, 
# as a geoprocessing script tool, or as a module imported in
# another script
if __name__ == '__main__':
    # Arguments are optional
    argv = tuple(arcpy.GetParameterAsText(i)
        for i in range(arcpy.GetArgumentCount()))
    #main_funct(*argv)
    main_funct(715, 18.752) #Test to show that it works

Answer 1

This isn't a direct answer to the question, but it may help. I'm not an engineer and could be getting some details wrong.

I'm surprised you're being given a distance and not a measured fault current. The way I understand locating electrical faults on a distribution system is that it's a function of measured fault current at a point (the breaker) and the impedance of the downline wire. For the same fault current, you'll get different distances based on the wire type, and the wire type will probably vary depending on the path you take from the breaker down through the branches to the possible fault locations.

The engineering, mapping, and outage management system we use (Milsoft's Windmil/Windmilmap) calculates the maximum available fault current for each segment of of wire and each protective device. The control on the breaker gives the measured fault current, and the Milsoft programs can predict possible fault locations based on the stored available fault current values.

It might work to pre-calculate a distance-from-the-source value and store it in your wire feature classes. You could consider weighting the distance calculation based on the wire type if you have some engineering support. Then instead of using the geometric distance, you'd use the pre-calculated (and regularly recalculated) distance value. I think that would simplify your approach at the cost of some maintenance overhead.