Shapely Split LineStrings at Intersections with other LineStrings

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?

Solution

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 http://stackoverflow.com/questions/6618515/sorting-list-based-on-values-from-another-list
    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 http://stackoverflow.com/questions/6618515/sorting-list-based-on-values-from-another-list    
    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