diff --git a/docs/mcp-tools.md b/docs/mcp-tools.md index 16369915..d1c2045c 100644 --- a/docs/mcp-tools.md +++ b/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. |
Schéma d’entrée brut @@ -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": [ diff --git a/package-lock.json b/package-lock.json index 64eca969..37da64d0 100644 --- a/package-lock.json +++ b/package-lock.json @@ -16,15 +16,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", "@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" }, @@ -3414,6 +3413,12 @@ "url": "https://github.com/sponsors/ljharb" } }, + "node_modules/geographiclib-geodesic": { + "version": "2.2.0", + "resolved": "https://registry.npmjs.org/geographiclib-geodesic/-/geographiclib-geodesic-2.2.0.tgz", + "integrity": "sha512-cIedo9VTYb0DFufodgibDmVfsWe9EASqb/kUByl09xc6PZYvLvlc89BHCThtGTPf2OII/zWJGxsR3Uz6O7QOVw==", + "license": "MIT" + }, "node_modules/get-east-asian-width": { "version": "1.7.0", "resolved": "https://registry.npmjs.org/get-east-asian-width/-/get-east-asian-width-1.7.0.tgz", @@ -4717,12 +4722,6 @@ "url": "https://opencollective.com/node-fetch" } }, - "node_modules/node-vincenty": { - "version": "0.0.6", - "resolved": "https://registry.npmjs.org/node-vincenty/-/node-vincenty-0.0.6.tgz", - "integrity": "sha512-oxiqnpfc9LHxm5SqH69WM+rkaIzifZGGvoer3AFyWEipNLJ4LhurwaDBhvgXglootehVs8iKWDMfLapPki+esg==", - "license": "BSD" - }, "node_modules/npm-run-path": { "version": "6.0.0", "resolved": "https://registry.npmjs.org/npm-run-path/-/npm-run-path-6.0.0.tgz", diff --git a/package.json b/package.json index 7d8d2a35..59d39955 100644 --- a/package.json +++ b/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" }, diff --git a/src/helpers/distance.ts b/src/helpers/distance.ts index c46f7447..5031108b 100644 --- a/src/helpers/distance.ts +++ b/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; diff --git a/src/tools/DistanceTool.ts b/src/tools/DistanceTool.ts index f8e9e209..b4a660a5 100644 --- a/src/tools/DistanceTool.ts +++ b/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 { }); 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] diff --git a/src/types/node-vincenty.d.ts b/src/types/node-vincenty.d.ts deleted file mode 100644 index 3884c090..00000000 --- a/src/types/node-vincenty.d.ts +++ /dev/null @@ -1,15 +0,0 @@ -declare module "node-vincenty" { - export interface VincentyDistanceResult { - distance: number; - initialBearing: number; - finalBearing: number; - } - - export function distVincenty( - lat1: number, - lon1: number, - lat2: number, - lon2: number, - callback?: (distance: number, initialBearing?: number, finalBearing?: number) => void, - ): VincentyDistanceResult | number | null; -} diff --git a/test/helpers/distance.test.ts b/test/helpers/distance.test.ts index 5e3e5e75..6227045d 100644 --- a/test/helpers/distance.test.ts +++ b/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}`); }); diff --git a/test/tools/distance.test.ts b/test/tools/distance.test.ts index fa749949..1f175d18 100644 --- a/test/tools/distance.test.ts +++ b/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: { diff --git a/tsconfig.test.json b/tsconfig.test.json index 65098b24..8118f6bb 100644 --- a/tsconfig.test.json +++ b/tsconfig.test.json @@ -12,7 +12,6 @@ }, "include": [ "./test/**/*", - "./src/types/**/*.d.ts", "./vitest*.mts" ], "exclude": [