I have a LineString defined by two points, so essentially a straight line segment, and I wanted to project a point on to it. I am aware of .project
and .interpolate
. However when the point is "outside" the segment, I don't want the closest point on the segment, but I want to extend the segment and draw a line going through the point and is orthogonal to the (extended) line segment. I want the coordinate of the projection.
For example if the point is "within" the segment
from shapely.geometry import Point
from shapely.geometry import LineString
point = Point(0.2, 0.5)
dist = LineString([(0, 1), (1, 1)]).project(point)
list(LineString([(0, 1), (1, 1)]).interpolate(dist).coords)
Anyone knows what to do when the point is outside of the segment?
It will be probably easiest to do this manually. If you denote the angle x - u - v
as alpha
, then
cos(alpha) = (v - u).(x - u) / (|x - u|*|v - u|)
where .
denotes the dot product, and | |
represents the Euclidean norm.
The distance d
of P(x)
from u
is therefore:
d = cos(alpha)*|x - u| = (v - u).(x - u) / |v - u|
Having calculated d
, the projected point P(x)
is then easily obtained as:
P(x) = u + d*(v - u)/|v - u|
The implementation:
import numpy as np
from shapely.geometry import Point
from shapely.geometry import LineString
point = Point(0.2, 0.5)
line = LineString([(0, 1), (1, 1)])
x = np.array(point.coords[0])
u = np.array(line.coords[0])
v = np.array(line.coords[len(line.coords)-1])
n = v - u
n /= np.linalg.norm(n, 2)
P = u + n*np.dot(x - u, n)
print(P) #0.2 1.