Search code examples
pythongispointprojectionshapely

Projection of a point to a line segment Python Shapely


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?


Solution

  • 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.