3

Is there a way to "sort" the points of a MultiLineString such that they follow a more sequential order? Let me illustrate what I mean with a graphical example:

Input

Suppose I have the following MultiLineString: MULTILINESTRING ((120 105, 140 98),(160 103, 140 98),(100 100, 120 105)) Input MultiLineString

Note how, to our human eyes, it seems out of order: the first (orange) segment is in the middle, the second (green) segment is coded backwards in relation to the other segments (i.e., pointing to the left instead of to the right).

Desired Output

The desired output would be the following LineString: LINESTRING (100 100, 120 105, 140 98 160 103)

Output LineString

Note how the output is now looks much more sensible. It's one continuous LineString, the segments are all nicely ordered from left to right and the green segment itself is now oriented from left to right.

Is there a simple way to perform this kind of "reordering" inside MultiLineStrings?

What I've done so far

I'm trying to develop a solution using Python's Shapely library. My approach has two different steps:

  • Extract unique coordinates
  • Sort unique coordinates

I've been able to do the first part quite easily, but it's not clear to me how to reorder the points.

I thought about a heuristic that would find the point that is cumulatively furthest to all other points, label that as the starting point, add all other points to a pool of unsorted points and then sequentially find the next closest point until I run out of unsorted points.

However, this sorting heuristic would fail in A BUNCH of cases: consider a geometry that describes a spiral. Or a strong zig-zagging pattern. Or a geometry with a bunch of self-intersections. This heuristic would just fail completely.

Therefore, I suspect this approach of extracting unique points and sorting them isn't the right one.

Back to the start

So, I'm back to my main question: how can I sort the coordinates of a MultiLineString such that they make more intuitive sense?

Reproducible example and code

Here's the Python code for my very incomplete solution thus far.

import shapely

# Write the MultiLineString WKT
mls_wkt = 'MULTILINESTRING ((120 105, 140 98),(160 103, 140 98),(100 100, 120 105))'

# Generate the shapely geometry
mls = shapely.wkt.loads(mls_wkt)

# Function that extracts unique coordinate values
def get_unique_points(input_multilinestring):
    # Gets a nested list of coordinates
    coords_nested = [list(this_geom.coords) for this_geom in input_multilinestring.geoms]
    
    # Flattens the list
    coords = [item for sub_list in coords_nested for item in sub_list]

    # Gets only the unique coords
    coords_set = set()
    coords_final = []
    for i,coord in enumerate(coords):
        if coord not in coords_set:
            coords_set.add(coord)
            coords_final.append(coord)
    return coords_final

2 Answers 2

3

You might be able to just use shapely.ops.linemerge. I haven't used it much, and cannot speak to how it handles potential floating point differences where the geometries in the MultiLineString meet.

from shapely import wkt
from shapely.ops import linemerge


mls = wkt.loads('MULTILINESTRING ((120 105, 140 98),(160 103, 140 98),(100 100, 120 105))')

line_string = linemerge(mls.geoms)
# LINESTRING (100 100, 120 105, 140 98, 160 103)
1
  • Thank you, this is great! I need to look into some more details about linemerge, but it seems to be exactly what I'm looking for. Thanks again!!!
    – Felipe D.
    Commented Nov 8, 2021 at 14:06
0

I had an issue with the order of points in a geojson LineString. One possible solution is to reorder all the points by finding the closest point to each next point.

from math import radians, sin, cos, sqrt, atan2

def haversine_distance(lat1, lon1, lat2, lon2):
    R = 6371 # radius of Earth in km
    lat1, lon1, lat2, lon2 = map(radians, [lat1, lon1, lat2, lon2])
    dlat = lat2 - lat1
    dlon = lon2 - lon1
    a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
    c = 2 * atan2(sqrt(a), sqrt(1-a))
    distance = R * c
    return distance

def sort_points_by_distance(points, reference_point):
    sorted_points = []
    remaining_points = points.copy()
    current_point = reference_point
    while remaining_points:
        closest_point = min(remaining_points, key=lambda x: haversine_distance(current_point[0], current_point[1], x[0], x[1]))
        sorted_points.append(closest_point)
        remaining_points.remove(closest_point)
        current_point = closest_point
    return sorted_points

# Example usage
points = [(37.7749, -122.4194), (51.5074, -0.1278), (40.7128, -74.0060)]
reference_point = (37.7749, -122.4194)
sorted_points = sort_points_by_distance(points, reference_point)
print(sorted_points)

reference_point = start point

2
  • Hi there! Thanks for the suggestion! I think that "closest distance" might be a bit tricky, because we might have some segments that are curved, which might make it look that they are closer but not exactly the next one "in line", you know?
    – Felipe D.
    Commented Feb 22, 2023 at 19:34
  • 1
    You're right, this won't work for all cases. I usually run this for manual fixing of some broken routes. Commented Mar 1, 2023 at 6:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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