Creating a line that goes through points

2370
12
Jump to solution
08-15-2023 01:32 PM
mtorfi
by
Occasional Contributor

Hello, 

I would like to create a line that goes through adjacent points. The condition that I would like to apply to this are the following:

1- if the line that goes through some points change direction, I want a new line to be created

2- I want to know the length of each linesubmittoESRI.PNG

In the picture above, I have shown what I would like to do as an example. I would appreciate the community's help in this matter. 

1 Solution

Accepted Solutions
JohannesLindner
MVP Frequent Contributor

Assuming your points are ordered:

  • Run Points To Line
  • Copy/Paste the script below into your Python Window (without lines 79 and after) and call divide_lines_into_straight_lines
    • alternatively, turn the script into a script tool and run that (also enclosed at end of this post)

 

import arcpy
import math
from pathlib import Path

# function to divide a list of vertices
def divide_path_into_straight_paths(path, angle_tolerance):
    """ Divides an arcpy line path into straight parts.

        path: [arcpy.Point], the vertices of the path
        angle_tolerance: float, difference of consecutive angles (in degrees)
            between vertices smaller than this value are considered to be straight

        Returns the split paths as [[arcpy.Point]]
    """
    segments = [
        [path[i-1], path[i]]
        for i in range(1, len(path))
        ]
    angles = [
        math.degrees(math.atan2(s[0].Y - s[1].Y, s[0].X - s[1].X)) % 180  # % 180 to get the same angle regardless of vertex order
        for s in segments
        ]
    split_paths = []
    split_path_angles = []
    for i, (segment, angle) in enumerate(zip(segments, angles)):
        if i == 0:
            split_paths.append(segment)
            split_path_angles.append([angle])
            continue
        angle_difference = abs(angle - angles[i-1])
        if angle_difference > 90:
            angle_difference = abs(angle_difference - 180)  # to make 179° comparable to 1°
        if angle_difference <= angle_tolerance:
            split_paths[-1].append(segment[-1])
            split_path_angles[-1].append(angle)
        else:
            split_paths.append(segment)
            split_path_angles.append([angle])
    return split_paths, split_path_angles


# function to divide a line geometry
def divide_line_into_straight_lines(line, angle_tolerance):
    """ Divides an arcpy line geometry into straight parts.

        line: arcpy.Polyline, the line geometry
        angle_tolerance: float, difference of consecutive angles (in degrees)
            between vertices smaller than this value are considered to be straight

        Returns the split lines as [arcpy.Polyline]
    """
    split_lines = []
    for part in line:
        split_paths, split_path_angles = divide_path_into_straight_paths(part, angle_tolerance)
        for path in split_paths:
            split_lines.append(arcpy.Polyline(arcpy.Array([path]), spatial_reference=line.spatialReference))
    return split_lines, split_path_angles


# function to divide all lines in a fc and write the output into a new fc
def divide_lines_into_straight_lines(input_lines, output_lines, angle_tolerance):
    """ Divides all polylines in a featureclass or layer (honors selection)
        into straight parts.

        input_lines: str, path to the input featureclass / name of the layer
        output_lines: str, path of the output featureclass
        angle_tolerance: float, difference of consecutive angles (in degrees)
            between vertices smaller than this value are considered to be straight
    """
    lines = [row for row in arcpy.da.SearchCursor(input_lines, ["OID@", "SHAPE@"])]
    out_path = Path(output_lines)
    out_fc = arcpy.management.CreateFeatureclass(str(out_path.parent), str(out_path.stem), "POLYLINE", spatial_reference=lines[0][1].spatialReference)
    arcpy.management.AddField(out_fc, "FID", "LONG")
    arcpy.management.AddField(out_fc, "LENGTH", "DOUBLE")
    arcpy.management.AddField(out_fc, "ANGLE_START", "DOUBLE")
    arcpy.management.AddField(out_fc, "ANGLE_MIN", "DOUBLE")
    arcpy.management.AddField(out_fc, "ANGLE_MAX", "DOUBLE")
    arcpy.management.AddField(out_fc, "ANGLE_RANGE", "DOUBLE")
    arcpy.management.AddField(out_fc, "VERTEX_COUNT", "LONG")
    with arcpy.da.InsertCursor(out_fc, ["FID", "SHAPE@", "LENGTH", "ANGLE_START", "ANGLE_MIN", "ANGLE_MAX", "ANGLE_RANGE", "VERTEX_COUNT"]) as cursor:
        for oid, shp in lines:
            split_shapes, split_angles = divide_line_into_straight_lines(shp, angle_tolerance)
            for shape, angles in zip(split_shapes, split_angles):
                a_start = angles[0]
                a_min = min(angles)
                a_max = max(angles)
                a_range = a_max - a_min
                if a_range > angle_tolerance:
                    a_range = abs(a_range - 180)
                v_count = len(angles) + 1
                cursor.insertRow([oid, shape, shape.length, a_start, a_min, a_max, a_range, v_count])


# you can use the script as tool, parameters:
    # Feature Layer, Required, Input, Filter for Polylines
    # Feature Layer, Required, Output, Filter for Polylines
    # Double, Required, Input
if __name__ == "__main__":
    input_lines = arcpy.GetParameterAsText(0)
    output_lines = arcpy.GetParameterAsText(1)
    angle_tolerance = arcpy.GetParameter(2)
    divide_lines_into_straight_lines(input_lines, output_lines, angle_tolerance)

 

 

Points:

JohannesLindner_0-1692178303041.png

 

Points To Line:

JohannesLindner_1-1692178331036.png

 

 

Straight Lines:

JohannesLindner_2-1692178361996.png

 

JohannesLindner_3-1692178404445.png

 


Have a great day!
Johannes

View solution in original post

0 Kudos
12 Replies
JohannesLindner
MVP Frequent Contributor

Assuming your points are ordered:

  • Run Points To Line
  • Copy/Paste the script below into your Python Window (without lines 79 and after) and call divide_lines_into_straight_lines
    • alternatively, turn the script into a script tool and run that (also enclosed at end of this post)

 

import arcpy
import math
from pathlib import Path

# function to divide a list of vertices
def divide_path_into_straight_paths(path, angle_tolerance):
    """ Divides an arcpy line path into straight parts.

        path: [arcpy.Point], the vertices of the path
        angle_tolerance: float, difference of consecutive angles (in degrees)
            between vertices smaller than this value are considered to be straight

        Returns the split paths as [[arcpy.Point]]
    """
    segments = [
        [path[i-1], path[i]]
        for i in range(1, len(path))
        ]
    angles = [
        math.degrees(math.atan2(s[0].Y - s[1].Y, s[0].X - s[1].X)) % 180  # % 180 to get the same angle regardless of vertex order
        for s in segments
        ]
    split_paths = []
    split_path_angles = []
    for i, (segment, angle) in enumerate(zip(segments, angles)):
        if i == 0:
            split_paths.append(segment)
            split_path_angles.append([angle])
            continue
        angle_difference = abs(angle - angles[i-1])
        if angle_difference > 90:
            angle_difference = abs(angle_difference - 180)  # to make 179° comparable to 1°
        if angle_difference <= angle_tolerance:
            split_paths[-1].append(segment[-1])
            split_path_angles[-1].append(angle)
        else:
            split_paths.append(segment)
            split_path_angles.append([angle])
    return split_paths, split_path_angles


# function to divide a line geometry
def divide_line_into_straight_lines(line, angle_tolerance):
    """ Divides an arcpy line geometry into straight parts.

        line: arcpy.Polyline, the line geometry
        angle_tolerance: float, difference of consecutive angles (in degrees)
            between vertices smaller than this value are considered to be straight

        Returns the split lines as [arcpy.Polyline]
    """
    split_lines = []
    for part in line:
        split_paths, split_path_angles = divide_path_into_straight_paths(part, angle_tolerance)
        for path in split_paths:
            split_lines.append(arcpy.Polyline(arcpy.Array([path]), spatial_reference=line.spatialReference))
    return split_lines, split_path_angles


# function to divide all lines in a fc and write the output into a new fc
def divide_lines_into_straight_lines(input_lines, output_lines, angle_tolerance):
    """ Divides all polylines in a featureclass or layer (honors selection)
        into straight parts.

        input_lines: str, path to the input featureclass / name of the layer
        output_lines: str, path of the output featureclass
        angle_tolerance: float, difference of consecutive angles (in degrees)
            between vertices smaller than this value are considered to be straight
    """
    lines = [row for row in arcpy.da.SearchCursor(input_lines, ["OID@", "SHAPE@"])]
    out_path = Path(output_lines)
    out_fc = arcpy.management.CreateFeatureclass(str(out_path.parent), str(out_path.stem), "POLYLINE", spatial_reference=lines[0][1].spatialReference)
    arcpy.management.AddField(out_fc, "FID", "LONG")
    arcpy.management.AddField(out_fc, "LENGTH", "DOUBLE")
    arcpy.management.AddField(out_fc, "ANGLE_START", "DOUBLE")
    arcpy.management.AddField(out_fc, "ANGLE_MIN", "DOUBLE")
    arcpy.management.AddField(out_fc, "ANGLE_MAX", "DOUBLE")
    arcpy.management.AddField(out_fc, "ANGLE_RANGE", "DOUBLE")
    arcpy.management.AddField(out_fc, "VERTEX_COUNT", "LONG")
    with arcpy.da.InsertCursor(out_fc, ["FID", "SHAPE@", "LENGTH", "ANGLE_START", "ANGLE_MIN", "ANGLE_MAX", "ANGLE_RANGE", "VERTEX_COUNT"]) as cursor:
        for oid, shp in lines:
            split_shapes, split_angles = divide_line_into_straight_lines(shp, angle_tolerance)
            for shape, angles in zip(split_shapes, split_angles):
                a_start = angles[0]
                a_min = min(angles)
                a_max = max(angles)
                a_range = a_max - a_min
                if a_range > angle_tolerance:
                    a_range = abs(a_range - 180)
                v_count = len(angles) + 1
                cursor.insertRow([oid, shape, shape.length, a_start, a_min, a_max, a_range, v_count])


# you can use the script as tool, parameters:
    # Feature Layer, Required, Input, Filter for Polylines
    # Feature Layer, Required, Output, Filter for Polylines
    # Double, Required, Input
if __name__ == "__main__":
    input_lines = arcpy.GetParameterAsText(0)
    output_lines = arcpy.GetParameterAsText(1)
    angle_tolerance = arcpy.GetParameter(2)
    divide_lines_into_straight_lines(input_lines, output_lines, angle_tolerance)

 

 

Points:

JohannesLindner_0-1692178303041.png

 

Points To Line:

JohannesLindner_1-1692178331036.png

 

 

Straight Lines:

JohannesLindner_2-1692178361996.png

 

JohannesLindner_3-1692178404445.png

 


Have a great day!
Johannes
0 Kudos
mtorfi
by
Occasional Contributor

This worked just the way I wanted it to. 

Thank you very much Johannes.

Majid

0 Kudos
mtorfi
by
Occasional Contributor

@JohannesLindner I tried using the tool for when a line goes in two/three/more different directions and it did not work. Is there a way to change the code to do the following: 

ESRI1.PNG

Appreciate your help Johannes. 

Majid

0 Kudos
JohannesLindner
MVP Frequent Contributor

Hmm... The tool has no problem dividing this structure, but it's programmed to work with lines, not points.

If you input this multipart line (it also works for multiple single-part lines):

JohannesLindner_0-1692201366988.png

 

It returns the correct splits:

JohannesLindner_1-1692201467234.png

 

 

The problem then becomes "How to get the points into the correct line shape(s)", which is easy for simple point chains (Points To Line), but might be more complicated here.

Points To Line has a Line Field parameter. If the points have an attribute that differs between the line parts, you can use that.

If you don't have such a field, then you might be stuck doing these corners manually. I could probably rewrite the tool to work with points instead of lines, but I have no idea how I would tackle this case, so it wouldn't do you any good...


Have a great day!
Johannes
0 Kudos
mtorfi
by
Occasional Contributor

What I noticed is that if all the lines are under 1 line, this tool does indeed work. 

But if I have multiple lines that have their own values and GPS points and shape lengths, then it does not. Is there a way to combine the lines with one another? I know relate and join are two tools. I think once we have the bigger line then all of these issues go away. 

Thank you Johannes,

Majid

0 Kudos
mtorfi
by
Occasional Contributor

Also, it seems as though once the problem gets bigger, it is not capable of handling the curves and nicks. 

ESRI2.PNG

I have lines that have different Line Field values but the code does not seem to work. 

Appreciate your help Johannes, 

Majid

0 Kudos
mtorfi
by
Occasional Contributor

Just to give you the sort of errors I am getting, here it is: 

ESRI3.PNG

0 Kudos
JohannesLindner
MVP Frequent Contributor

Just to give you the sort of errors I am getting, here it is: 

Haha, I thought about this possible error when writing the code but decided against writing checks that would make the code more complicated because I thought it wasn't very probable that the error actually occurs... I will fix it later.

 

 Is there a way to combine the lines with one another?

Dissolve without specifying a Dissolve Field should combine all lines into one huge multipart line.

 

it seems as though once the problem gets bigger, it is not capable of handling the curves and nicks

This is something I can't test without actual data. Can you send me your point and/or line feature classes here or in a private message?


Have a great day!
Johannes
JohannesLindner
MVP Frequent Contributor
  • fixed the ZeroDivisionError
  • fixed the behavior for nearly horizontal lines
    • here, the angle can switch from 180°/359° to 0°. this is a small angle difference, but the tool would split the line, because it saw a difference of 179°... now it works correctly.
  • Added some fields in the output fc, mostly for debugging, but maybe it's useful...
  • fixed in script and reuploaded the tool

 

it seems as though once the problem gets bigger, it is not capable of handling the curves and nicks

This is something I can't test without actual data. Can you send me your point and/or line feature classes here or in a private message?

Thanks for the data. After fixing the error, the tool works fine for me. It takes some time (3 to 4 minutes), but that's to be expected, because it does it does expensive geometry computations on 39.000 features...

I took a look at some of the more complex regions, but I couldn't find anywhere where it wasn't "capable of handling the curves and nicks". Can you try again and send me coordinates of places here it doesn't work?

Just to be sure: The Angle Tolerance parameter controls how the "curves and nicks" are handled. If the  difference between angles of two consecutive segments of an input line is bigger than that value (in degrees), then these segments will be in different output lines. If the angle difference is smaller than the tolerance, they will be part of the same output line. So if you want greater differentiation, set the tolerance to a smaller value.

For example, this is the split with an Angle Tolerance of 10°:

JohannesLindner_0-1692302658006.png

 

When I use 50°, there are less splits:

JohannesLindner_1-1692302824454.png

 


Have a great day!
Johannes
0 Kudos