Fetch continuous section between two points on a MultiLineString using shapely

Question

Let's say I have a MultiLineString like this:

ml_string = MultiLineString([
   [[0,0], [1,1], [2,2], [3,3]],
   [[0,0], [0,1], [0,2], [0,3]]
])

Given the following two points:

a = Point([1,1])
b = Point([3,3])

Using shapely, wow can I get the first (if any) continuous section in the MultiLineString?

In this case that would be ([1,1], [2,2], [3,3])

Answer 1

You can use Shapely's split function for this:

from shapely.geometry import MultiLineString, Point
from shapely.ops import split

ml_string = MultiLineString([
   [[0,0], [1,1], [2,2], [3,3]],
   [[0,0], [0,1], [0,2], [0,3]]
])

a = Point([1,1])
b = Point([3,3])

section = None
for ls in ml_string:
    if a.intersects(ls) and b.intersects(ls):
        # Sort points by distance along the LineString
        if ls.project(a) < ls.project(b):
            start = a
            end = b
        else:
            start = b
            end = a

        section = split(ls, start)[1]
        section = split(section, end)[0]

print(section)  # LINESTRING (1 1, 2 2, 3 3)

Answer 2

After trying out various methods, I found a solution. The algorithm works as follows:

• Find the first line string that has both points
• Find the distance of both points
• Cut the line at distance A and take the second part
• Cut the second part at distance B - distance A and take the first part

This is the code I ended up using:

def find_line(multi_line, point_a, point_b):
    for line in multi_line:
        has_point_a = line.distance(point_a) < 5
        has_point_b = line.distance(point_b) < 5
        if has_point_a and has_point_b:
            return line

def cut(line, distance):
    if distance <= 0.0 or distance >= line.length:
        return [LineString(line)]
    coords = list(line.coords)
    for i, p in enumerate(coords):
        pd = line.project(Point(p))
        if pd == distance:
            return [
                LineString(coords[:i+1]),
                LineString(coords[i:])]
        if pd > distance:
            cp = line.interpolate(distance)
            return [
                LineString(coords[:i] + [(cp.x, cp.y)]),
                LineString([(cp.x, cp.y)] + coords[i:])]

def get_section_between_distances(line_to_cut, distance_a, distance_b):
    first_cut = cut(line_to_cut, distance_a)[1]
    return cut(first_cut, distance_b - distance_a)[0]

def find_section(multi_line, point_a, point_b):
    line_to_cut = find_line(multi_line, point_a, b)

    distance_a = line_to_cut.project(point_a)
    distance_b = line_to_cut.project(point_b)

    if distance_a > distance_b:
        return get_section_between_distances(line_to_cut, point_a, point_b)
    else:
        return get_section_between_distances(line_to_cut, point_b, point_a)

section = find_section(ml_string, a, b)