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
Thanks to all!
