Proposal on top of distance-tool

docs/mcp-tools.md

@@ -2234,7 +2234,7 @@ Renvoie la distance (en mètres) entre deux points à partir de leur longitude e
 | --- | --- | --- | --- |
 | `arrival` | object | oui | Le point d'arrivée |
 | `departure` | object | oui | Le point de départ |
-| `profile` | string (enum) | non | Le type de chemin suivi : `direct` distance à vol d'oiseau (Terre ronde, précision à 0.5%), `vincenty` distance à vol d'oiseau (Terre ellipsoïde, plus précise et coûteuse, précision à 1mm). Par défaut : `direct`. Valeurs : direct, vincenty. Valeur par défaut : direct. |
+| `profile` | string (enum) | non | Le type de chemin suivi : `spherical` distance à vol d'oiseau (Terre ronde, précision à 0.5%), `ellipsoidal` distance à vol d'oiseau (Terre ellipsoïde, plus précise, précision à 0.5cm). Valeurs : spherical, ellipsoidal. Valeur par défaut : spherical. |
 
 <details>
 <summary>Schéma d’entrée brut</summary>
@@ -2292,11 +2292,11 @@ Renvoie la distance (en mètres) entre deux points à partir de leur longitude e
     "profile": {
       "type": "string",
       "enum": [
-        "direct",
-        "vincenty"
+        "spherical",
+        "ellipsoidal"
       ],
-      "default": "direct",
-      "description": "Le type de chemin suivi : `direct` distance à vol d'oiseau (Terre ronde, précision à 0.5%), `vincenty` distance à vol d'oiseau (Terre ellipsoïde, plus précise et coûteuse, précision à 1mm). Par défaut : `direct`."
+      "default": "spherical",
+      "description": "Le type de chemin suivi : `spherical` distance à vol d'oiseau (Terre ronde, précision à 0.5%), `ellipsoidal` distance à vol d'oiseau (Terre ellipsoïde, plus précise, précision à 0.5cm)."
     }
   },
   "required": [

package.json

@@ -62,15 +62,14 @@
     "@turf/bbox-polygon": "^7.4.0",
     "@turf/centroid": "^7.4.0",
     "@turf/circle": "^7.4.0",
-    "@turf/distance": "^7.4.0",
     "@turf/helpers": "^7.4.0",
     "@turf/length": "^7.4.0",
-    "earcut": "^3.2.4",
     "@turf/midpoint": "^7.4.0",
+    "earcut": "^3.2.4",
+    "geographiclib-geodesic": "^2.2.0",
     "jsts": "^2.12.1",
     "mcp-framework": "^0.2.22",
     "node-fetch": "^3.3.2",
-    "node-vincenty": "^0.0.6",
     "winston": "^3.19.0",
     "zod": "^3.25.76"
   },

src/helpers/distance.ts

@@ -3,7 +3,7 @@ import DistanceOp from "jsts/org/locationtech/jts/operation/distance/DistanceOp.
 import IndexedFacetDistance from "jsts/org/locationtech/jts/operation/distance/IndexedFacetDistance.js";
 import RelateOp from "jsts/org/locationtech/jts/operation/relate/RelateOp.js";
 import type { Geometry, Position } from "geojson";
-import { distVincenty } from "node-vincenty";
+import geodesic from "geographiclib-geodesic";
 import { GeometryFactory } from "jsts/org/locationtech/jts/geom.js";
 import GeometryLocation from "jsts/org/locationtech/jts/operation/distance/GeometryLocation.js";
 import { midpoint } from "@turf/midpoint";
@@ -43,7 +43,7 @@ export interface DistanceResult {
 }
 
 /** Point-to-point metric: spherical great-circle, or geodesic on WGS84. */
-export type Metric = "haversine" | "vincenty";
+export type Metric = "spherical" | "ellipsoidal";
 
 type JstsCoord = { x: number; y: number };
 
@@ -62,7 +62,7 @@ type JstsGeometry = {
 const EARTH_RADIUS_M = 6_371_000;
 const DEG_TO_RAD = Math.PI / 180;
 
-/** WGS84 ellipsoid parameters — the ellipsoid `distVincenty` solves on. */
+/** WGS84 ellipsoid parameters, the ellipsoid `ellipsoidalDistance` solves on. */
 const WGS84_SEMI_MAJOR_M = 6_378_137;
 const WGS84_FLATTENING = 1 / 298.257223563;
 const WGS84_ECCENTRICITY_SQ = WGS84_FLATTENING * (2 - WGS84_FLATTENING);
@@ -204,7 +204,7 @@ function nearestOnSegment(p: Position, a: Position, b: Position, pointDistance:
  * that sit within a fraction of a percent of each other.
  */
 function curvatureRadii(lat: number, metric: Metric): { meridional: number; primeVertical: number } {
-  if (metric === "haversine") return { meridional: EARTH_RADIUS_M, primeVertical: EARTH_RADIUS_M };
+  if (metric === "spherical") return { meridional: EARTH_RADIUS_M, primeVertical: EARTH_RADIUS_M };
   const sinLat = Math.sin(lat * DEG_TO_RAD);
   const w = 1 - WGS84_ECCENTRICITY_SQ * sinLat * sinLat;
   return {
@@ -230,7 +230,7 @@ function makeAzimuthalEquidistant(center: Position, metric: Metric) {
   /**
    * Euler's radius of curvature in the normal section of azimuth α, given the
    * unit direction (east, north) = (sin α, cos α) in the tangent plane. Both
-   * radii are equal under "haversine", where this collapses to the sphere.
+   * radii are equal under "spherical", where this collapses to the sphere.
    */
   function radiusInDirection(east: number, north: number): number {
     return 1 / (north * north / meridional + east * east / primeVertical);
@@ -338,7 +338,10 @@ function actualClosestOnGeometryLocation(locA: GeometryLocation, locB: GeometryL
   function getSegment(loc: GeometryLocation) {
     const coords: JstsCoord[] = loc.getGeometryComponent().getCoordinates();
     if (coords.length > 1) {
-      const idx: number = loc.getSegmentIndex();
+      // jsts indexes facets in chunks of 6 segments, so a line of 6k+1 vertices
+      // ends with a single-vertex chunk whose location carries the last vertex
+      // index: its segment is then the last one.
+      const idx = Math.min(loc.getSegmentIndex(), coords.length - 2);
       return { start: unproj(coords[idx]), stop: unproj(coords[idx + 1]) };
     }
   }
@@ -376,10 +379,10 @@ function actualClosestOnGeometryLocation(locA: GeometryLocation, locB: GeometryL
  *
  * @param {object} gA GeoJSON Geometry
  * @param {object} gB GeoJSON Geometry
- * @param metric Point-to-point metric: "haversine" (default) or "vincenty".
+ * @param metric Point-to-point metric: "spherical" (default, haversine) or "ellipsoidal" (geographiclib).
  */
-export function distance(gA: Geometry, gB: Geometry, metric: Metric = "haversine"): DistanceResult {
-  const pointDistance = metric === "vincenty" ? distanceVincenty : haversine;
+export function distance(gA: Geometry, gB: Geometry, metric: Metric = "spherical"): DistanceResult {
+  const pointDistance = metric === "ellipsoidal" ? ellipsoidalDistance : haversine;
   if (gA.type == "Point" && gB.type == "Point") // fast-path for the common case
     return { point1: gA.coordinates, point2: gB.coordinates, distance: Math.round(pointDistance(gA.coordinates, gB.coordinates)*100)/100 };
 
@@ -443,14 +446,14 @@ export function distance(gA: Geometry, gB: Geometry, metric: Metric = "haversine
   return { ...best, distance: Math.round(best.distance * 100) / 100 };
 }
 
-/** Computes Vincenty distance between two lat/lon points through node-vincenty. */
-export function distanceVincenty(a: Position, b: Position) {
-  // distVincenty takes lat1, lon1, lat2, lon2
-  const result = distVincenty(a[1], a[0], b[1], b[0]);
-  if (typeof result !== "object" || result === null) {
-    throw new Error("Vincenty formula failed to converge (antipodal or near-antipodal points)");
-  }
-  return result.distance;
+/**
+ * Computes the geodesic distance on WGS84 between two lon/lat points, with
+ * Karney's algorithm (geographiclib). Unlike Vincenty's iteration, it is
+ * defined for every pair, antipodal points included.
+ */
+export function ellipsoidalDistance(a: Position, b: Position) {
+  // Inverse takes lat1, lon1, lat2, lon2; asking for DISTANCE always sets s12.
+  return geodesic.Geodesic.WGS84.Inverse(a[1], a[0], b[1], b[0], geodesic.Geodesic.DISTANCE).s12!;
 }
 
 export default distance;

src/tools/DistanceTool.ts

@@ -1,5 +1,5 @@
 /**
- * MCP tool exposing distance lookup for a single geographic position.
+ * MCP tool exposing the distance between two geographic positions.
  */
 
 import BaseTool from "./BaseTool.js";
@@ -9,7 +9,7 @@ import { READ_ONLY_OPEN_WORLD_TOOL_ANNOTATIONS } from "../helpers/toolAnnotation
 import { lonSchema, latSchema } from "../helpers/schemas.js";
 import { generatePublishedInputSchema } from "../helpers/jsonSchema.js";
 import logger from "../logger.js";
-import { distanceVincenty, haversine } from "../helpers/distance.js";
+import { ellipsoidalDistance, haversine } from "../helpers/distance.js";
 
 // --- Schemas ---
 
@@ -23,12 +23,11 @@ const distanceInputSchema = z.object({
     lat: latSchema.describe("La latitude du point d'arrivée."),
   }).describe("Le point d'arrivée"),
   profile: z
-    .enum(["direct", "vincenty"])
-    .default("direct")
+    .enum(["spherical", "ellipsoidal"])
+    .default("spherical")
     .describe(["Le type de chemin suivi :",
-      " `direct` distance à vol d'oiseau (Terre ronde, précision à 0.5%),",
-      " `vincenty` distance à vol d'oiseau (Terre ellipsoïde, plus précise et coûteuse, précision à 1mm)",
-      ". Par défaut : `direct`."
+      " `spherical` distance à vol d'oiseau (Terre ronde, précision à 0.5%),",
+      " `ellipsoidal` distance à vol d'oiseau (Terre ellipsoïde, plus précise, précision à 0.5cm).",
     ].join("")),
 }).strict();
 
@@ -71,9 +70,9 @@ class DistanceTool extends BaseTool<DistanceInput> {
     });
 
     switch (input.profile) {
-      case "direct":
-      case "vincenty": {
-        const pointDistance = input.profile == "direct" ? haversine : distanceVincenty;
+      case "spherical":
+      case "ellipsoidal": {
+        const pointDistance = input.profile == "spherical" ? haversine : ellipsoidalDistance;
         const raw = pointDistance(
           [input.departure.lon, input.departure.lat],
           [input.arrival.lon, input.arrival.lat]

test/helpers/distance.test.ts

@@ -1,6 +1,6 @@
 import { describe, expect, it } from "vitest";
 
-import distance, { distanceVincenty, haversine } from "../../src/helpers/distance.js";
+import distance, { ellipsoidalDistance, haversine } from "../../src/helpers/distance.js";
 import { besancon, chamonix, marseille, paris, parisMarseille } from "../samples";
 import type {
   Geometry,
@@ -22,7 +22,7 @@ function ensureSymmetricDistance(
   a: Geometry,
   b: Geometry,
   label: string,
-  metric?: "haversine" | "vincenty",
+  metric?: "spherical" | "ellipsoidal",
 ) {
   const ab = distance(a, b, metric).distance;
   const ba = distance(b, a, metric).distance;
@@ -41,18 +41,24 @@ function expectThrowBothWays(a: Geometry, b: Geometry, error: RegExp) {
 
 describe("distance helper", () => {
   describe("baseline real-world checks", () => {
-    it("supports a Vincenty-backed point metric", () => {
+    it("supports an ellipsoidal point metric", () => {
       const defaultDistance = distance(paris, marseille).distance;
-      const vincentyDistance = ensureSymmetricDistance(paris, marseille, "Paris-Marseille Vincenty", "vincenty");
-      const result = distance(paris, marseille, "vincenty");
-      const expected = Math.round(distanceVincenty(paris.coordinates, marseille.coordinates) * 100) / 100;
+      const ellipsoidal = ensureSymmetricDistance(paris, marseille, "Paris-Marseille ellipsoidal", "ellipsoidal");
+      const result = distance(paris, marseille, "ellipsoidal");
+      const expected = Math.round(ellipsoidalDistance(paris.coordinates, marseille.coordinates) * 100) / 100;
       expect(result.distance).toBe(expected);
-      expect(result.distance).toBe(vincentyDistance);
+      expect(result.distance).toBe(ellipsoidal);
       expect(result.distance).not.toBe(defaultDistance);
       expect(result.point1).toEqual(paris.coordinates);
       expect(result.point2).toEqual(marseille.coordinates);
     });
 
+    it("computes the ellipsoidal distance between antipodal points", () => {
+      // Vincenty's iteration has no answer here; Karney's algorithm goes over the pole.
+      const antipode = [paris.coordinates[0] - 180, -paris.coordinates[1]];
+      expect(ellipsoidalDistance(paris.coordinates, antipode)).toBeCloseTo(20003931.46, 2);
+    });
+
     it("computes Paris->Marseille", () => {
       const result = distance(paris, marseille);
       expectCloseRatio(result.distance, 662_488.38, 0.0005, "Paris-Marseille distance");
@@ -521,6 +527,19 @@ describe("distance helper", () => {
   });
 
   describe("line/segment regression scenarios", () => {
+    it("handles a 7-vertex line whose nearest point is its last vertex", () => {
+      // jsts indexes a line in chunks of 6 segments: 7 vertices leave a
+      // single-vertex chunk at the end, located with the last vertex index.
+      const line: LineString = {
+        type: "LineString",
+        coordinates: Array.from({ length: 7 }, (_, i) => [2 + 0.01 * i, 48.85]),
+      };
+      const end = line.coordinates[6];
+      const point: Point = { type: "Point", coordinates: [end[0] + 0.005, end[1]] };
+      const expected = Math.round(haversine(point.coordinates, end) * 100) / 100;
+      expect(ensureSymmetricDistance(point, line, "point beyond the end of a 7-vertex line")).toBe(expected);
+    });
+
     it("returns zero when line crosses polygon edge", () => {
       const polygon: Polygon = {
         type: "Polygon",
@@ -754,7 +773,7 @@ describe("distance helper", () => {
   // Candidate edges have to be ranked in the same metric the answer is
   // reported in. The planar searches inside the helper run in a projection, and
   // if that projection carries a sphere's local scale while the caller asked
-  // for Vincenty, an edge that is genuinely nearer on WGS84 loses to one that
+  // for the ellipsoid, an edge that is genuinely nearer on WGS84 loses to one that
   // is only nearer on a sphere. The pairs below sit inside the narrow window
   // where the two metrics disagree, so each one fails if the projection and the
   // metric ever drift apart again.
@@ -771,7 +790,7 @@ describe("distance helper", () => {
       { lat: 60, dLon: 1.998318 }, // window (1.996636, 2.000000)
     ];
 
-    it.each(cases)("prefers the parallel edge under Vincenty at lat $lat", ({ lat, dLon }) => {
+    it.each(cases)("prefers the parallel edge on the ellipsoid at lat $lat", ({ lat, dLon }) => {
       const p: Geometry = { type: "Point", coordinates: [0, lat] };
       const edges: Geometry = {
         type: "MultiLineString",
@@ -780,15 +799,15 @@ describe("distance helper", () => {
           [[-dLon, lat + 1], [dLon, lat + 1]],  // parallel edge: nearer on WGS84
         ],
       };
-      const parallelEdge = distanceVincenty([0, lat], [0, lat + 1]);
-      const meridianEdge = distanceVincenty([0, lat], [dLon, lat]);
+      const parallelEdge = ellipsoidalDistance([0, lat], [0, lat + 1]);
+      const meridianEdge = ellipsoidalDistance([0, lat], [dLon, lat]);
       expect(parallelEdge, `lat ${lat}: the case only bites if WGS84 prefers the parallel edge`).toBeLessThan(meridianEdge);
 
-      const d = ensureSymmetricDistance(p, edges, `Vincenty edge ranking at lat ${lat}`, "vincenty");
+      const d = ensureSymmetricDistance(p, edges, `ellipsoidal edge ranking at lat ${lat}`, "ellipsoidal");
       expect(d, `lat ${lat}: expected the parallel edge at ~${parallelEdge.toFixed(2)}, got ${d}`).toBeCloseTo(parallelEdge, 1);
     });
 
-    it.each(cases)("prefers the meridian edge under haversine at lat $lat", ({ lat, dLon }) => {
+    it.each(cases)("prefers the meridian edge on the sphere at lat $lat", ({ lat, dLon }) => {
       const p: Geometry = { type: "Point", coordinates: [0, lat] };
       const edges: Geometry = {
         type: "MultiLineString",
@@ -804,7 +823,7 @@ describe("distance helper", () => {
       // The mirror of the test above: the fix has to follow the requested
       // metric, not hardcode the ellipsoid. A geodesic ranking here would
       // return the parallel edge instead.
-      const d = ensureSymmetricDistance(p, edges, `haversine edge ranking at lat ${lat}`, "haversine");
+      const d = ensureSymmetricDistance(p, edges, `spherical edge ranking at lat ${lat}`, "spherical");
       expect(d, `lat ${lat}: expected the meridian edge at ~${meridianEdge.toFixed(2)}, got ${d}`).toBeLessThan(parallelEdge);
       expectCloseRatio(d, meridianEdge, 0.001, `haversine edge ranking at lat ${lat}`);
     });

test/tools/distance.test.ts

@@ -13,13 +13,13 @@ describe("Test DistanceTool", () => {
         expect(tool.toolDefinition.title).toEqual("Distance entre deux points");
         expect(tool.toolDefinition.inputSchema.required).not.toContain("profile");
         expect(tool.toolDefinition.inputSchema.properties?.profile).toMatchObject({
-            enum: ["direct", "vincenty"],
-            default: "direct",
+            enum: ["spherical", "ellipsoidal"],
+            default: "spherical",
         });
         expect(tool.toolDefinition.outputSchema).toBeDefined();
     });
 
-    it.each([undefined, "direct", "vincenty"])("should return a structured distance for profile %s", async (profile) => {
+    it.each([undefined, "spherical", "ellipsoidal"])("should return a structured distance for profile %s", async (profile) => {
         const tool = new DistanceTool();
         const response = await tool.toolCall({
             params: {

tsconfig.test.json

@@ -12,7 +12,6 @@
   },
   "include": [
     "./test/**/*",
-    "./src/types/**/*.d.ts",
     "./vitest*.mts"
   ],
   "exclude": [