# How to calculate the nearest point on a line to a given point

Iβve a line string that has long/lat coordinates as follows:

[
[
146.54294815408372,
-36.04985956921628
],
[
146.5457684472937,
-36.04924988190973
],
[
146.54706548053542,
-36.04879871026288
],
[
146.54919201177802,
-36.047176909745176
],
[
146.55074543531066,
-36.04585993436855
]
]


And a point input like this:

[
146.54376257030356,
-36.04616479173944
]


So, if we given these two input to a function, I need to calculate the nearest point on that line. As shown in the following picture, Black dots shows the line input and orange dot shows the expected output from the function.

Calculate the Euclidean distance - Wikipedia

1. Parse your input to get a list of line segments (\vec{s}_i, \vec{t}_i).

2. Project your point \vec{p} onto each line segment to get \vec{q}_i:

Let \vec{v}_i = \vec{t}_i - \vec{s}_i and \vec{p}_i = \vec{p} - \vec{s}_i. Then \vec{s}_i + t \vec{v}_i, where t = \frac{\vec{v}^T \vec{p}_i}{\vec{v}_i^T \vec{v}_i}, is the orthogonal projection onto line (\vec{s}_i, \vec{t}_i). You are interested in line segments, so we need to clamp the parameter t. So, \vec{q}_i = \vec{s}_i + \min(1, \max(0, t)) \vec{v}_i.

1. Use \vec{q}_i with the minimum distance d_i from \vec{p}:

You can use the squared euclidean distance, which is simply a dot product: d_i = (\vec{p} - \vec{q}_i)^T (\vec{p} - \vec{q}_i).

In julia, the second point can be implemented like this

using LinearAlgebra
# given a point p and a line segment starting at point s and ending at point t,
# find a point on that line that is closest to p
function nearest_point_on_line(p, s, t)
@assert ndim(s) == ndim(t) == ndim(p)
v = (t - s)
t = clamp((v β (p - s)) / (v β v), 0, 1)
return s + t * v
end

6 Likes

Just to add a solution using GEOS nearestPoints, which works between a larger set of geometries and is battle-tested:

import LibGEOS

linestring = LibGEOS.LineString([
[146.54294815408372, -36.04985956921628],
[146.5457684472937, -36.04924988190973],
[146.54706548053542, -36.04879871026288],
[146.54919201177802, -36.047176909745176],
[146.55074543531066, -36.04585993436855],
])
point = LibGEOS.Point(146.54376257030356, -36.04616479173944)

on_linestring, on_point = LibGEOS.nearestPoints(linestring, point)


We can check the result visually:

using CairoMakie, GeoInterfaceMakie
# enable plotting LibGEOS geometries
GeoInterfaceMakie.@enable LibGEOS.AbstractGeometry

fig, ax, p = plot(linestring; axis=(aspect=DataAspect(),))
plot!(ax, point)
plot!(ax, on_linestring)
fig


5 Likes

Amazing @visr , Youβre a lifesaver. This is exactly what I needed.

Just a small question along the same, Since these line strings are on the earthβs surface, how would it affect the accuracy of the returning on_linestring value?

1 Like

The problem of the GEOS solution is β¦ itβs not correct because GEOS is a Cartesian library. When working with geographical coordinates you want the point to be along a great circle. In this case the difference is small, but it is no zero. GMT computes it though I find some issues that need to be ironed up. Namely, to compute this one need to save the lines segments in a file first

using GMT

gmtwrite("line.dat", [146.54294815408372 -36.04985956921628; 146.5457684472937 -36.04924988190973; 146.54706548053542 -36.04879871026288; 146.54919201177802 -36.047176909745176; 146.55074543531066 -36.04585993436855])

pt2 = mapproject([146.54376257030356 -36.04616479173944], dist2line="line.dat")

1Γ5 GMTdataset{Float64, 2}
Row β   col.1     col.2    col.3    col.4     col.5
ββββββΌβββββββββββββββββββββββββββββββββββββββββββββββ
1 β 146.544  -36.0462  376.851  146.545  -36.0495


The the result is in columns 4 & 5 (column 3 is the distance in meters).
Also PrettyTables is not printing enough decimals in this case, so we need to

julia> pt2[4]
146.544839684871

julia> pt2[5]
-36.049450668819865


For reference, the GEOS solution is

julia> on_linestring
POINT (146.54448928065796 -36.04952641048415)


If the line segment was longer, the difference between Cartesian and Geographic would be larger.

6 Likes

Indeed itβs good to stress that GEOS assumes a plane. I added the answer by @joa-quim to my plot as a red dot to show the difference.

You could also project your points first using e.g. Proj.jl. In time it would also be nice to be able to use S2 Geometry, which assumes a sphere.

4 Likes

Oh nice @joa-quim , This is exactly what I needed. Appreciate the help mate. Just wondering are you by any chance familiar with this error ? Iβm using Mac with M1 processor.

Failed to precompile GMT [5752ebe1-31b9-557e-87aa-f909b540aa54] to /Users/cbandara/.julia/compiled/v1.6/GMT/jl_rcWYt7.
ERROR: LoadError: SystemError: opening file β/Users/cbandara/.julia/packages/GMT/1kBBX/deps/deps.jlβ: No such file or directory

When I tried to install this libaray, it keeps giving me this error. Iβm using Mac with M1 processor.

1 Like

Thanks @visr, Appreciate the help on this.

Btw, by any chance are you familiar with this error ? I got this when I try to add the library ?

Failed to precompile GMT [5752ebe1-31b9-557e-87aa-f909b540aa54] to /Users/cbandara/.julia/compiled/v1.6/GMT/jl_rcWYt7.
ERROR: LoadError: SystemError: opening file β/Users/cbandara/.julia/packages/GMT/1kBBX/deps/deps.jlβ: No such file or directory

Hmm, that file (deps.jl) is generated by the build process. Donβt know why it was not in your case.
Try to run import Pkg; Pkg.build("GMT"). If it doesnβt solve it better to move this into a GMT.jl issue

There is the issue of multiple results with the same (or almost the same distance). Perhaps an eps parameter will help to determine closest points, and a vector/set returned. Or a findfirst logic, or return the closest point and the gap to the next closest point.
The simplest example to think of is the middle of an equilatral triangle.

1 Like

You mean using the mapproject function from GMT ?

I mean, when the input is something like this:

The result could be any of the red points. And it isnβt clear how to choose which one, or to return them all.

Tried a case like that (a rectangle to make it simple) in GMT and it returned only one point. Presumably the first it found.

To output multiple nearest points within a specified tolerance:

Code
# Assumes that the input line nodes are spatially ordered

distance(P1, P2) = hypot(P1[1] - P2[1], P1[2] - P2[2])

dot(x,y) = sum(a*b for (a,b) in zip(x,y))

# Given a line segment starting at point S[1] and ending at point S[2],
# find a point on that line that is closest to P0
function nearest_point_on_segment(S, P0)
V = S[2] - S[1]
t = clamp(dot(V, P0 - S[1]) / dot(V, V), 0., 1.)
return S[1] + t * V
end

# Input data
P = [ [0.0, 0.0], [1.0, 0.0], [0.5, sqrt(3)/2], [0.0, 0.0] ]
P0 =  [0.5, sqrt(3)/6]

# Find the closest points on each of the segments
S = [ [P[i], P[i+1]] for i in 1:length(P)-1]
Pn = nearest_point_on_segment.(S, Ref(P0))
d = distance.(Pn, Ref(P0))
dmin = minimum(d)
Ο΅ = 1e-3 * dmin   # tolerance on distance
ix = findall(x -> abs(x - dmin) < Ο΅, d)

# Plot results
using Plots
p = plot(first.(P), last.(P), m=:o, ms=3, ratio=1, lims=(-0.1,1.1))
scatter!((P0[1], P0[2]), c=:red, label="P0")
for i in ix; scatter!((Pn[i][1], Pn[i][2]), c=:red, m=:x, label="P\$i"); end
p

1 Like

@Dan Ahh, Yes that makes sense. This is the issue with LibGEOS library, based on the answer provided by @visr , does that mean on_linestring can have multiple values when the line string is almost a triangle?

@Chamika_Kasun if you ever feel the need to modify these algorithms or use a different combination of functions without reaching out to external libraries not written in Julia, consider reading the source code of Meshes.jl.

Regarding distance calculations, we plan to add more methods here after we finish the implementation of a related set of features. Regarding nearest neighbors, take a look at the examples here.