Python: point on a line closest to third point

Question

I have a line/vector between two XY points (p1 and p2) and a third XY point (p3) that is outside the line. According to this post I know how to get the distance of that point to the line. But what I'm actually looking for is a point (p4) on that line that is in a minimum distance (d) to the third point (p3). I found this post, but I feel it's not the correct solution. Maybe there's something included in Numpy or Python?

According to @allo I tried the following. You can download my code as Python file or Jupyter Notebook (both Python3).

points = [[1, 1], [3, 1], [2.5, 2], [2.5, 1]]
import matplotlib.pyplot as plt
%matplotlib inline

fig, ax = plt.subplots()
fig.set_size_inches(6,6)

x, y = zip(*points[:2])
l1, = ax.plot(x,y, color='blue')
scatter1 = ax.scatter(x=x,y=y, color='blue', marker='x', s=80, alpha=1.0)

x, y = zip(*points[2:])
l2, = ax.plot(x,y, color='red')
scatter2 = ax.scatter(x=x,y=y, color='red', marker='x', s=80, alpha=1.0)

p1 = Vector2D(*points[0])
p2 = Vector2D(*points[1])
p3 = Vector2D(*points[2])

p1p2 = p2.sub_vector(p1)
p1p3 = p3.sub_vector(p1)

angle_p1p2_p1p3 = p1p2.get_angle_radians(p1p3)
length_p1p3 = p1p3.get_length()
length_p1p2 = p1p2.get_length()

p4 = p1.add_vector(p1p2.multiply(p1p3.get_length()/p1p2.get_length()).multiply(math.cos(p1p2.get_angle_radians(p1p3))))

#p4 = p1 + p1p2 * length(p1p3)/length(p1p2)*cos(angle(p1p2, p1p3))

p4 = p1.add_vector(p1p2.multiply(length_p1p3/length_p1p2*math.cos(angle_p1p2_p1p3)))
p4

Which results in p4 = (1.8062257748298551, 1.0) but should obviously be (2.5, 1.0).

Answer 1

Analytical Geometry

Let's start with the assigned line, we define the line in terms of two points on it (x1, y1) and (x2, y2).

With dx = x2-x1 and dy = y2-y1 we can formally write every point on the line as (x12, y12) = (x1, y1) + a*(dx, dy) where a is a real number.

Using an analogous notation a point on the line passing in (x3, y3) and perpendicular to the assigned one is (x34, y34) = (x3, y3) + b*(-dy, +dx).

To find the intersection we have to impose (x12, y12) = (x34, y34) or (x1, y1) + a*(dx, dy) = (x3, y3) + b*(-dy, +dx).

Writing separately the equations for x and y

y1 + a dy - y3 - b dx = 0
x1 + a dx + b dy - x3 = 0

it is a linear system in a and b whose solution is

a = (dy y3 - dy y1 + dx x3 - dx x1) / (dy^2 + dx^2)
b = (dy x3 - dy x1 - dx y3 + dx y1) / (dy^2 + dx^2)

The coordinates of the closest point to (x3, y3) lying on the line are (x1+a*dx, y1+a*dy) — you need to compute only the coefficient a.

Numerically speaking, the determinant of the linear system is dx**2+dy**2 so you have problems only when the two initial points are extremely close to each other with respect to their distance w/r to the third point.

Python Implementation

We use a 2-uple of floats to represent a 2D point and we define a function whose arguments are 3 2-uples representing the points that define the line (p1, p2) and the point (p3) that determines the position of p4 on said line.

In [16]: def p4(p1, p2, p3):
    ...:     x1, y1 = p1
    ...:     x2, y2 = p2
    ...:     x3, y3 = p3
    ...:     dx, dy = x2-x1, y2-y1
    ...:     det = dx*dx + dy*dy
    ...:     a = (dy*(y3-y1)+dx*(x3-x1))/det
    ...:     return x1+a*dx, y1+a*dy

To test the implementation I am using the three points used by the OP to demonstrate their issues with this problem:

In [17]: p4((1.0, 1.0), (3.0, 1.0), (2.5, 2))
Out[17]: (2.5, 1.0)

It seems that the result of p4(...) coincides with the OP expectation.

Answer 2

Shapely's distance() function returns the minimum distance:

>>> from shapely.geometry import LineString as shLs
>>> from shapely.geometry import Point as shPt
>>> l = shLs([ (1,1), (3,1)])
>>> p = shPt(2,2)
>>> dist = p.distance(l)
1.0
>>> l.interpolate(dist).wkt
'POINT (2 1)'