1

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])

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)
| improve this answer | |
  • That does indeed work for the sample, but it doesn't when I use real data. I suppose this is because of floating point calculations – Kenneth Mar 26 at 9:48
  • I've tried to time the code as well and it seems my approach is a bit faster as well. Thanks for the answer though! – Kenneth Mar 26 at 9:54
  • @Kenneth: Good work finding a solution! I'd be surprised if that approach is faster in the general case (mainly because of the loop executing line.project for all coordinates between the start of the line and the first cut) but if it works for your application that's great. – Jake Mar 26 at 10:09
1

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(multi_line, 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(multi_line, point_a, point_b)
    else:
        return get_section_between_distances(multi_line, point_b, point_a)

section = find_section(ml_string, a, b)
| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.