Fixed multiple modules that were performing scalar instead of spherical calculations

packages/turf-line-slice/test.ts

@@ -5,7 +5,7 @@ import { fileURLToPath } from "url";
 import { loadJsonFileSync } from "load-json-file";
 import { writeJsonFileSync } from "write-json-file";
 import { truncate } from "@turf/truncate";
-import { featureCollection } from "@turf/helpers";
+import { featureCollection, point, lineString } from "@turf/helpers";
 import { lineSlice } from "./index.js";
 
 const __dirname = path.dirname(fileURLToPath(import.meta.url));
@@ -38,3 +38,24 @@ test("turf-line-slice", (t) => {
   }
   t.end();
 });
+
+test("turf-nearest-point-on-line -- issue 2023", (t) => {
+  const ptStart = point([3.69140625, 51.72702815704774]);
+  const ptEnd = point([0.31936718356317106, 47.93913163509963]);
+  const line = lineString([
+    [3.69140625, 51.72702815704774],
+    [-5.3173828125, 41.60722821271717],
+  ]);
+
+  const slice = lineSlice(ptStart, ptEnd, line);
+
+  t.deepEqual(
+    truncate(slice, { precision: 8 }).geometry.coordinates,
+    [
+      [3.69140625, 51.72702816],
+      [-0.03079923, 48.08596086],
+    ],
+    "slice should be [[3.69140625, 51.72702816], [-0.03079923, 48.08596086]]"
+  );
+  t.end();
+});

packages/turf-line-slice/test/out/line1.geojson

@@ -38,9 +38,9 @@
       "geometry": {
         "type": "LineString",
         "coordinates": [
-          [-97.835729, 22.247393],
+          [-97.835747, 22.247595],
           [-97.820892, 22.17596],
-          [-97.738467, 22.051207]
+          [-97.738477, 22.051413]
         ]
       }
     }

packages/turf-line-slice/test/out/line2.geojson

@@ -37,8 +37,8 @@
       "geometry": {
         "type": "LineString",
         "coordinates": [
-          [0.049987, 0.049987],
-          [0.849858, 0.849858]
+          [0.05, 0.050008],
+          [0.849981, 0.850017]
         ]
       }
     }

packages/turf-line-slice/test/out/route2.geojson

@@ -3796,7 +3796,7 @@
       "geometry": {
         "type": "LineString",
         "coordinates": [
-          [-111.895792, 48.877416],
+          [-111.895843, 48.877468],
           [-111.878176, 48.860393],
           [-111.867242, 48.849753],
           [-111.866486, 48.849013],

packages/turf-nearest-point-on-line/bench.ts

@@ -2,7 +2,7 @@ import fs from "fs";
 import path from "path";
 import { fileURLToPath } from "url";
 import { loadJsonFileSync } from "load-json-file";
-import Benchmark from "benchmark";
+import Benchmark, { Event } from "benchmark";
 import { nearestPointOnLine } from "./index.js";
 
 const __dirname = path.dirname(fileURLToPath(import.meta.url));
@@ -29,10 +29,19 @@ const fixtures = fs.readdirSync(directory).map((filename) => {
  * route1 x 195 ops/sec ±2.23% (80 runs sampled)
  * route2 x 218 ops/sec ±2.42% (78 runs sampled)
  */
+let totalTime = 0.0;
 const suite = new Benchmark.Suite("turf-nearest-point-on-line");
 for (const { name, geojson } of fixtures) {
   const [line, point] = geojson.features;
-  suite.add(name, () => nearestPointOnLine(line, point));
+  suite.add(name, () => nearestPointOnLine(line, point), {
+    onComplete: (e: Event) =>
+      (totalTime = totalTime += e.target.times?.elapsed),
+  });
 }
 
-suite.on("cycle", (e) => console.log(String(e.target))).run();
+suite
+  .on("cycle", (e: Event) => console.log(String(e.target)))
+  .on("complete", () =>
+    console.log(`completed in ${totalTime.toFixed(2)} seconds`)
+  )
+  .run();

packages/turf-nearest-point-on-line/index.ts

@@ -1,11 +1,8 @@
-import { Feature, Point, LineString, MultiLineString } from "geojson";
-import { bearing } from "@turf/bearing";
+import { Feature, Point, Position, LineString, MultiLineString } from "geojson";
 import { distance } from "@turf/distance";
-import { destination } from "@turf/destination";
-import { lineIntersect as lineIntersects } from "@turf/line-intersect";
 import { flattenEach } from "@turf/meta";
-import { point, lineString, Coord, Units } from "@turf/helpers";
-import { getCoords } from "@turf/invariant";
+import { point, Coord, Units } from "@turf/helpers";
+import { getCoord, getCoords } from "@turf/invariant";
 
 /**
  * Takes a {@link Point} and a {@link LineString} and calculates the closest Point on the (Multi)LineString.
@@ -51,6 +48,8 @@ function nearestPointOnLine<G extends LineString | MultiLineString>(
     throw new Error("lines and pt are required arguments");
   }
 
+  const ptPos = getCoord(pt);
+
   let closestPt: Feature<
     Point,
     { dist: number; index: number; multiFeatureIndex: number; location: number }
@@ -68,80 +67,48 @@ function nearestPointOnLine<G extends LineString | MultiLineString>(
       const coords: any = getCoords(line);
 
       for (let i = 0; i < coords.length - 1; i++) {
-        //start
+        //start - start of current line section
         const start: Feature<Point, { dist: number }> = point(coords[i]);
         start.properties.dist = distance(pt, start, options);
-        //stop
+        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
         const sectionLength = distance(start, stop, options);
-        //perpendicular
-        const heightDistance = Math.max(
-          start.properties.dist,
-          stop.properties.dist
-        );
-        const direction = bearing(start, stop);
-        const perpendicularPt1 = destination(
-          pt,
-          heightDistance,
-          direction + 90,
-          options
-        );
-        const perpendicularPt2 = destination(
-          pt,
-          heightDistance,
-          direction - 90,
-          options
-        );
-        const intersect = lineIntersects(
-          lineString([
-            perpendicularPt1.geometry.coordinates,
-            perpendicularPt2.geometry.coordinates,
-          ]),
-          lineString([start.geometry.coordinates, stop.geometry.coordinates])
-        );
+        let intersectPos: Position;
+        let wasEnd: boolean;
+
+        // Short circuit if snap point is start or end position of the line
+        // segment.
+        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 {
+          // Otherwise, find the nearest point the hard way.
+          [intersectPos, , wasEnd] = nearestPointOnSegment(
+            start.geometry.coordinates,
+            stop.geometry.coordinates,
+            getCoord(pt)
+          );
+        }
         let intersectPt:
           | Feature<
               Point,
               { dist: number; multiFeatureIndex: number; location: number }
             >
           | undefined;
 
-        if (intersect.features.length > 0 && intersect.features[0]) {
-          intersectPt = {
-            ...intersect.features[0],
-            properties: {
-              dist: distance(pt, intersect.features[0], options),
-              multiFeatureIndex: multiFeatureIndex,
-              location:
-                length + distance(start, intersect.features[0], options),
-            },
-          };
-        }
-
-        if (start.properties.dist < closestPt.properties.dist) {
-          closestPt = {
-            ...start,
-            properties: {
-              ...start.properties,
-              index: i,
-              multiFeatureIndex: multiFeatureIndex,
-              location: length,
-            },
-          };
-        }
-
-        if (stop.properties.dist < closestPt.properties.dist) {
-          closestPt = {
-            ...stop,
-            properties: {
-              ...stop.properties,
-              index: i + 1,
-              multiFeatureIndex: multiFeatureIndex,
-              location: length + sectionLength,
-            },
-          };
+        if (intersectPos) {
+          intersectPt = point(intersectPos, {
+            dist: distance(pt, intersectPos, options),
+            multiFeatureIndex: multiFeatureIndex,
+            location: length + distance(start, intersectPos, options),
+          });
         }
 
         if (
@@ -150,9 +117,15 @@ function nearestPointOnLine<G extends LineString | MultiLineString>(
         ) {
           closestPt = {
             ...intersectPt,
-            properties: { ...intersectPt.properties, index: i },
+            properties: {
+              ...intersectPt.properties,
+              // Legacy behaviour where index progresses to next segment # if we
+              // went with the end point this iteration.
+              index: wasEnd ? i + 1 : i,
+            },
           };
         }
+
         // update length
         length += sectionLength;
       }
@@ -162,5 +135,126 @@ function nearestPointOnLine<G extends LineString | MultiLineString>(
   return closestPt;
 }
 
+type Vector = [number, number, number];
+
+function dot(v1: Vector, v2: Vector): number {
+  const [v1x, v1y, v1z] = v1;
+  const [v2x, v2y, v2z] = v2;
+  return v1x * v2x + v1y * v2y + v1z * v2z;
+}
+
+// https://en.wikipedia.org/wiki/Cross_product
+function cross(v1: Vector, v2: Vector): Vector {
+  const [v1x, v1y, v1z] = v1;
+  const [v2x, v2y, v2z] = v2;
+  return [v1y * v2z - v1z * v2y, v1z * v2x - v1x * v2z, v1x * v2y - v1y * v2x];
+}
+
+function magnitude(v: Vector) {
+  return Math.sqrt(Math.pow(v[0], 2) + Math.pow(v[1], 2) + Math.pow(v[2], 2));
+}
+
+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));
+}
+
+function toRadians(degrees: number) {
+  return degrees * (Math.PI / 180);
+}
+
+function toDegrees(radians: number) {
+  return radians * (180 / Math.PI);
+}
+
+function lngLatToVector(a: Position): Vector {
+  const lat = toRadians(a[1]);
+  const lng = toRadians(a[0]);
+  return [
+    Math.cos(lat) * Math.cos(lng),
+    Math.cos(lat) * Math.sin(lng),
+    Math.sin(lat),
+  ];
+}
+
+function vectorToLngLat(v: Vector): Position {
+  const [x, y, z] = v;
+  const lat = toDegrees(Math.asin(z));
+  const lng = toDegrees(Math.atan2(y, x));
+
+  return [lng, lat];
+}
+
+function nearestPointOnSegment(
+  posA: Position,
+  posB: Position,
+  posC: Position
+): [Position, boolean, 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
+
+  // Convert lng/lat to spherical coords
+  const A = lngLatToVector(posA); // the vector from 0,0,0 to posA
+  const B = lngLatToVector(posB);
+  const C = lngLatToVector(posC);
+
+  // Components of target point.
+  const [Cx, Cy, Cz] = C;
+
+  // 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;
+
+  const f = c * E - b * F;
+  const g = a * F - c * D;
+  const h = b * D - a * E;
+
+  const t = 1 / Math.sqrt(Math.pow(f, 2) + Math.pow(g, 2) + Math.pow(h, 2));
+
+  // 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];
+
+  // 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;
+  }
+
+  // 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
+  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];
+    } else {
+      return [vectorToLngLat(B), false, true];
+    }
+  }
+
+  // As angleAI nor angleBI don't exceed angleAB, I is on the segment
+  return [vectorToLngLat(I), false, false];
+}
+
 export { nearestPointOnLine };
 export default nearestPointOnLine;

packages/turf-nearest-point-on-line/package.json

@@ -61,12 +61,9 @@
     "write-json-file": "^5.0.0"
   },
   "dependencies": {
-    "@turf/bearing": "workspace:^",
-    "@turf/destination": "workspace:^",
     "@turf/distance": "workspace:^",
     "@turf/helpers": "workspace:^",
     "@turf/invariant": "workspace:^",
-    "@turf/line-intersect": "workspace:^",
     "@turf/meta": "workspace:^",
     "@types/geojson": "^7946.0.10",
     "tslib": "^2.6.2"