Skip to content

Fix @turf/nearest-point-on-line endpoint selection and degenerate input cases #2940

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 5 commits into from
Nov 11, 2025
+189 −93
Merged

Fix @turf/nearest-point-on-line endpoint selection and degenerate input cases #2940

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
1c5f047
@turf/nearest-point-on-line code simplifications
bratter Oct 30, 2025
a80a828
Fix @turf nearest-point-on-line incorrect endpoints
bratter Nov 3, 2025
7511687
Fix: @turf/nearest-point-on-line failing test
bratter Nov 10, 2025
fcc30b8
Merge branch 'master' into fix-nearest-point-on-line-incorrect-endpoi…
bratter Nov 11, 2025
cf33d55
Fix: @turf/nearest-point-on-line faililng test
bratter Nov 11, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
158 changes: 85 additions & 73 deletions packages/turf-nearest-point-on-line/index.ts
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -75,12 +75,10 @@ function nearestPointOnLine<G extends LineString | MultiLineString>(
for (let i = 0; i < coords.length - 1; i++) {
//start - start of current line section
const start: Feature<Point, { dist: number }> = point(coords[i]);
start.properties.dist = distance(pt, start, options);
const startPos = getCoord(start);

//stop - end of current line section
const stop: Feature<Point, { dist: number }> = point(coords[i + 1]);
stop.properties.dist = distance(pt, stop, options);
const stopPos = getCoord(stop);

// sectionLength
Expand All @@ -89,40 +87,28 @@ function nearestPointOnLine<G extends LineString | MultiLineString>(
let wasEnd: boolean;

// Short circuit if snap point is start or end position of the line
// segment or if start is equal to stop position.
if (startPos[0] === ptPos[0] && startPos[1] === ptPos[1]) {
[intersectPos, , wasEnd] = [startPos, undefined, false];
} else if (stopPos[0] === ptPos[0] && stopPos[1] === ptPos[1]) {
[intersectPos, , wasEnd] = [stopPos, undefined, true];
} else if (startPos[0] === stopPos[0] && startPos[1] === stopPos[1]) {
[intersectPos, , wasEnd] = [stopPos, undefined, true];
// Test the end position first for consistency in case they are
// coincident
if (stopPos[0] === ptPos[0] && stopPos[1] === ptPos[1]) {
[intersectPos, wasEnd] = [stopPos, true];
} else if (startPos[0] === ptPos[0] && startPos[1] === ptPos[1]) {
[intersectPos, wasEnd] = [startPos, false];
} else {
// Otherwise, find the nearest point the hard way.
[intersectPos, , wasEnd] = nearestPointOnSegment(
start.geometry.coordinates,
stop.geometry.coordinates,
getCoord(pt)
[intersectPos, wasEnd] = nearestPointOnSegment(
startPos,
stopPos,
ptPos
);
}
let intersectPt:
| Feature<
Point,
{ dist: number; multiFeatureIndex: number; location: number }
>
| undefined;

if (intersectPos) {
intersectPt = point(intersectPos, {
dist: distance(pt, intersectPos, options),
multiFeatureIndex: multiFeatureIndex,
location: length + distance(start, intersectPos, options),
});
}

if (
intersectPt &&
intersectPt.properties.dist < closestPt.properties.dist
) {
const intersectPt = point(intersectPos, {
dist: distance(pt, intersectPos, options),
multiFeatureIndex: multiFeatureIndex,
location: length + distance(start, intersectPos, options),
});

if (intersectPt.properties.dist < closestPt.properties.dist) {
closestPt = {
...intersectPt,
properties: {
Expand Down Expand Up @@ -164,10 +150,15 @@ function cross(v1: Vector, v2: Vector): Vector {
return [v1y * v2z - v1z * v2y, v1z * v2x - v1x * v2z, v1x * v2y - v1y * v2x];
}

function magnitude(v: Vector) {
function magnitude(v: Vector): number {
return Math.sqrt(Math.pow(v[0], 2) + Math.pow(v[1], 2) + Math.pow(v[2], 2));
}

function normalize(v: Vector): Vector {
const mag = magnitude(v);
return [v[0] / mag, v[1] / mag, v[2] / mag];
}

function angle(v1: Vector, v2: Vector): number {
const theta = dot(v1, v2) / (magnitude(v1) * magnitude(v2));
return Math.acos(Math.min(Math.max(theta, -1), 1));
Expand All @@ -185,7 +176,10 @@ function lngLatToVector(a: Position): Vector {

function vectorToLngLat(v: Vector): Position {
const [x, y, z] = v;
const lat = radiansToDegrees(Math.asin(z));
// Clamp the z-value to ensure that is inside the [-1, 1] domain as required
// by asin. Note that
Copy link
Collaborator

@mfedderly mfedderly Nov 8, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hanging Note that in the comment

Copy link
Contributor Author

@bratter bratter Nov 8, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah whoops. I'll fix that when I get a chance and also look at the point-on-polygon tests that failed. These were succeeding when I tried earlier, but are worthwhile checking into due to the floating point funniness.

Copy link
Contributor Author

@bratter bratter Nov 10, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I was able to confirm that the failed tests in turf-point-to-polygon-distance are due to a single bit difference in the latitude of the nearest point between the old version and the new version. See below for an example from the first failing test:

Latitude from the nearest failing point
Old: 52.47479359036524 => 0 10000000100 1010 00111100 11000110 00001001 01001111 01111100 00001100
New: 52.47479359036523 => 0 10000000100 1010 00111100 11000110 00001001 01001111 01111100 00001011

Then as the result is in meters, the single bit difference creeps up. Yay for floating point :-).

Given that a single bit floating point difference is inside even the boundaries of trig functions in different JS engines I see no issue in just changing the expected distance, although it does beg a larger questions as to whether those tests should test with a tolerance like some others do... problem for another time.

I will amend and update the PR.

Copy link
Collaborator

@mfedderly mfedderly Nov 10, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for the detailed readout from this investigation!

const zClamp = Math.min(Math.max(z, -1), 1);
const lat = radiansToDegrees(Math.asin(zClamp));
const lng = radiansToDegrees(Math.atan2(y, x));

return [lng, lat];
Expand All @@ -195,7 +189,7 @@ function nearestPointOnSegment(
posA: Position, // start point of segment to measure to
posB: Position, // end point of segment to measure to
posC: Position // point to measure from
): [Position, boolean, boolean] {
): [Position, boolean] {
// Based heavily on this article on finding cross track distance to an arc:
// https://gis.stackexchange.com/questions/209540/projecting-cross-track-distance-on-great-circle

Expand All @@ -206,62 +200,80 @@ function nearestPointOnSegment(
const B = lngLatToVector(posB); // ... to posB
const C = lngLatToVector(posC); // ... to posC

// Components of target point.
const [Cx, Cy, Cz] = C;
// The axis (normal vector) of the great circle plane containing the line segment
const segmentAxis = cross(A, B);

// Two degenerate cases exist for the segment axis cross product. The first is
// when vectors are aligned (within the bounds of floating point tolerance).
// The second is where vectors are antipodal (again within the bounds of
// tolerance. Both cases produce a [0, 0, 0] cross product which invalidates
// the rest of the algorithm, but each case must be handled separately:
// - The aligned case indicates coincidence of A and B. therefore this can be
// an early return assuming the closest point is the end (for consistency).
// - The antipodal case is truly degenerate - an infinte number of great
// circles are possible and one will always pass through C. We assume this
// is undefined behavior and therefore throw. Callers can catch this and
// return 0 if they wish.
if (segmentAxis[0] === 0 && segmentAxis[1] === 0 && segmentAxis[2] === 0) {
if (dot(A, B) > 0) {
return [[...posB], true];
} else {
throw new Error(
Copy link

@HarelM HarelM Nov 2, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure I fully understand this part.
the bottom line from my perspective is that I don't expect nearest-point-on-line to ever throw.
So if this is caused by ambiguity (if I understand what you wrote above) then I would expect this to return the first one, even if it's not deterministic, I would still prefer that over a method that throws...
I might have completely misunderstood what you wrote, but I also couldn't see a test to over this code flow, so do let me know what I'm missing.

Copy link
Contributor Author

@bratter bratter Nov 2, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, I totally get what you are saying and could go either way. In case its helpful, I'll lay out all the conflicting thought process:

This case should only be hit when the input line segment is made up of two points that are antipodal. Theoretically, for instance, [0, 0] and [180, 0], but practically doesn't quite work like that because of floating point. This is mathematically undefined because the is an infinite number of valid great circles that pass through antipodal points.

So the two options are:

  1. Throw for undefined behavior:
    1. This give stricter mathematical correctness, which is good in some cases but may not be necessary when the focus is pragmatism.
    2. Gives the caller an opportunity to handle in the best way for the application - for instance us picking the right answer may not be the best in animation scenarios where it might cause snapping etc.
    3. Downside is of course that it throws and throwing is always less than ideal, people won't handle it, and might crash the application if not caught.
  2. Assume the point lies on the line:
    1. More permissively, we can simply choose one of the possible great circles. The only one that really makes sense is the one that also passes through the test point, making the point exactly on the line. I don't think there are any other options that would possibly make sense.
    2. This avoids any exceptions in nearest-point-on-line, which as you point out seems like a reasonable expectation.
    3. Determining that the point on the line in this case is not entirely logically unsound, and would likely be a reasonable result in almost all circumstances, especially given this antipodes case is unlikely to occur in practice.
    4. Antipodal points in this kind of method are kind of crazy anyway, and I doubt are seen at all in practice - and if it mattered, you'd think the geometry would include an intermediate point to disambiguate.

So I can see advantages each way. Wanted to drop the more correct version in for the first pass, but definitely open to changing.

Copy link

@HarelM HarelM Nov 2, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If someone needs to know that this might be problematic you could export a method called isAntiPodal to allow testing it if someone is concerned about the results of this method.
This avoid throwing and allows to decide what to do in this case if this interested you.
Another option is simply to call it out in the method docs so that people might know this is a "know limitation".

Copy link
Contributor Author

@bratter bratter Nov 3, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Makes sense, adding to the docs would work.. I like not throwing. I'll do in the function signature for here now, and then can do in the website docs later.

Noticed another issue with one of the comparisons that I didn't catch before so need to make a further update, so will amend this PR with both fixes shortly.

Copy link

@HarelM HarelM Nov 3, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good!

`Undefined arc segment, line segment endpoints [[${posA}], [${posB}]] are antipodes`
);
}
}

// Calculate coefficients.
const [D, E, F] = cross(A, B);
const a = E * Cz - F * Cy;
const b = F * Cx - D * Cz;
const c = D * Cy - E * Cx;
// The axis of the great circle passing through the segment's axis and the
// target point
const targetAxis = cross(segmentAxis, C);

const f = c * E - b * F;
const g = a * F - c * D;
const h = b * D - a * E;
// This cross product also has a degenerate case where the segment axis is
// coincident with or antipodal to the target point. In this case the point
// is equidistant to the entire segment. For consistency, we early return the
// endpoint as the matching point.
if (targetAxis[0] === 0 && targetAxis[1] === 0 && targetAxis[2] === 0) {
return [[...posB], true];
}

const t = 1 / Math.sqrt(Math.pow(f, 2) + Math.pow(g, 2) + Math.pow(h, 2));
// The line of intersection between the two great circle planes
const intersectionAxis = cross(targetAxis, segmentAxis);

// Vectors to the two points these great circles intersect.
const I1: Vector = [f * t, g * t, h * t];
const I2: Vector = [-1 * f * t, -1 * g * t, -1 * h * t];
// Vectors to the two points these great circles intersect are the normalized
// intersection and its antipodes
const I1 = normalize(intersectionAxis);
const I2: Vector = [-I1[0], -I1[1], -I1[2]];

// Figure out which is the closest intersection to this segment of the great
// circle.
const angleAB = angle(A, B);
const angleAI1 = angle(A, I1);
const angleBI1 = angle(B, I1);
const angleAI2 = angle(A, I2);
const angleBI2 = angle(B, I2);

let I: Vector;

if (
(angleAI1 < angleAI2 && angleAI1 < angleBI2) ||
(angleBI1 < angleAI2 && angleBI1 < angleBI2)
) {
I = I1;
} else {
I = I2;
}
// Figure out which is the closest intersection to this segment of the great circle
// Note that for points on a unit sphere, the dot product represents the
// cosine of the angle between the two vectors which monotonically increases
// the closer the two points are together and therefore determines proximity
const I = dot(C, I1) > dot(C, I2) ? I1 : I2;

// I is the closest intersection to the segment, though might not actually be
// ON the segment.

// If angle AI or BI is greater than angleAB, I lies on the circle *beyond* A
// and B so use the closest of A or B as the intersection
const angleAB = angle(A, B);
const lngLatI = vectorToLngLat(I);

if (angle(A, I) > angleAB || angle(B, I) > angleAB) {
if (
distance(vectorToLngLat(I), vectorToLngLat(A)) <=
distance(vectorToLngLat(I), vectorToLngLat(B))
) {
return [vectorToLngLat(A), true, false];
// Similar to the usage above, we use the larger dot product to determine
// which endpoint is closer to the test coordinates
// Note that the > means we defer to the endpoint when equidistant,
// following the segment tracking logic in the caller
if (dot(A, C) > dot(B, C)) {
// Clone position when returning as it is reasonable to not expect structural
// sharing on the returned Position in all return cases
return [[...posA], false];
} else {
return [vectorToLngLat(B), false, true];
return [[...posB], true];
}
}

// As angleAI nor angleBI don't exceed angleAB, I is on the segment
return [vectorToLngLat(I), false, false];
return [lngLatI, false];
}

export { nearestPointOnLine };
Expand Down
66 changes: 66 additions & 0 deletions packages/turf-nearest-point-on-line/test.ts
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ import {
point,
featureCollection,
round,
degreesToRadians,
} from "@turf/helpers";
import { nearestPointOnLine } from "./index.js";

Expand Down Expand Up @@ -517,3 +518,68 @@ test("turf-nearest-point-on-line -- issue 2808 redundant point support", (t) =>

t.end();
});

test("turf-nearest-point-on-line -- issue 2934 correct endpoint chosen when in opposite hemisphere", (t) => {
// create a long line where the southern end point should be chosen, but the
// northern endpoint is closer to the projected great circles
const line = lineString([
[25, 88],
[25, -70],
]);
const pt = point([-45, -88]);
const nearest = nearestPointOnLine(line, pt);

t.deepEqual(
truncate(nearest, { precision: 8 }).geometry.coordinates,
[25, -70],
"nearest point is in the Southern Hemisphere"
);

t.end();
});

test("turf-nearest-point-on-line -- issue 2939 nearestPointOnSegment handles tiny segments", (t) => {
// create a test case where the line segment passed is small enough such that
// a cross product used to generate its normal in nearestPointOnSegment ends
// up as [0, 0, 0] but large enough so that the two points are not ===
//
// In v7.2.0 this test case was failing in browser (Chrome 141.0) but
// succeeding on Node (v22.20.0) due to a bit difference in the cosine
// implementation. To reproduce the same effect in Node, we monkey patch
// Math.cos to temporarily return the browser value as it was too difficult to
// find a repro case in node. This monkey patch version fails on v7.2.0 but
// passes with degenerate point fixes.

const lngA = 35.000102519989014;
const lngB = 35.00010251998902;
const line = lineString([
[lngA, 32.00010141921075],
[lngB, 32.00010141921075],
]);

// WARN: Override cos function to give specific results for this test
const originalCos = Math.cos;
Math.cos = (x) => {
// The cosine of lngA gives a different result in the tested browsers than
// it does in node, therefore capturing here for this test run only.
if (x === degreesToRadians(lngA)) {
return originalCos(degreesToRadians(lngB));
} else {
return originalCos(x);
}
};

const pt = point([35.0005, 32.0005]);
const nearest = nearestPointOnLine(line, pt);

// WARN: Reset cos function
Math.cos = originalCos;

t.deepEqual(
nearest.geometry.coordinates,
[35.00010251998902, 32.00010141921075],
"nearest point should be the end point of the line string"
);

t.end();
});