Hi,
I have two Vector layers of lines. I am trying to find all the lines that are with in a set range of up to 200 meters, from one layer and are Perpendicular to the lines from the other layer (angles from 70-110 degrees). The Idea is to be able to sort the relevant lines by distance and angles.
Unfortunately I am not familiar with coding, if there is a solution that involves coding I will need guidance for it. All help is greatly appreciated.
Thank you, Peter
Solved! Go to Solution.
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:
To run it:
You won't get around a little (might also be a lot) Python for this problem...
I can see what I can whip up when I get some time.
But first, I need you to clarify some things:
Dear Johannes,
Thank you so much for you kind assistance I truly appreciate it. I sincerely apologize for my delayed response.
Regarding your questions:
Both layers are Polylines.
1. We need all close perpendicular lines.
2. We have a polyline, so in the case of a segment with more than a single line, we would like to take the angel of the largest line within the segment for the calculation of perpendicularity.
3. We would like that in such a case it will be present multiple times in the table.
We would also like that the output will contain a new column which will represent the FID of the segment in fc2 the segment in fc1 is perpendicular to.
Once again thank you,
Peter
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:
To run it:
This Is amazing, Thank you so much.