4

I was trying to create river cross-section profiles based on the point terrestical measurements. When trying to create a Shapely LineString from a Series of points with the common id, I realized that the order of given points really matters as the LineString would just connect given points 'indexwise' (connect points in the list-given order). The below code illustrates the default behaviour:

from shapely.geometry import Point, LineString
import geopandas as gpd
import numpy as np
import matplotlib.pyplot as plt

# Generate random points
x=np.random.randint(0,100,10)
y=np.random.randint(0,50,10)
data = zip(x,y)

# Create Point and default LineString GeoSeries
gdf_point = gpd.GeoSeries([Point(j,k) for j,k in data])
gdf_line = gpd.GeoSeries(LineString(zip(x,y)))

# plot the points and "default" LineString
ax = gdf_line.plot(color='red')
gdf_point.plot(marker='*', color='green', markersize=5,ax=ax)

That would produce the image:

Default LineString

Question: Is there any built-in method within Shapely that would automatically create the most logical (a.k.a.: the shortest, the least complicated, the least criss-cross,...) line through the given list of random 2D points?

Below can you find the desired line (green) compared to the default (red).

Desired LineString

Marjan Moderc
  • 2,747
  • 23
  • 44
  • 1
    Assuming you don't know the order or neighbours ahead of time, you could try building a graph that connects every node to every other node, then search for 'simple paths' and select the paths with the same number of steps as the number of nodes, then select the shortest of these? This would require something like the [`all_simple_paths`](https://networkx.readthedocs.io/en/stable/reference/generated/networkx.algorithms.simple_paths.all_simple_paths.html#networkx.algorithms.simple_paths.all_simple_paths) in networkX. – songololo Jan 06 '17 at 11:16
  • Wow, looks promising! Will take a look into this. – Marjan Moderc Jan 06 '17 at 11:22
  • small correction: the path length would be nodes - 1 – songololo Jan 06 '17 at 11:32

3 Answers3

3

Here is what solved my cross-section LineString simplification problem. However, my solution doesn't correctly address computationally more complex task of finding the ultimately shortest path through the given points. As the commenters suggested, there are many libraries and scripts available to solve that particulal problem, but in case anyone want to keep it simple, you can use what did the trick for me. Feel free to use and comment!

def simplify_LineString(linestring):

    '''
    Function reorders LineString vertices in a way that they each vertix is followed by the nearest remaining vertix.
    Caution: This doesn't calculate the shortest possible path (travelling postman problem!) This function performs badly
    on very random points since it doesn't see the bigger picture.
    It is tested only with the positive cartesic coordinates. Feel free to upgrade and share a better function!

    Input must be Shapely LineString and function returns Shapely Linestring.

    '''

    from shapely.geometry import Point, LineString
    import math

    if not isinstance(linestring,LineString):
        raise IOError("Argument must be a LineString object!")

    #create a point lit
    points_list = list(linestring.coords)

    ####
    # DECIDE WHICH POINT TO START WITH - THE WESTMOST OR SOUTHMOST? (IT DEPENDS ON GENERAL DIRECTION OF ALL POINTS)
    ####
    points_we = sorted(points_list, key=lambda x: x[0])
    points_sn = sorted(points_list, key=lambda x: x[1])

    # calculate the the azimuth of general diretction
    westmost_point = points_we[0]
    eastmost_point = points_we[-1]

    deltay = eastmost_point[1] - westmost_point[1]
    deltax = eastmost_point[0] - westmost_point[0]

    alfa = math.degrees(math.atan2(deltay, deltax))
    azimut = (90 - alfa) % 360

    if (azimut > 45 and azimut < 135):
        #General direction is west-east
        points_list = points_we
    else:
        #general direction is south-north
        points_list = points_sn

    ####
    # ITERATIVELY FIND THE NEAREST VERTIX FOR THE EACH REMAINING VERTEX
    ####

    # Create a new, ordered points list, starting with the east or southmost point.
    ordered_points_list = points_list[:1]

    for iteration in range(0, len(points_list[1:])):

        current_point = ordered_points_list[-1]  # current point that we are looking the nearest neighour to
        possible_candidates = [i for i in points_list if i not in ordered_points_list]  # remaining (not yet sortet) points

        distance = 10000000000000000000000
        best_candidate = None
        for candidate in possible_candidates:
            current_distance = Point(current_point).distance(Point(candidate))
            if current_distance < distance:
                best_candidate = candidate
                distance = current_distance

        ordered_points_list.append(best_candidate)

    return LineString(ordered_points_list)
Marjan Moderc
  • 2,747
  • 23
  • 44
  • Adding an `if ordered_point_list[-1] is None: ordered_point_list.pop()` at the end should allow this to handle LinearRings as well, via `return type(linestring)(ordered_points_list)`. – Dominik Stańczak Dec 04 '19 at 06:57
1

There is no built in function, but shapely has a distance function.

You could easily iterate over the points and calculate the shortest distance between them and construct the 'shortest' path.

There are some examples in the offical github repo.

Maurice Meyer
  • 17,279
  • 4
  • 30
  • 47
  • I have already made a similar function myself, but it only starts with the westmost or southmost point and iteratively finds the closest neighbour for every remaining point. But this in fact doesn't mean the shortest path all together. I will post more if I come up with a decent solution. Thanks anyways! – Marjan Moderc Jan 06 '17 at 10:37
0

Google's OR-Tools offer a nice and efficient way for solving the Travelling Salesman Problem: https://developers.google.com/optimization/routing/tsp.

Following the tutorial on their website would give you a solution from this (based on your example code):

enter image description here

to this:

enter image description here

Markus
  • 165
  • 11