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]

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:

1. Do you need all close perpendicular lines or only the closest perpendicular line?
2. Do you have lines (only 2 vertices) or polylines (multiple vertices)? If polylines:
   1. What angle are you looking at? The angle between start and end point? The mean of the angles of every segment? The mean of angles in the 200 meter buffer? Something else?
3. If a line of fc1 is close and perpendicular to multiple lines of fc2, should it be present in the output table multiple times or only once? if only once, which distance and angle should be used?

[PeterBillig]

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

[JohannesLindner · Accepted 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:

• 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.

[PeterBillig]

This Is amazing, Thank you so much.