Shapely Split LineStrings at Intersections with other LineStrings

Question

I have a set of LineStrings which are intersected by other LineStrings and I want to split the LineString into separate segments at these intersection points. I have a solution but I don't think it's the best approach.

Let's say we are dealing with one LineString:

>>> import shapely
>>> from shapely.geometry import *
>>> import geopandas as gpd
>>> 
>>> MyLine=LineString([(0,0),(5,0),(10,3)])
>>> MyLine
<shapely.geometry.linestring.LineString object at 0x1277EEB0>
>>> 

And 2 lines that intersect this LineString:

>>> IntersectionLines=gpd.GeoSeries([LineString([(2,1.5),(3,-.5)]), LineString([(5,.5),(7,.5)])])
>>> IntersectionLines
0    LINESTRING (2 1.5, 3 -0.5)
1     LINESTRING (5 0.5, 7 0.5)
dtype: object
>>> 

It looks like this:

I can get the intersection points as follows:

>>> IntPoints=MyLine.intersection(IntersectionLines.unary_union)
>>> IntPointCoords=[x.coords[:][0] for x in IntPoints]
>>> IntPointCoords
[(2.75, 0.0), (5.833333333333333, 0.5)]
>>> 

And then I get extract the start and end points and create pairs with these and the intersection points that will be used to form line segments:

>>> StartPoint=MyLine.coords[0]
>>> EndPoint=MyLine.coords[-1]
>>> SplitCoords=[StartPoint]+IntPointCoords+[EndPoint]
>>> 
>>> def pair(list):
...     for i in range(1, len(list)):
...         yield list[i-1], list[i]
... 
>>> 
>>> SplitSegments=[LineString(p) for p in pair(SplitCoords)]     
>>> SplitSegments
[<shapely.geometry.linestring.LineString object at 0x127F7A90>, <shapely.geometry.linestring.LineString object at 0x127F7570>, <shapely.geometry.linestring.LineString object at 0x127F7930>]
>>> 

>>> gpd.GeoSeries(SplitSegments)
0                      LINESTRING (0 0, 2.75 0)
1    LINESTRING (2.75 0, 5.833333333333333 0.5)
2      LINESTRING (5.833333333333333 0.5, 10 3)
dtype: object
>>> 

However, one issue I have with my approach is that I know which intersection point should be joined with the LineString start point and which intersection point should be paired with the LineString end point. This program works if the intersection points are listed along the line in the same order as the start and the ending point. I imagine there would be situations where this would not always be the case? Is there a way to generalize this approach or is there a better approach entirely?

Lutenist

Answer 1

Here is a more general way: calculating the distance along the line for all points (start and end point of the line + points where you want to split), sort by these points and then generate the line segments in the right order. Together in a funtion:

def cut_line_at_points(line, points):

    # First coords of line (start + end)
    coords = [line.coords[0], line.coords[-1]]

    # Add the coords from the points
    coords += [list(p.coords)[0] for p in points]

    # Calculate the distance along the line for each point
    dists = [line.project(Point(p)) for p in coords]

    # sort the coords based on the distances
    # see https://mcmap.net/q/63925/-sorting-list-according-to-corresponding-values-from-a-parallel-list-duplicate
    coords = [p for (d, p) in sorted(zip(dists, coords))]

    # generate the Lines
    lines = [LineString([coords[i], coords[i+1]]) for i in range(len(coords)-1)]

    return lines

Applying this function on your example:

In [13]: SplitSegments = cut_line_at_points(MyLine, IntPoints)

In [14]: gpd.GeoSeries(SplitSegments)
Out[14]:
0                      LINESTRING (0 0, 2.75 0)
1    LINESTRING (2.75 0, 5.833333333333333 0.5)
2      LINESTRING (5.833333333333333 0.5, 10 3)
dtype: object

The only thing is that this does not preserve the corner from your original line (but your example in the question does this neither, so I don't know if this is a requirement. It would be possible but make it a bit more complex)

---

Update A version that keeps the corners in the original line intact (my approach is to also keep a list of 0/1 that indicates if a coord is to be split or not):

def cut_line_at_points(line, points):

    # First coords of line
    coords = list(line.coords)

    # Keep list coords where to cut (cuts = 1)
    cuts = [0] * len(coords)
    cuts[0] = 1
    cuts[-1] = 1

    # Add the coords from the points
    coords += [list(p.coords)[0] for p in points]    
    cuts += [1] * len(points)        

    # Calculate the distance along the line for each point    
    dists = [line.project(Point(p)) for p in coords]    
​
    # sort the coords/cuts based on the distances    
    # see https://mcmap.net/q/63925/-sorting-list-according-to-corresponding-values-from-a-parallel-list-duplicate    
    coords = [p for (d, p) in sorted(zip(dists, coords))]    
    cuts = [p for (d, p) in sorted(zip(dists, cuts))]          

    # generate the Lines    
    #lines = [LineString([coords[i], coords[i+1]]) for i in range(len(coords)-1)]    
    lines = []        

    for i in range(len(coords)-1):    
        if cuts[i] == 1:    
            # find next element in cuts == 1 starting from index i + 1   
            j = cuts.index(1, i + 1)    
            lines.append(LineString(coords[i:j+1]))            

    return lines

Applied on the example:

In [3]: SplitSegments = cut_line_at_points(MyLine, IntPoints)

In [4]: gpd.GeoSeries(SplitSegments)
Out[4]:
0                           LINESTRING (0 0, 2.75 0)
1    LINESTRING (2.75 0, 5 0, 5.833333333333333 0.5)
2           LINESTRING (5.833333333333333 0.5, 10 3)
dtype: object

Coimbatore

Thank, you this is great. You raise a good point about preserving the corner from the original line - I had not thought about that, but you're right it should be a requirement. I took a shot at adapting your code so that the corner is included. It's easy to get the corner included if it's treated as two separate segments (just make it so that the coords variable includes all coordinates rather than just the end points), but keeping that as one segment together seems more complicated. My approach seems a bit kludgy to me - do you mind having a look? – Lutenist

If I understand the question correctly, the same can be achieved by the split function: split(line, intersection_lines) where intersection_lines is a MultiLineString. This was added only in Shapely v1.6.0 (2017-08-21). Would you like to update your answer? – Stratosphere

Answer 2

I love joris's approach. Unfortunately, I ran into a critical difficulty when trying to use it: if the linestring has two points at the same coordinates, their projections are ambiguous. Both will get the same projection value and be sorted together.

This is especially obvious if you have a path that begins and ends at the same point. The ending point gets a projection of 0 and gets sorted at the start, and that throws the entire algorithm off since it expects a "cuts" value of "1" at the end.

Here's a solution that works in shapely 1.6.1:

import shapely.ops
from shapely.geometry import MultiPoint

def cut_linestring_at_points(linestring, points):
    return shapely.ops.split(linestring, MultiPoint(points))

Yes, it really is that simple. The catch here is that the points must be exactly on the line. If they're not, snap them to the line as in this answer.

The return value is a MultiLineString, and you can get at the component LineStrings using its geoms method.

Nevadanevai

Answer 3

And here is my attempt to adapt the function by joris so that the corner of the line segment is included as well. This is does not yet work perfectly because in addition to including the merged segment that includes the corner, it also includes the original unmerged segment.

def cut_line_at_points(line, points):

    #make the coordinate list all of the coords that define the line
    coords=line.coords[:]
    coords += [list(p.coords)[0] for p in points]

    dists = [line.project(Point(p)) for p in coords]

    coords = [p for (d, p) in sorted(zip(dists, coords))]

    lines = [LineString([coords[i], coords[i+1]]) for i in range(len(coords)-1)]

    #Now go through the lines and merge together as one segment if there is no point interrupting it
    CleanedLines=[]      
    for i,line in enumerate(lines):
        if i<>len(lines)-1:
            LinePair=[line,lines[i+1]] 
            IntPoint= LinePair[0].intersection(LinePair[1])
            if IntPoint not in points:
                CleanedLine=shapely.ops.linemerge(LinePair)
            else:
                CleanedLine=line
        else:
            CleanedLine=line


        CleanedLines.append(CleanedLine)
    return CleanedLines

>>> SplitSegments = cut_line_at_points(MyLine, IntPoints)
>>> gpd.GeoSeries(SplitSegments)
0                           LINESTRING (0 0, 2.75 0)
1    LINESTRING (2.75 0, 5 0, 5.833333333333333 0.5)
2            LINESTRING (5 0, 5.833333333333333 0.5)
3           LINESTRING (5.833333333333333 0.5, 10 3)
dtype: object
>>> 

Lutenist

I adopted my version slightly to include the corner point correctly, see my answer. – Coimbatore

Answer 4

@joris's approach is very nice but it bugs out if you try to pass it a list of points, some of which not actually intersect the line, which in my case happens because I pre-calculate a list of intersection points from a list of many lines.

I was able to fix it by pre-filtering the input point list down to only those that actually intersected before proceeding with the function. It is not going to be efficient for large lists of points, but in my case, my lists are always rather small so it was good enough for me. It also works if there are no points that intersect the line, and will just short-circuit to returning the original line as a list in that case (for the sake of consistency)

I initially used line.intersects(point) but it always returned False, probably due to interpolation precision.

def cut_line_at_points(line, points):

    # Filter out any points that are not on the line
    # 0.01 is arbitrary, make it smaller for more precision
    points = [point for point in points if line.distance(point) < 0.01]
    if not points:
        return [line]

    # First coords of line
    coords = list(line.coords)

    # Keep list coords where to cut (cuts = 1)
    cuts = [0] * len(coords)
    cuts[0] = 1
    cuts[-1] = 1

    # Add the coords from the points
    coords += [list(p.coords)[0] for p in points]
    cuts += [1] * len(points)

    # Calculate the distance along the line for each point
    dists = [line.project(Point(p)) for p in coords]

    # sort the coords/cuts based on the distances
    # see https://mcmap.net/q/63925/-sorting-list-according-to-corresponding-values-from-a-parallel-list-duplicate
    coords = [p for (d, p) in sorted(zip(dists, coords))]
    cuts = [p for (d, p) in sorted(zip(dists, cuts))]

    # generate the Lines
    # lines = [LineString([coords[i], coords[i+1]]) for i in range(len(coords)-1)]
    lines = []

    for i in range(len(coords) - 1):
        if cuts[i] == 1:
            # find next element in cuts == 1 starting from index i + 1
            j = cuts.index(1, i + 1)
            lines.append(LineString(coords[i:j + 1]))

    return lines

Staff