Python: Closest Point to a line

Question

I have the following question. I have a box full of coordinates and three points, which build up a line. Now i want to calculate the shortest distance of all box coordinates to that line. I have three methods to do that and the vtk and numpy version always have the same result but not the distance method of shapely. But i need the shapely version, because i want to measure the closest distance from a point to the hwole line and not to the seperate line segments. Here is an example code so far. So the problem is the "pdist":

from shapely.geometry import LineString, Point
import vtk, numpy as np
import itertools

import math

from numpy.linalg import norm

x1=np.arange(4,21)
y1=np.arange(4,21)
z1=np.arange(-7,6)

linepoints = np.array([[1,10,0],[10,10,0],[15,15,0]])


for i in itertools.product(x1,y1,z1):

    for m in range(len(linepoints)-1):

        line3 = LineString([linepoints[m],linepoints[m+1]])

        p = Point(i)

        d = norm(np.cross(linepoints[m]-linepoints[m+1], linepoints[m]-i))/norm(linepoints[m+1]-linepoints[m])

        dist=math.sqrt(vtk.vtkLine().DistanceToLine(i,linepoints[m],linepoints[m+1]))

        pdist = p.distance(line3)

        print(d,dist,pdist)

Answer 1

The problem is that with cross-product you are calculating orthogonal distance to the line spanned by the segment defined by points linepoints[m] and linepoints[m+1]. On the other hand, Shapely calculates distance to the segment, i.e., it returns the distance either to the orthogonal projection or to one of the boundary points should the orthogonal projection fall "outside" of the segment.

To get consistent results, you could calculate the orthogonal projection yourself and then invoke the Shapely distance method:

import numpy as np
from shapely.geometry import Point, LineString


A = np.array([1,0])
B = np.array([3,0])
C = np.array([0,1])


l = LineString([A, B])
p = Point(C)


d = np.linalg.norm(np.cross(B - A, C - A))/np.linalg.norm(B - A)

n = B - A
v = C - A

z = A + n*(np.dot(v, n)/np.dot(n, n))

print(l.distance(p), d, Point(z).distance(p))
#1.4142135623730951 1.0 1.0

However, note that Shapely effectively ignores the z-coordinate. Thus for example:

import numpy as np
from shapely.geometry import Point, LineString

A = np.array([1,0,1])
B = np.array([0,0,0])

print(Point([1,0,1]).distance(Point([0,0,0])))

return as distance merely 1.

EDIT: based on your comment, here would be a version which calculates the distance (for arbitrary number of dimensions) to the segment:

from shapely.geometry import LineString, Point
import numpy as np
import itertools

import math

from numpy.linalg import norm

x1=np.arange(4,21)
y1=np.arange(4,21)
z1=np.arange(-7,6)

linepoints = np.array([[1,10,0],[10,10,0],[15,15,0]])

def dist(A, B, C):
    """Calculate the distance of point C to line segment spanned by points A, B.

    """

    a = np.asarray(A)
    b = np.asarray(B)
    c = np.asarray(C)

    #project c onto line spanned by a,b but consider the end points
    #should the projection fall "outside" of the segment    
    n, v = b - a, c - a

    #the projection q of c onto the infinite line defined by points a,b
    #can be parametrized as q = a + t*(b - a). In terms of dot-products,
    #the coefficient t is (c - a).(b - a)/( (b-a).(b-a) ). If we want
    #to restrict the "projected" point to belong to the finite segment
    #connecting points a and b, it's sufficient to "clip" it into
    #interval [0,1] - 0 corresponds to a, 1 corresponds to b.

    t = max(0, min(np.dot(v, n)/np.dot(n, n), 1))
    return np.linalg.norm(c - (a + t*n)) #or np.linalg.norm(v - t*n)


for coords in itertools.product(x1,y1,z1):
    for m in range(len(linepoints)-1):

        line3 = LineString([linepoints[m],linepoints[m+1]])
        d = dist(linepoints[m], linepoints[m+1], coords)
        print(coords, d)