Snap a point to the closest part of a line?

Question

I have a work order management system that has Python 2.7.

• In the system, it is only possible to use the libraries that are included in the standard Python 2.7 library (it's not possible to import other libraries).

---

The system has lines and points:

Line Vertex 1: 676561.00, 4860927.00
Line Vertex 2: 676557.00, 4860939.00

Point 100: 676551.00, 4860935.00
Point 200: 676558.00, 4860922.00

---

I want to snap the points to the nearest part of the line.

Snapping can be defined as:

Moving a point so that it coincides exactly with the closest part of a line.

---

There appear to be two different mathematical scenarios at play:

A. The closest part of a line is the position along the line where the point is perpendicular to the line.

B. Or, the closest part of the line is simply the closest vertex.

---

Is it possible to snap a point to the closest part of a line using Python 2.7 (and the standard 2.7 library)?

Answer 1

Another option:

# snap a point to a 2d line
# parameters: 
#   A,B: the endpoints of the line
#   C: the point we want to snap to the line AB
# all parameters must be a tuple/list of float numbers
def snap_to_line(A,B,C):
    Ax,Ay = A
    Bx,By = B
    Cx,Cy = C

    # special case: A,B are the same point: just return A
    eps = 0.0000001
    if abs(Ax-Bx) < eps and abs(Ay-By) < eps: return [Ax,Ay] 

    # any point X on the line can be represented by the equation
    #   X = A + t * (B-A)
    # so for point C we compute its perpendicular D on the line 
    # and the parameter t for D
    # if t is between 0..1 then D is on the line so the snap point is D
    # if t < 0 then the snap point is A
    # if t > 1 then the snap point is B
    #
    # explanation of the formula for distance from point to line:
    #    http://paulbourke.net/geometry/pointlineplane/
    #
    dx = Bx-Ax
    dy = By-Ay
    d2 = dx*dx + dy*dy
    t = ((Cx-Ax)*dx + (Cy-Ay)*dy) / d2
    if t <= 0: return A
    if t >= 1: return B
    return [dx*t + Ax, dy*t + Ay]


if __name__=="__main__":
    A,B = [676561.00, 4860927.00],[676557.00, 4860939.00]
    C = [676551.00, 4860935.00]
    print(snap_to_line(A,B,C))
    C = [676558.00, 4860922.00]
    print(snap_to_line(A,B,C))

Answer 2

You can calculate the line equation and then the distance between each point and the line.

For instance:

import collections
import math

Point = collections.namedtuple('Point', "x, y")


def distance(pt, a, b, c):
    # line eq: ax + by + c = 0
    return math.fabs(a * pt.x + b * pt.y + c) / math.sqrt(a**2 + b**2)


l1 = Point(676561.00, 4860927.00)
l2 = Point(676557.00, 4860939.00)

# line equation
a = l2.y - l1.y
b = l1.x - l2.x
c = l2.x * l1.y - l2.y * l1.x

assert a * l1.x + b * l1.y + c == 0
assert a * l2.x + b * l2.y + c == 0

p100 = Point(676551.00, 4860935.00)
p200 = Point(676558.00, 4860922.00)

print(distance(p100, a, b, c))
print(distance(p200, a, b, c))

You get:

6.957010852370434
4.427188724235731

---

Edit1: calculating the orthographic projection

What you want is the coordinates of the orthographic projection of p100 and p200 on the line (l1, l2).

You can calculate that as follow:

import collections
import math

Point = collections.namedtuple('Point', "x, y")


def snap(pt, pt1, pt2):
    # type: (Point, Point, Point) -> Point
    v = Point(pt2.x - pt1.x, pt2.y - pt1.y)
    dv = math.sqrt(v.x ** 2 + v.y ** 2)
    bh = ((pt.x - pt1.x) * v.x + (pt.y - pt1.y) * pt2.y) / dv
    h = Point(
        pt1.x + bh * v.x / dv,
        pt1.y + bh * v.y / dv
    )
    return h


l1 = Point(676561.00, 4860927.00)
l2 = Point(676557.00, 4860939.00)

p100 = Point(676551.00, 4860935.00)
p200 = Point(676558.00, 4860922.00)

s100 = snap(p100, l1, l2)
s200 = snap(p200, l1, l2)

print(s100)
print(p100)

You get:

Point(x=-295627.7999999998, y=7777493.4)
Point(x=676551.0, y=4860935.0)

You can check that the snapped points are on the line:

# line equation
a = l2.y - l1.y
b = l1.x - l2.x
c = l2.x * l1.y - l2.y * l1.x

assert math.fabs(a * s100.x + b * s100.y + c) < 1e-6
assert math.fabs(a * s200.x + b * s200.y + c) < 1e-6

---

Edit2: snap to the line segment

If you want to snap to a line segment, you need to check if the orthographic projection is inside the line segment or not.

• If the orthographic projection is inside the line segment: it is the solution,
• If it is near an extremity of the segment, this extremity is the solution.

You can do that as bellow:

def distance_pts(pt1, pt2):
    v = Point(pt2.x - pt1.x, pt2.y - pt1.y)
    dv = math.sqrt(v.x ** 2 + v.y ** 2)
    return dv


def snap(pt, pt1, pt2):
    # type: (Point, Point, Point) -> Point
    v = Point(pt2.x - pt1.x, pt2.y - pt1.y)
    dv = distance_pts(pt1, pt2)
    bh = ((pt.x - pt1.x) * v.x + (pt.y - pt1.y) * pt2.y) / dv
    h = Point(pt1.x + bh * v.x / dv, pt1.y + bh * v.y / dv)
    if 0 <= (pt1.x - h.x) / (pt2.x - h.y) < 1:
        # in the line segment
        return h
    elif distance_pts(h, pt1) < distance_pts(h, pt2):
        # near pt1
        return pt1
    else:
        # near pt2
        return pt2

The solutions for p100 and p200 are:

Point(x=676557.0, y=4860939.0)
Point(x=676551.0, y=4860935.0)

Answer 3

This is an extremely helpful resource.

import collections
import math

Line = collections.namedtuple('Line', 'x1 y1 x2 y2')
Point = collections.namedtuple('Point', 'x y')


def lineLength(line):
  dist = math.sqrt((line.x2 - line.x1)**2 + (line.y2 - line.y1)**2)
  return dist


## See http://paulbourke.net/geometry/pointlineplane/
## for helpful formulas

## Inputs
line = Line(0.0, 0.0, 100.0, 0.0)
point = Point(50.0, 1500)

## Calculations
print('Inputs:')
print('Line defined by: ({}, {}) and ({}, {})'.format(line.x1, line.y1, line.x2, line.y2))
print('Point "P": ({}, {})'.format(point.x, point.y))

len = lineLength(line)
if (len == 0):
  raise Exception('The points on input line must not be identical')

print('\nResults:')
print('Length of line (calculated): {}'.format(len))

u = ((point.x - line.x1) * (line.x2 - line.x1) + (point.y - line.y1) * (line.y2 - line.y1)) / (
    len**2)

# restrict to line boundary
if u > 1:
  u = 1
elif u < 0:
  u = 0

nearestPointOnLine = Point(line.x1 + u * (line.x2 - line.x1), line.y1 + u * (line.y2 - line.y1))
shortestLine = Line(nearestPointOnLine.x, nearestPointOnLine.y, point.x, point.y)

print('Nearest point "N" on line: ({}, {})'.format(nearestPointOnLine.x, nearestPointOnLine.y))
print('Length from "P" to "N": {}'.format(lineLength(shortestLine)))