How to find the 2nd nearest point of a LineString in Shapely

Given a certain LineString and point p:

from shapely.ops import nearest_points
from shapely.geometry import Point

p = Point(51.21745162000732, 4.41871738126533)
linestring = LineString([(51.2176008, 4.4177154), (51.21758, 4.4178548), (51.2175729, 4.4179023), (51.21745162000732, 4.41871738126533)])

The nearest point to p is calculated by:

n_p = nearest_points(linestring, p)[0]

Conclusion it's the exact same point, which is normal since the exact same value is also in the linestring, but I need to know the nearest point, apart from the point itself. So how can I find the second nearest point?

Solution

In the general case, the simplest solution would be to construct a new geometric object from your LineString but without the nearest point, and then get the nearest point with this new geometry.:

from shapely.geometry import LineString, MultiPoint, Point
from shapely.ops import nearest_points

point = Point(51.21745162000732, 4.41871738126533)
line = LineString([(51.2176008, 4.4177154), (51.21758, 4.4178548), 
                   (51.2175729, 4.4179023), (51.21745162000732, 4.41871738126533)])

nearest_point = nearest_points(line, point)[0]
line_points_except_nearest = MultiPoint([point for point in linestring.coords 
                                         if point != (nearest_point.x, nearest_point.y)])
second_nearest = nearest_points(line_points_except_nearest, point)[0]

Alternatively, if you don't want to construct a new object because of, for example, memory constraints, you could run over all the points in the LineString with heapq.nsmallest:

import heapq

line_points = map(Point, line.coords)
nearest, second_nearest = heapq.nsmallest(2, line_points, key=point.distance)

In your specific case, when all the points are collinear, you can also calculate distances with the neighboring points of the nearest point:

index = list(line.coords).index((point.x, point.y))
if index == 0:
    second_nearest = Point(line.coords[1])
elif index == len(line.coords) - 1:
    second_nearest = Point(line.coords[-2])
else:
    second_nearest = min(Point(line.coords[index - 1]),
                         Point(line.coords[index + 1]),
                         key=point.distance)