Coordinates of the closest points of two geometries in Shapely
Solution 1
The GIS term you are describing is linear referencing, and Shapely has these methods.
# Length along line that is closest to the point
print(line.project(p))
# Now combine with interpolated point on line
p2 = line.interpolate(line.project(p))
print(p2) # POINT (5 7)
An alternative method is to use nearest_points
:
from shapely.ops import nearest_points
p2 = nearest_points(line, p)[0]
print(p2) # POINT (5 7)
which provides the same answer as the linear referencing technique does, but can determine the nearest pair of points from more complicated geometry inputs, like two polygons.
Solution 2
In case you have a single segment (e.g.: a line, as referring to the title) rather than a list of segments, here is what I did, and with a passing test case. Please consider that some users on this page are looking just for that from looking at the title, coming from a Google search.
Python code:
def sq_shortest_dist_to_point(self, other_point):
dx = self.b.x - self.a.x
dy = self.b.y - self.a.y
dr2 = float(dx ** 2 + dy ** 2)
lerp = ((other_point.x - self.a.x) * dx + (other_point.y - self.a.y) * dy) / dr2
if lerp < 0:
lerp = 0
elif lerp > 1:
lerp = 1
x = lerp * dx + self.a.x
y = lerp * dy + self.a.y
_dx = x - other_point.x
_dy = y - other_point.y
square_dist = _dx ** 2 + _dy ** 2
return square_dist
def shortest_dist_to_point(self, other_point):
return math.sqrt(self.sq_shortest_dist_to_point(other_point))
A test case:
def test_distance_to_other_point(self):
# Parametrize test with multiple cases:
segments_and_point_and_answer = [
[Segment(Point(1.0, 1.0), Point(1.0, 3.0)), Point(2.0, 4.0), math.sqrt(2.0)],
[Segment(Point(1.0, 1.0), Point(1.0, 3.0)), Point(2.0, 3.0), 1.0],
[Segment(Point(0.0, 0.0), Point(0.0, 3.0)), Point(1.0, 1.0), 1.0],
[Segment(Point(1.0, 1.0), Point(3.0, 3.0)), Point(2.0, 2.0), 0.0],
[Segment(Point(-1.0, -1.0), Point(3.0, 3.0)), Point(2.0, 2.0), 0.0],
[Segment(Point(1.0, 1.0), Point(1.0, 3.0)), Point(2.0, 3.0), 1.0],
[Segment(Point(1.0, 1.0), Point(1.0, 3.0)), Point(2.0, 4.0), math.sqrt(2.0)],
[Segment(Point(1.0, 1.0), Point(-3.0, -3.0)), Point(-3.0, -4.0), 1],
[Segment(Point(1.0, 1.0), Point(-3.0, -3.0)), Point(-4.0, -3.0), 1],
[Segment(Point(1.0, 1.0), Point(-3.0, -3.0)), Point(1, 2), 1],
[Segment(Point(1.0, 1.0), Point(-3.0, -3.0)), Point(2, 1), 1],
[Segment(Point(1.0, 1.0), Point(-3.0, -3.0)), Point(-3, -1), math.sqrt(2.0)],
[Segment(Point(1.0, 1.0), Point(-3.0, -3.0)), Point(-1, -3), math.sqrt(2.0)],
[Segment(Point(-1.0, -1.0), Point(3.0, 3.0)), Point(3, 1), math.sqrt(2.0)],
[Segment(Point(-1.0, -1.0), Point(3.0, 3.0)), Point(1, 3), math.sqrt(2.0)],
[Segment(Point(1.0, 1.0), Point(3.0, 3.0)), Point(3, 1), math.sqrt(2.0)],
[Segment(Point(1.0, 1.0), Point(3.0, 3.0)), Point(1, 3), math.sqrt(2.0)]
]
for i, (segment, point, answer) in enumerate(segments_and_point_and_answer):
result = segment.shortest_dist_to_point(point)
self.assertAlmostEqual(result, answer, delta=0.001, msg=str((i, segment, point, answer)))
Note: I assume this function is inside a Segment
class.
In case your line is infinite, don't limit the lerp
from 0 to 1 only, but still at least provide two distinct a
and b
points.
Asif Rehan
Updated on July 14, 2022Comments
-
Asif Rehan almost 2 years
There is a polyline with a list of coordinates of the vertices = [(x1,y1), (x2,y2), (x3,y3),...] and a point(x,y). In Shapely,
geometry1.distance(geometry2)
returns the shortest distance between the two geometries.>>> from shapely.geometry import LineString, Point >>> line = LineString([(0, 0), (5, 7), (12, 6)]) # geometry2 >>> list(line.coords) [(0.0, 0.0), (5.0, 7.0), (12.0, 6.0)] >>> p = Point(4,8) # geometry1 >>> list(p.coords) [(4.0, 8.0)] >>> p.distance(line) 1.4142135623730951
But I also need to find the coordinate of the point on the line that is closest to the point(x,y). In the example above, this is the coordinate of the point on the
LineString
object that is 1.4142135623730951 unit distant fromPoint(4,8)
. The methoddistance()
should have the coordinates when calculating the distance. Is there any way to get it returned from this method? -
Asif Rehan almost 10 yearsquick note for anyone who might be following the post:
line.project(p)
measures the projection point along the line from the starting node-coordinate given in the argument of LineString() method. Butnp = line.interpolate(line.project(p))
gives the POINT (5 7) as in global coordinates, not from the start of the line. Here the line started from (0,0). So it might be confusing. -
Dmitri Chubarov over 4 yearsThis is a perfect solution. How about using the "precise" lerp to compute
x
andy
as inx = (1 - lerp) * self.a.x + lerp * self.b.x
? See for instance stackoverflow.com/a/23716956/1328439 for more details and a comparison. -
user14241 almost 4 yearsBe careful declaring a variable 'np' as this maybe conflict with a numpy import.
-
SchemeSonic over 2 yearsHe is asking for the coordinate of the point on the line if he refers to a vertex on the line, Shapely documents says that: Note that the nearest points may not be existing vertices in the geometries