Projection of a point to a line segment Python Shapely

sofaruntitled picture sofaruntitled · Mar 2, 2018 · Viewed 7.7k times · Source

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.

enter image description here

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?

Answer

ewcz picture ewcz · Mar 2, 2018

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.