Selection by Angel and Distance

PeterBillig

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

JohannesLindner

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.