selection_lines = "SelectLines" # path or layer name of the lines used to select
test_lines = "TestLines" # path or layer name of the lines to test for perpendicularity
angle_bounds = [70, 110] # what angles do you consider to be perpendicular?
distance = 200 # what's your search range?
output_location = "memory" # where do you want to save the result?
output_name = "PerpendicularLines" # what do you want to call the result?
import arcpy, math
def get_longest_segment(polyline):
"""Returns an arcpy.Polyline object representing the longest segment
of the input arcpy.Polyline
"""
# get line segments
sr = polyline.spatialReference
segments = []
for part in polyline:
for v1, v2 in zip(part, part[1:]):
segment = arcpy.Polyline(arcpy.Array([v1, v2]), spatial_reference=sr)
segments.append(segment)
# sort by length
segments.sort(key=lambda s: s.length)
# return last (longest) segment
return segments[-1]
def get_angle_between_lines(line1, line2):
"""Returns the angle between two arcpy.Polylines.
angle in degrees between 0° and 180°
"""
# get start and end points of both lines
fp1, lp1 = line1.firstPoint, line1.lastPoint
fp2, lp2 = line2.firstPoint, line2.lastPoint
# convert lines to vectors [dx, dy]
v1 = [lp1.X - fp1.X, lp1.Y - fp1.Y]
v2 = [lp2.X - fp2.X, lp2.Y - fp2.Y]
# get dot product and magnitudes
dp = v1[0] * v2[0] + v1[1] * v2[1]
m1 = math.sqrt(v1[0]**2 + v1[1]**2)
m2 = math.sqrt(v2[0]**2 + v2[1]**2)
# get the angle in radians
a = math.acos(dp/(m1 * m2))
# convert to degrees, cap to (0, 180) and return
return math.degrees(a) % 180
# create output table
arcpy.env.addOutputsToMap = True
output_table = arcpy.management.CreateTable(output_location, output_name)
arcpy.management.AddField(output_table, "FID_Selection", "LONG")
arcpy.management.AddField(output_table, "FID_Test", "LONG")
arcpy.management.AddField(output_table, "Angle", "FLOAT")
arcpy.env.addOutputsToMap = False
# create buffer around the selection lines
selection_buffer = arcpy.analysis.Buffer(selection_lines, "memory/SelectionBuffer", distance)
# intersect that buffer with the test lines
selection_test_intersect = arcpy.analysis.Intersect([selection_buffer, test_lines], "memory/SelectionTestIntersect", "ONLY_FID")
# read the pairs of (selection_OID, test_OID) from that intersection
selection_test_pairs = list({row[-2:] for row in arcpy.da.SearchCursor(selection_test_intersect, ["*"])})
arcpy.env.addOutputsToMap = True
# read lines as dictionaries {oid: geometry}
selection_shapes = {oid: shp for oid, shp in arcpy.da.SearchCursor(selection_lines, ["OID@", "SHAPE@"])}
test_shapes = {oid: shp for oid, shp in arcpy.da.SearchCursor(test_lines, ["OID@", "SHAPE@"])}
# start writing into the output table
with arcpy.da.InsertCursor(output_table, ["FID_Selection", "FID_Test", "Angle"]) as cursor:
# loop through the selection-test pairs
for sel_oid, test_oid in selection_test_pairs:
# get longest segments of both lines
sel_segment = get_longest_segment(selection_shapes[sel_oid])
test_segment = get_longest_segment(test_shapes[test_oid])
# get angle between segments
angle = get_angle_between_lines(sel_segment, test_segment)
# if perpendicular, write into the output table
if angle_bounds[0] <= angle <= angle_bounds[1]:
cursor.insertRow([sel_oid, test_oid, angle])
This script will output a table with the following fields:
- FID_Selection: the OBJECTID of the line in your selection layer
- FID_Test: the OBJECTID of the line in your test layer
- Angle: the angle between the longest segments of those lines
To run it:
- Open the Python Window
- Copy and paste the script into the Python Window
- Edit the variables at the start
- If you use layers for selection_lines and test_lines, make sure there is no selection.
- Use forward slashes for paths
- output_location = "C:/SomeFolder/Database.gdb"
- using "memory" writes into RAM, which is faster (especially for large datasets), but the data is lost when you close ArcGIS. Don't forget to export the table!
- Hit Enter twice.
Have a great day!
Johannes