Determine which intersection result further along a line in Shapely

Question

Context

(minimum working example code at the end)

I am trying to create a line smoothing algorithm: From the first line vertex, a circle of a specified radius is created and this circle is intersected with the line (LineString). This results into 1 to many intersection results (Points). The farthest intersection point along the line is then taken, a new circle is created, intersected with the line and so on until the line endpoint is approached.

Visualized result with an oddly-passing parameters is in the image below. The blue line is the input, the green line is the smoothed result. The line is in the UTM Projection, but it doesn't matter.

The issue

Now I roughly know how to design most of the algorithm (I can pass it on full here if successful).

What I don't know is how to always choose the farthest intersection point along the line from the current intersection results. Shapely utilizes Point or MultiPoint as an interserciton result in this case. If I reach MultiPoint parts through geoms, they're probably ordered by ascending x coordinate (not sure though).

Attempted solution

I first thought I'd always take the last member from the circle_intersections.geoms. As you probably guess, if I apply this to the algorithm, this wouldn't work once the line changes its direction against the ascending x coordinate. Then, the intersection point I need is usually the first (or other) from the list. And thus, it gets into infinite loop, passing the intersection points between two circles.

In the image above, it worked just because the radius was set such that the searched intersection point was always further in x coordinate.

Secondly, I tried to iterate through the line and filter out the intersection point which splits the line farthest along. Here I must be making some kind of mistake I don't see.

You find both attempts in the find_farthest_intersection function.

Code

Imports:

import shapely
from shapely.geometry import Point, LineString
from shapely.ops import split
import matplotlib.pyplot as plt

The algorithm:

def smooth(point, geometry, r, first_pass=True):
    circle = point.buffer(r)
    circle_intersections = geometry.intersection(circle.exterior)

    # It's always a single point at the first iteration
    if isinstance(circle_intersections, Point):
        next_point = circle_intersections
    else:
        next_point = find_farthest_intersection(
            geometry,
            circle_intersections
        )

    # If you want to visualize in Jupyter
    # plt.plot(
    #     *point.buffer(6).exterior.xy,
    #     *circle.exterior.xy,
    #     *next_point.buffer(5).exterior.xy,
    # )

    if not first_pass and Point(geometry.coords[-1]).within(circle):
        return [geometry.coords[-1]]

    return [point.coords[0]] + smooth(
        next_point,
        geometry,
        r,
        first_pass=False
    )


def find_farthest_intersection(geometry, intersections):
    # How to return the right intersection furthest along the line?

    # Naive and wrong solution
    # When the line heads against X this becomes the opposite point than I want
    # return circle_intersections.geoms[-1]

    # Other attempted solution
    # Not working, results in the same problem as with the solution above.

    # Find each intersection point's distance from the
    # line beginning
    intersection_pnt_dsts_from_line_start = [
        split(geometry, point).geoms[0].length
        for point in intersections.geoms
    ]

    # Return the intersection point with the maximum
    # distance from the line beginning
    return max(
        zip(
            intersection_pnt_dsts_from_line_start,
            intersections.geoms
        ),
        key=lambda item: item[0]
    )[1]

If you want to visualize in eg. Jupyter (also uncomment pyplots in the smooth function).

# it's in UTM btw
line = shapely.wkt.loads('LINESTRING (478098.5964211893 5442495.543688663, 478229.0423690955 5442525.981076508, 478242.0869638861 5442558.592563484, 478198.6049812507 5442595.552248725, 478129.033809034 5442649.904727018, 478168.1675934059 5442691.212610522, 478209.4754769095 5442665.123420941, 478209.4754769095 5442628.163735701, 478213.8236751731 5442615.11914091, 478231.2164682273 5442599.900446988, 478252.957459545 5442593.378149592, 478305.1358387075 5442602.07454612, 478322.5286317617 5442652.078826151, 478329.0509291569 5442732.520494026, 478383.4034074512 5442762.957881871, 478353.5095443893 5442726.541721413, 478421.4501422573 5442696.37609596, 478471.7261846794 5442718.117087278, 478434.7664994393 5442728.987582936, 478399.9809133308 5442732.248731634, 478458.6815898888 5442757.250871649, 478526.0786629738 5442780.078912533, 478532.6009603692 5442848.563035184, 478579.3440917023 5442883.348621292)')

r = 100  # radius parameter

plt.plot(*line.xy)

result = smooth(
    Point(line.coords[0]),  # first time it's the first line node
    line,
    r # smaller number or more complicated lines cause fail
)

# Would then do
# plt.plot(*LineString(result).xy)

score 0 · Answer 1 · answered Nov 13 '21 at 16:27

So I played around more with geometric operations and with help of #50194077 found the solution is rather simple. I can split the line by each intersection point and measure the distance from the line beginning of the first split result. The longest distance points to the farthest intersection point, i.e. the intersection point I want.

The key to facilitate the split is to create a buffer of the intersection point. Its size should addresses shapely precision issues and be as low as possible (eg 1Oe-9).

Revised function:

def find_farthest_intersection(geometry, intersection_points):
    intersection_points_from_geom_distance = []

    for intersection_point in intersection_points:
        # Create a little buffer
        precision_buffer = intersection_point.buffer(10e-9)

        # Split line on buffer and take the first line fragment
        intersection_points_from_geom_distance.append([
            intersection_point,
            split(line, precision_buffer)[0].length]
        )

    # Return the intersection point with the maximum
    # distance from the line beginning
    return max(
        intersection_points_from_geom_distance,
        key=lambda item: item[1]
    )[0]

Determine which intersection result further along a line in Shapely

Context

The issue

Attempted solution

Code

1 Answers1