From bcca3cba0486249296d9e2c7f742999cbb601995 Mon Sep 17 00:00:00 2001 From: Lionel Zoubritzky Date: Thu, 18 Jun 2026 09:23:41 +0200 Subject: [PATCH 01/11] fix: Improve distance computation --- package-lock.json | 50 ++- package.json | 3 + src/gpf/parcellaire-express.ts | 2 +- src/gpf/urbanisme.ts | 4 +- src/helpers/distance.ts | 395 +++++++++++++++++-- src/types/node-vincenty.d.ts | 15 + src/wfs/spatialExtras.ts | 2 +- test/helpers/distance.test.ts | 677 +++++++++++++++++++++++++++++++-- test/wfs/response.test.ts | 4 +- tsconfig.test.json | 1 + 10 files changed, 1073 insertions(+), 80 deletions(-) create mode 100644 src/types/node-vincenty.d.ts diff --git a/package-lock.json b/package-lock.json index 3b03e050..eb88b2ed 100644 --- a/package-lock.json +++ b/package-lock.json @@ -23,6 +23,8 @@ "jsts": "^2.12.1", "mcp-framework": "^0.2.22", "node-fetch": "^3.3.2", + "node-vincenty": "^0.0.6", + "rbush": "^4.0.1", "winston": "^3.19.0", "zod": "^3.25.76" }, @@ -41,6 +43,7 @@ "@turf/intersect": "^7.4.0", "@types/geojson": "^7946.0.16", "@types/node": "^26.6.3", + "@types/rbush": "^4.0.0", "@types/supertest": "^7.2.1", "@vitest/coverage-v8": "^5.0.3", "ajv": "^8.20.0", @@ -992,22 +995,25 @@ } }, "node_modules/@napi-rs/wasm-runtime": { - "version": "1.1.6", - "resolved": "https://registry.npmjs.org/@napi-rs/wasm-runtime/-/wasm-runtime-1.1.6.tgz", - "integrity": "sha512-ZLv/JdUfkvOy9eCnnBaGfiO+XimbjebAeO+MRQqD/B+FR1tnRN0tpKSJHRbE8sFfS6aqsXZ67TQjfwfsxULVbg==", + "version": "1.2.4", + "resolved": "https://registry.npmjs.org/@napi-rs/wasm-runtime/-/wasm-runtime-1.2.4.tgz", + "integrity": "sha512-AJxoUD2/15ESHbvpcyjU274nsAPLuOtPHCk0vKJM5pj//Fg/B1FXNWjPnXTT9PymCYYiHo4zPj0ZomXBKhoy7g==", "dev": true, "license": "MIT", "optional": true, "dependencies": { "@tybys/wasm-util": "^0.10.3" }, + "engines": { + "node": "^20.19.0 || ^22.13.0 || >=23.5.0" + }, "funding": { "type": "github", "url": "https://github.com/sponsors/Brooooooklyn" }, "peerDependencies": { - "@emnapi/core": "^1.7.1", - "@emnapi/runtime": "^1.7.1" + "@emnapi/core": "^1.7.1 || ^2.0.0-alpha.4", + "@emnapi/runtime": "^1.7.1 || ^2.0.0-alpha.4" } }, "node_modules/@noble/hashes": { @@ -1547,9 +1553,9 @@ } }, "node_modules/@tybys/wasm-util": { - "version": "0.10.3", - "resolved": "https://registry.npmjs.org/@tybys/wasm-util/-/wasm-util-0.10.3.tgz", - "integrity": "sha512-F3fo1MYrRJYL3zER0OUOmkutjr1Vp23m7OsSgp7nq4SP6OqX6C/56XFIPAl5bt3zaBRjmW7SGz3u/6LwFpYcOg==", + "version": "0.10.4", + "resolved": "https://registry.npmjs.org/@tybys/wasm-util/-/wasm-util-0.10.4.tgz", + "integrity": "sha512-W3c4gRigFS0T/Ma4qIYF3GDAc5AQdHb1yL5znJT1Zv1YaD9Kitx656wBjvr19qbiosmZT8lWDM5BEMynUqX65A==", "dev": true, "license": "MIT", "optional": true, @@ -1644,6 +1650,13 @@ "kleur": "^3.0.3" } }, + "node_modules/@types/rbush": { + "version": "4.0.0", + "resolved": "https://registry.npmjs.org/@types/rbush/-/rbush-4.0.0.tgz", + "integrity": "sha512-+N+2H39P8X+Hy1I5mC6awlTX54k3FhiUmvt7HWzGJZvF+syUAAxP/stwppS8JE84YHqFgRMv6fCy31202CMFxQ==", + "dev": true, + "license": "MIT" + }, "node_modules/@types/superagent": { "version": "8.1.11", "resolved": "https://registry.npmjs.org/@types/superagent/-/superagent-8.1.11.tgz", @@ -4672,6 +4685,12 @@ "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", @@ -5186,6 +5205,12 @@ "dev": true, "license": "MIT" }, + "node_modules/quickselect": { + "version": "3.0.0", + "resolved": "https://registry.npmjs.org/quickselect/-/quickselect-3.0.0.tgz", + "integrity": "sha512-XdjUArbK4Bm5fLLvlm5KpTFOiOThgfWWI4axAZDWg4E/0mKdZyI9tNEfds27qCi1ze/vwTR16kvmmGhRra3c2g==", + "license": "ISC" + }, "node_modules/range-parser": { "version": "1.2.1", "resolved": "https://registry.npmjs.org/range-parser/-/range-parser-1.2.1.tgz", @@ -5210,6 +5235,15 @@ "node": ">= 0.10" } }, + "node_modules/rbush": { + "version": "4.0.1", + "resolved": "https://registry.npmjs.org/rbush/-/rbush-4.0.1.tgz", + "integrity": "sha512-IP0UpfeWQujYC8Jg162rMNc01Rf0gWMMAb2Uxus/Q0qOFw4lCcq6ZnQEZwUoJqWyUGJ9th7JjwI4yIWo+uvoAQ==", + "license": "MIT", + "dependencies": { + "quickselect": "^3.0.0" + } + }, "node_modules/react": { "version": "19.3.0", "resolved": "https://registry.npmjs.org/react/-/react-19.3.0.tgz", diff --git a/package.json b/package.json index df08b4d9..a9dd89a4 100644 --- a/package.json +++ b/package.json @@ -69,6 +69,8 @@ "jsts": "^2.12.1", "mcp-framework": "^0.2.22", "node-fetch": "^3.3.2", + "node-vincenty": "^0.0.6", + "rbush": "^4.0.1", "winston": "^3.19.0", "zod": "^3.25.76" }, @@ -83,6 +85,7 @@ "@turf/intersect": "^7.4.0", "@types/geojson": "^7946.0.16", "@types/node": "^26.6.3", + "@types/rbush": "^4.0.0", "@types/supertest": "^7.2.1", "@vitest/coverage-v8": "^5.0.3", "ajv": "^8.20.0", diff --git a/src/gpf/parcellaire-express.ts b/src/gpf/parcellaire-express.ts index 90a11411..d063e5e9 100644 --- a/src/gpf/parcellaire-express.ts +++ b/src/gpf/parcellaire-express.ts @@ -104,7 +104,7 @@ export async function getParcellaireExpress(lon: number, lat: number): Promise

180°) are rejected: per RFC 7946 they SHOULD be split before being + * passed here (Polygon → MultiPolygon, LineString → MultiLineString). + * Antimeridian-adjacent geometry *pairs* (individually valid, but facing each + * other across ±180°) are rejected too: the planar nearest-point search below + * cannot resolve them correctly. + * + * The implementation follows three steps: + * 1. ask JTS for an exact planar nearest-point seed and an early zero-distance check, + * 2. decompose both geometries into vertices and segments with an R-tree index, + * 3. refine the best geodesic answer by projecting vertices onto candidate segments, + * pruning most edges through cheap bounds and the R-tree index. + */ + +export interface DistanceResult { + distance: number; + point1: Position; + point2: Position; +} + +const EARTH_RADIUS_M = 6_371_000; +const DEG_TO_RAD = Math.PI / 180; +// Smallest WGS84 meridional meters per latitude degree (near the equator). +// We use this conservative floor where bounds must never overestimate distance: +// pruning lower-bounds and latitude expansion in search boxes. +const MIN_LAT_METERS_PER_DEGREE = 110_574; +// Spherical meters per degree from EARTH_RADIUS_M. +// Used for spherical approximations (haversine and lon/lat scaling heuristics). +const SPHERE_METERS_PER_DEGREE = EARTH_RADIUS_M * DEG_TO_RAD; + +interface Edge { + a: Position; + b: Position; + minX: number; + maxX: number; + minY: number; + maxY: number; + lonMetersPerDegree: number; +} + +interface Shape { + vertices: Position[]; + edges: Edge[]; + edgeIndex: RBush | null; +} + +interface NearestPoint { + distance: number; + point: Position; +} + +/** Returns the great-circle distance between two lon/lat positions in meters. */ +export function haversine(a: Position, b: Position): number { + const phi1 = a[1]*DEG_TO_RAD; + const phi2 = b[1]*DEG_TO_RAD; + const dPhi = phi2 - phi1; + const dLambda = (b[0] - a[0])*DEG_TO_RAD; + const h = Math.sin(dPhi / 2) ** 2 + Math.cos(phi1) * Math.cos(phi2) * Math.sin(dLambda / 2) ** 2; + return EARTH_RADIUS_M * 2 * Math.atan2(Math.sqrt(h), Math.sqrt(1 - h)); +} + +/** + * Throws if any ring (Polygon/MultiPolygon) or segment (LineString/MultiLineString) + * has a longitude jump > 180°, indicating an antimeridian-crossing geometry that + * SHOULD be split per RFC 7946 before being used here. + */ +function assertNoAntimeridian(geom: Geometry): void { + function checkCoords(coords: Position[]): void { + for (let i = 0; i < coords.length - 1; i++) { + if (Math.abs(coords[i + 1][0] - coords[i][0]) > 180) + throw new Error( + "Antimeridian-crossing geometries SHOULD be split per RFC 7946 " + + "(Polygon → MultiPolygon, LineString → MultiLineString)", + ); + } + } + switch (geom.type) { + case "LineString": checkCoords(geom.coordinates); break; + case "MultiLineString": geom.coordinates.forEach(checkCoords); break; + case "Polygon": geom.coordinates.forEach(checkCoords); break; + case "MultiPolygon": geom.coordinates.forEach(poly => poly.forEach(checkCoords)); break; + case "GeometryCollection": geom.geometries.forEach(assertNoAntimeridian); break; + } +} + +/** Returns the [min, max] longitude spanned by a set of vertices. */ +function lonExtent(vertices: Position[]): [number, number] { + let min = Infinity, max = -Infinity; + for (const [lon] of vertices) { + if (lon < min) min = lon; + if (lon > max) max = lon; + } + return [min, max]; +} + +/** + * Throws if gA and gB's vertex extents are closer going the other way around + * ±180° than directly: JTS's planar nearest-point search and the R-tree pruning + * below assume a flat lon/lat plane and would then silently pick the wrong, + * far-side pairing (see the module docstring). + */ +function assertNoAntimeridianPair(verticesA: Position[], verticesB: Position[]): void { + const [aMin, aMax] = lonExtent(verticesA); + const [bMin, bMax] = lonExtent(verticesB); + const rawGap = Math.max(0, bMin - aMax, aMin - bMax); + const wrappedGap = 360 - Math.max(aMax, bMax) + Math.min(aMin, bMin); + if (wrappedGap < rawGap) + throw new Error( + "Antimeridian-adjacent geometry pair: the shortest path appears to cross ±180° " + + "longitude, which this planar-nearest-point implementation cannot resolve correctly", + ); +} + +/** Precomputes one segment with its bbox and longitude scaling metadata. */ +function makeEdge(a: Position, b: Position): Edge { + const minX = Math.min(a[0], b[0]); + const maxX = Math.max(a[0], b[0]); + const minY = Math.min(a[1], b[1]); + const maxY = Math.max(a[1], b[1]); + const maxAbsLat = Math.max(Math.abs(minY), Math.abs(maxY)); + return { a, b, minX, maxX, minY, maxY, lonMetersPerDegree: SPHERE_METERS_PER_DEGREE * Math.cos(maxAbsLat*DEG_TO_RAD) }; +} /** - * Compute approximative distance in meters between gA and gB. - * - * TODO: replace the lon/lat planar nearest-point step with a geodesic geometry distance. + * Decomposes a geometry into its vertices and edges, then builds a search Shape. + * An R-tree index is created for shapes with more than 64 edges. + */ +function buildShape(geom: Geometry): Shape { + const vertices: Position[] = []; + const edges: Edge[] = []; + + function process(g: Geometry): void { + switch (g.type) { + case "Point": vertices.push(g.coordinates); break; + case "MultiPoint": vertices.push(...g.coordinates); break; + case "LineString": + vertices.push(...g.coordinates); + for (let i = 0; i < g.coordinates.length - 1; i++) edges.push(makeEdge(g.coordinates[i], g.coordinates[i + 1])); + break; + case "MultiLineString": + for (const line of g.coordinates) { + vertices.push(...line); + for (let i = 0; i < line.length - 1; i++) edges.push(makeEdge(line[i], line[i + 1])); + } + break; + case "Polygon": + for (const ring of g.coordinates) { + vertices.push(...ring); + for (let i = 0; i < ring.length - 1; i++) edges.push(makeEdge(ring[i], ring[i + 1])); + } + break; + case "MultiPolygon": + for (const poly of g.coordinates) + for (const ring of poly) { + vertices.push(...ring); + for (let i = 0; i < ring.length - 1; i++) edges.push(makeEdge(ring[i], ring[i + 1])); + } + break; + case "GeometryCollection": + for (const child of g.geometries) process(child); + break; + } + } + + process(geom); + const edgeIndex = edges.length > 64 // R-tree build cost is ~O(n log n); only worthwhile beyond this edge count + ? new RBush().load(edges) : null; + return { vertices, edges, edgeIndex }; +} + +/** Interpolates a point at parameter t along a lon/lat segment. */ +function pointOnSegment(a: Position, b: Position, t: number): Position { + return [a[0] + (b[0] - a[0]) * t, a[1] + (b[1] - a[1]) * t]; +} + +/** Maps one lon/lat position to its unit vector on the Earth sphere. */ +function toUnitVector([lon, lat]: Position): [number, number, number] { + const phi = lat*DEG_TO_RAD; + const lambda = lon*DEG_TO_RAD; + return [Math.cos(phi) * Math.cos(lambda), Math.cos(phi) * Math.sin(lambda), Math.sin(phi)]; +} + +/** Projects a point onto a segment in 3D chord space and clamps the result to [0,1]. */ +function chordProjection(p: Position, a: Position, b: Position): number { + const [px, py, pz] = toUnitVector(p); + const [ax, ay, az] = toUnitVector(a); + const [bx, by, bz] = toUnitVector(b); + const dx = bx - ax, dy = by - ay, dz = bz - az; + const length2 = dx * dx + dy * dy + dz * dz; + if (length2 < 1e-18) return 0; // squared chord length ≈ 0 → degenerate segment, return t=0 + const t = ((px - ax) * dx + (py - ay) * dy + (pz - az) * dz) / length2; + return Math.min(1, Math.max(0, t)); +} + +/** Refines the nearest point from p onto one segment until cm / 0.5% precision is reached. */ +function nearestOnSegment(p: Position, a: Position, b: Position, pointDistance: (a: Position, b: Position) => number): NearestPoint { + const distanceAt = (t: number) => pointDistance(p, pointOnSegment(a, b, t)); + + // Start from the orthogonal projection in 3D chord space, then refine locally + // on the actual lon/lat segment with progressively smaller steps. + let t = chordProjection(p, a, b); + let best = distanceAt(t); + const segmentLengthMeters = pointDistance(a, b); + const arcAngle = segmentLengthMeters / EARTH_RADIUS_M; + const requiredPrecisionMeters = Math.max(0.01, 0.005 * best); // stop when step < 0.5% of best distance, floored at 1 cm + + for ( + let step = Math.min(0.5, 0.2 * arcAngle); // initial step: 0.2 of the arc angle (radians), capped at 0.5 t-units (half segment) + step * segmentLengthMeters > requiredPrecisionMeters; + step /= 2 + ) { + let improved = true; + while (improved) { + improved = false; + for (const dt of [step, -step]) { + const candidate = Math.min(1, Math.max(0, t + dt)); + const d = distanceAt(candidate); + if (d < best) { + best = d; + t = candidate; + improved = true; + } + } + } + } + + return { distance: best, point: pointOnSegment(a, b, t) }; +} + +/** Computes a cheap metric lower bound from a point to an edge bbox. */ +function lowerBoundToEdge(p: Position, edge: Edge, pLonMetersPerDegree: number = SPHERE_METERS_PER_DEGREE*Math.cos(p[1]*DEG_TO_RAD)): number { + const latGap = Math.max(0, edge.minY - p[1], p[1] - edge.maxY); + const lonGap = Math.max(0, edge.minX - p[0], p[0] - edge.maxX); + const lonMetersPerDegree = Math.min(edge.lonMetersPerDegree, pLonMetersPerDegree); + return Math.max(latGap * MIN_LAT_METERS_PER_DEGREE, lonGap * lonMetersPerDegree); +} + +/** + * Expands a metric radius around a point into a conservative lon/lat search box. + * + * latDegrees uses the same MIN_LAT_METERS_PER_DEGREE floor as lowerBoundToEdge, so + * it never underestimates the latitude span needed to reach `meters` away. + * + * lonDegrees solves the exact haversine equation for the worst-case latitude in + * [p.lat - latDegrees, p.lat + latDegrees] (the one closest to a pole, where a + * degree of longitude is shortest). The approximation can under-estimate + * the needed longitude span as soon as the box reaches into higher latitudes + * than p itself, letting the true nearest edge fall outside the search box. + */ +function safeSearchBox(p: Position, meters: number): Pick { + const latDegrees = meters / MIN_LAT_METERS_PER_DEGREE; + const maxAbsLatDeg = Math.min(89.999, Math.max(Math.abs(p[1] - latDegrees), Math.abs(p[1] + latDegrees))); + const cosProduct = Math.cos(p[1]*DEG_TO_RAD) * Math.cos(maxAbsLatDeg*DEG_TO_RAD); + + let lonDegrees = 180; + if (cosProduct > 1e-9) { + const halfAngle = meters / (2 * EARTH_RADIUS_M); + const sinSqHalfDeltaLambda = Math.sin(halfAngle) ** 2 / cosProduct; + if (sinSqHalfDeltaLambda >= 1) + lonDegrees = Math.min(180, (2 * Math.asin(Math.sqrt(sinSqHalfDeltaLambda))) / DEG_TO_RAD); + } + + return { + minX: p[0] - lonDegrees, + minY: p[1] - latDegrees, + maxX: p[0] + lonDegrees, + maxY: p[1] + latDegrees, + }; +} + + +/** Finds the nearest point from p to a shape using vertex fallback and edge pruning. Returns null if nothing beats the upper bound. */ +function nearestOnShape(p: Position, shape: Shape, pointDistance: (a: Position, b: Position) => number, upperBound = Infinity): NearestPoint | null { + let best: NearestPoint | null = null; + + if (shape.edges.length === 0) { + for (const v of shape.vertices) { + const d = pointDistance(p, v); + if (d < upperBound && (best === null || d < best.distance)) { + best = { distance: d, point: v }; + } + } + return best; + } + + const pLonMetersPerDegree = SPHERE_METERS_PER_DEGREE * Math.cos(p[1]*DEG_TO_RAD); + const candidates = shape.edgeIndex + ? shape.edgeIndex.search(safeSearchBox(p, Math.min(upperBound, 20_000_000))) // 20 Mm ≈ half-Earth circumference: caps the bbox so safeSearchBox never overflows ±180°/±90° + : shape.edges; + + // Each edge first gets a cheap lower-bound test from its bbox in meters. + // Only the survivors pay for the more expensive point-to-segment refinement. + for (const edge of candidates) { + const bound = best === null ? upperBound : best.distance; + if (lowerBoundToEdge(p, edge, pLonMetersPerDegree) >= bound) continue; + const nearest = nearestOnSegment(p, edge.a, edge.b, pointDistance); + if (nearest.distance < (best === null ? upperBound : best.distance)) best = nearest; + } + + return best; +} + +/** + * Compute geodesic distance in meters between gA and gB, with closest points. + * + * JTS gives us a robust planar nearest-point seed and detects all exact touches, + * overlaps and containments. If the planar answer is not already zero, the final + * geodesic result is refined by scanning both directions: vertices of A against + * shape B, then vertices of B against shape A. * * @param {object} gA GeoJSON Geometry * @param {object} gB GeoJSON Geometry + * @param metric Point-to-point metric: "haversine" (default) or "vincenty". */ -export default function distance(gA: Geometry, gB: Geometry): number { - const geojsonReader = new GeoJSONReader(new GeometryFactory()) - const a = geojsonReader.read(gA); - const b = geojsonReader.read(gB); - - /* - * Get the 2 nearest points between a and b - * - * Note that it will project according to longitude and latitude axis, - * so it is not really accurate, but it is a good approximation - */ - const nearestPoints = DistanceOp.nearestPoints(a, b); - if ( nearestPoints.length !== 2 ) { - throw new Error('DistanceOp.nearestPoints should return 2 points'); +export function distance(gA: Geometry, gB: Geometry, metric: "haversine" | "vincenty" = "haversine"): DistanceResult { + const pointDistance = metric === "vincenty" ? distanceVincenty : haversine; + + assertNoAntimeridian(gA); + assertNoAntimeridian(gB); + + const shapeA = buildShape(gA); + const shapeB = buildShape(gB); + assertNoAntimeridianPair(shapeA.vertices, shapeB.vertices); + + const reader = new (GeoJSONReader as unknown as new () => { read(g: Geometry): unknown })(); + type JstsCoord = { x: number; y: number }; + const [n1, n2] = DistanceOp.nearestPoints(reader.read(gA), reader.read(gB)) as [JstsCoord, JstsCoord]; + const p1: Position = [n1.x, n1.y]; + const p2: Position = [n2.x, n2.y]; + + if (Math.hypot(n1.x - n2.x, n1.y - n2.y) < 1e-9) { // ~0.1 mm in degrees: below JTS floating-point noise → exact touch + return { distance: 0, point1: p1, point2: p1 }; + } + + let best: DistanceResult = { distance: pointDistance(p1, p2), point1: p1, point2: p2 }; + + for (const v of shapeA.vertices) { + const nearest = nearestOnShape(v, shapeB, pointDistance, best.distance); + if (nearest && nearest.distance < best.distance) { + best = { distance: nearest.distance, point1: v, point2: nearest.point }; } + } - /* - * haversine distance between the 2 nearest points (see https://turfjs.org/docs/api/distance) - */ - return turfDistance( - turfPoint([nearestPoints[0].x, nearestPoints[0].y]), - turfPoint([nearestPoints[1].x, nearestPoints[1].y]), - { units: 'meters' } - ); + for (const v of shapeB.vertices) { + const nearest = nearestOnShape(v, shapeA, pointDistance, best.distance); + if (nearest && nearest.distance < best.distance) { + best = { distance: nearest.distance, point1: nearest.point, point2: v }; + } + } + + 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; } + +export default distance; diff --git a/src/types/node-vincenty.d.ts b/src/types/node-vincenty.d.ts new file mode 100644 index 00000000..3884c090 --- /dev/null +++ b/src/types/node-vincenty.d.ts @@ -0,0 +1,15 @@ +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/src/wfs/spatialExtras.ts b/src/wfs/spatialExtras.ts index dd29ff55..2f226161 100644 --- a/src/wfs/spatialExtras.ts +++ b/src/wfs/spatialExtras.ts @@ -284,7 +284,7 @@ export function deriveFromGeometry(geometry: unknown, input: FeatureCollectionPo if (requires_distance_to_filter_center) { try { const filterCentroid = context.filterCentroid!; - ret.distance_to_filter_center = distance(geo, filterCentroid); + ret.distance_to_filter_center = distance(geo, filterCentroid).distance; } catch { ret.distance_to_filter_center = null; } diff --git a/test/helpers/distance.test.ts b/test/helpers/distance.test.ts index 09f9586f..124f9853 100644 --- a/test/helpers/distance.test.ts +++ b/test/helpers/distance.test.ts @@ -1,39 +1,646 @@ import { describe, expect, it } from "vitest"; -import distance from "../../src/helpers/distance.js"; -import type { Polygon } from "geojson"; -import {paris, marseille, besancon, parisMarseille} from '../samples'; - -describe("Test distance",() => { - describe("Test distance(Point,Point)", () => { - it("should return 662489.3m from Paris to Marseille",() => { - const result = distance(paris,marseille); - expect(result).toBeCloseTo(662489.3,1); - }); - }); - describe("Test distance(Point,LineString)", () => { - it("should return 209731.2m from Besançon to [Paris,Marseille]",() => { - const result = distance(besancon,parisMarseille); - expect(result).toBeCloseTo(209731.2,1); - }); - }); - - describe("Test distance(Point,Polygon)", () => { - it("should return 0m from Paris point to a polygon containing Paris",() => { - const polygonContainingParis: Polygon = { - "type": "Polygon", - "coordinates": [ - [ - [2.0, 48.0], - [3.0, 48.0], - [3.0, 49.0], - [2.0, 49.0], - [2.0, 48.0] - ] - ] - }; - const result = distance(paris, polygonContainingParis); - expect(result).toEqual(0); - }); +import distance, { distanceVincenty, haversine } from "../../src/helpers/distance.js"; +import { besancon, chamonix, marseille, paris, parisMarseille } from "../samples"; +import type { + Geometry, + GeometryCollection, + LineString, + MultiLineString, + MultiPoint, + MultiPolygon, + Point, + Polygon, +} from "geojson"; + +function expectCloseRatio(actual: number, expected: number, ratio: number, label: string) { + const tolerance = Math.abs(expected) * ratio; + expect(Math.abs(actual - expected), `${label}: expected ~${expected}, got ${actual}`).toBeLessThanOrEqual(tolerance); +} + +function ensureSymmetricDistance( + a: Geometry, + b: Geometry, + label: string, + metric?: "haversine" | "vincenty", +) { + const ab = distance(a, b, metric).distance; + const ba = distance(b, a, metric).distance; + expect(ab, `${label}: A->B/B->A mismatch`).toBeCloseTo(ba, 6); + return ab; +} + +function expectZeroBothWays(a: Geometry, b: Geometry, label: string) { + expect(ensureSymmetricDistance(a, b, label)).toBe(0); +} + +function expectThrowBothWays(a: Geometry, b: Geometry, error: RegExp) { + expect(() => distance(a, b)).toThrow(error); + expect(() => distance(b, a)).toThrow(error); +} + +describe("distance helper", () => { + describe("baseline real-world checks", () => { + it("supports a Vincenty-backed 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; + expect(result.distance).toBe(expected); + expect(result.distance).toBe(vincentyDistance); + expect(result.distance).not.toBe(defaultDistance); + expect(result.point1).toEqual(paris.coordinates); + expect(result.point2).toEqual(marseille.coordinates); + }); + + it("computes Paris->Marseille", () => { + const result = distance(paris, marseille); + expectCloseRatio(result.distance, 662_488.38, 0.0005, "Paris-Marseille distance"); + expect(result.point1).toEqual(paris.coordinates); + expect(result.point2).toEqual(marseille.coordinates); + }); + + it("computes Besancon->(Paris,Marseille) with point fixed on Besancon", () => { + const result = distance(besancon, parisMarseille); + expectCloseRatio(result.distance, 198_521.68, 0.002, "Besancon-line distance"); + expect(result.point1).toEqual(besancon.coordinates); + }); + + it("computes Paris->Lyon", () => { + const lyon: Geometry = { type: "Point", coordinates: [4.835, 45.764] }; + expectCloseRatio(distance(paris, lyon).distance, 392_000, 0.02, "Paris-Lyon distance"); + }); + + it("computes Eiffel Tower->Notre-Dame", () => { + const eiffelTower: Geometry = { type: "Point", coordinates: [2.2945, 48.8584] }; + const notreDame: Geometry = { type: "Point", coordinates: [2.3499, 48.853] }; + expectCloseRatio(distance(eiffelTower, notreDame).distance, 4_100, 0.05, "Eiffel-Notre-Dame distance"); + }); + + it("computes point near Loire to simplified Loire path", () => { + const loire: Geometry = { + type: "LineString", + coordinates: [ + [1.9, 47.9], + [1.33, 47.59], + [0.68, 47.39], + ], + }; + const nearbyPoint: Geometry = { type: "Point", coordinates: [1.9, 48.0] }; + const result = distance(nearbyPoint, loire); + expect(result.distance).toBeGreaterThan(0); + expectCloseRatio(result.distance, 11_100, 0.15, "point-to-Loire distance"); + }); + + it("keeps Corsica-mainland Var distance in plausible range", () => { + const corsica: Geometry = { + type: "Polygon", + coordinates: [[[8.5, 41.3], [9.6, 41.3], [9.6, 43.0], [8.5, 43.0], [8.5, 41.3]]], + }; + const mainlandVar: Geometry = { + type: "Polygon", + coordinates: [[[6.0, 43.0], [6.9, 43.0], [6.9, 43.3], [6.0, 43.3], [6.0, 43.0]]], + }; + const result = distance(corsica, mainlandVar); + expect(result.distance).toBeGreaterThan(100_000); + expect(result.distance).toBeLessThan(220_000); + }); + }); + + describe("containment and intersections", () => { + it("returns zero for point inside polygon", () => { + const point: Geometry = { type: "Point", coordinates: [5, 5] }; + const polygon: Geometry = { + type: "Polygon", + coordinates: [[[0, 0], [10, 0], [10, 10], [0, 10], [0, 0]]], + }; + expect(distance(point, polygon).distance).toBe(0); + }); + + it("returns positive for point in polygon hole", () => { + const point: Geometry = { type: "Point", coordinates: [5, 5] }; + const polygonWithHole: Geometry = { + type: "Polygon", + coordinates: [ + [[0, 0], [10, 0], [10, 10], [0, 10], [0, 0]], + [[4, 4], [4, 6], [6, 6], [6, 4], [4, 4]], + ], + }; + expect(distance(point, polygonWithHole).distance).toBeGreaterThan(0); + }); + + it("returns zero for crossing lines with shared closest point", () => { + const l1: Geometry = { type: "LineString", coordinates: [[0, 0], [10, 10]] }; + const l2: Geometry = { type: "LineString", coordinates: [[0, 10], [10, 0]] }; + const result = distance(l1, l2); + expectZeroBothWays(l1, l2, "crossing lines"); + expect(result.point1).toEqual([5, 5]); + expect(result.point2).toEqual([5, 5]); + }); + + it("returns zero for GeometryCollections that intersect", () => { + const containingPolygon: Polygon = { + type: "Polygon", + coordinates: [[[2.0, 48.0], [3.0, 48.0], [3.0, 49.0], [2.0, 49.0], [2.0, 48.0]]], + }; + + const gcA: GeometryCollection = { + type: "GeometryCollection", + geometries: [paris, { type: "LineString", coordinates: [[0, 0], [0.5, 0.5]] }], + }; + const gcB: GeometryCollection = { type: "GeometryCollection", geometries: [containingPolygon] }; + expectZeroBothWays(gcA, gcB, "intersecting geometry collections"); + }); + + it("does not misclassify outside point for polygon with duplicate vertex", () => { + const polygonWithDuplicateVertex: Polygon = { + type: "Polygon", + coordinates: [[[0, 0], [2, 0], [2, 0], [2, 2], [0, 2], [0, 0]]], + }; + const pointOutside: Point = { type: "Point", coordinates: [10, 10] }; + const result = distance(pointOutside, polygonWithDuplicateVertex); + expect(ensureSymmetricDistance(pointOutside, polygonWithDuplicateVertex, "duplicate vertex polygon outside point")).toBeGreaterThan(0); + expect(result.point1).toEqual([10, 10]); + }); + }); + + describe("antimeridian and poles", () => { + // Antimeridian-crossing geometries are rejected by design (RFC 7946 SHOULD split them). + const capVertexCount = 36; + const ring = Array.from({ length: capVertexCount }, (_, i) => + [-180 + (i * 360) / capVertexCount, 80] as [number, number], + ); + ring.push(ring[0]); + + it.each([ + { + name: "throws for a polygon that crosses the antimeridian", + g1: { type: "Point", coordinates: [172, 5] } as Geometry, + g2: { + type: "Polygon", + coordinates: [ + [[170, 0], [170, 10], [-170, 10], [-170, 0], [170, 0]], + [[175, 3], [175, 7], [178, 7], [178, 3], [175, 3]], + ], + } as Geometry, + error: /RFC 7946/, + }, + { + name: "throws for a polar-cap polygon whose closing edge spans > 180° of longitude", + g1: { type: "Point", coordinates: [0, 85] } as Geometry, + g2: { type: "Polygon", coordinates: [ring] } as Geometry, + error: /RFC 7946/, + }, + { + name: "throws for a pair of individually valid geometries facing each other across the antimeridian", + g1: { + type: "Polygon", + coordinates: [[[179.8, -0.1], [179.9, -0.1], [179.9, 0.1], [179.8, 0.1], [179.8, -0.1]]], + } as Geometry, + g2: { + type: "Polygon", + coordinates: [[[-179.9, -0.1], [-179.8, -0.1], [-179.8, 0.1], [-179.9, 0.1], [-179.9, -0.1]]], + } as Geometry, + error: /Antimeridian-adjacent/, + }, + ])("$name", ({ g1, g2, error }) => expectThrowBothWays(g1, g2, error)); + }); + + describe("antimeridian-safe Multi* geometries (RFC 7946 conforming)", () => { + it.each([ + { + name: "accepts a MultiPolygon with one sub-polygon on each side of the antimeridian", + g1: { type: "Point", coordinates: [179, 2.5] } as Geometry, + g2: { + type: "MultiPolygon", + coordinates: [ + [[[178, 0], [180, 0], [180, 5], [178, 5], [178, 0]]], + [[[-180, 0], [-178, 0], [-178, 5], [-180, 5], [-180, 0]]], + ], + } as Geometry, + }, + { + name: "accepts a MultiLineString representing a line cut at the antimeridian", + g1: { type: "Point", coordinates: [179, 2] } as Geometry, + g2: { + type: "MultiLineString", + coordinates: [ + [[178, 2], [180, 2]], + [[-180, 2], [-175, 2]], + ], + } as Geometry, + }, + ])("$name", ({ g1, g2 }) => { + expectZeroBothWays(g1, g2, "antimeridian-safe multi geometry"); + }); + }); + + describe("composite geometries", () => { + it("is symmetric for MultiLineString vs MultiPolygon", () => { + const multiLine: MultiLineString = { + type: "MultiLineString", + coordinates: [parisMarseille.coordinates, [[7.5, 46.0], [7.8, 46.3]]], + }; + const multiPolygon: MultiPolygon = { + type: "MultiPolygon", + coordinates: [ + [[[1.8, 48.2], [2.8, 48.2], [2.8, 49.1], [1.8, 49.1], [1.8, 48.2]]], + [[[6.6, 45.6], [7.1, 45.6], [7.1, 46.0], [6.6, 46.0], [6.6, 45.6]]], + ], + }; + + ensureSymmetricDistance(multiLine, multiPolygon, "MultiLineString/MultiPolygon"); + }); + + it("handles MultiPoint fallback", () => { + const a: MultiPoint = { + type: "MultiPoint", + coordinates: [paris.coordinates, marseille.coordinates], + }; + const b: MultiPoint = { + type: "MultiPoint", + coordinates: [chamonix.coordinates, paris.coordinates], + }; + expectZeroBothWays(a, b, "MultiPoint fallback"); + }); + + it("uses nearest piece for MultiPolygon", () => { + const multiPolygon: Geometry = { + type: "MultiPolygon", + coordinates: [ + [[[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]]], + [[[5, 5], [6, 5], [6, 6], [5, 6], [5, 5]]], + ], + }; + const p: Geometry = { type: "Point", coordinates: [5.5, 5.5] }; + expectZeroBothWays(p, multiPolygon, "nearest piece for MultiPolygon"); + }); + + it("uses nearest member for GeometryCollection", () => { + const collection: Geometry = { + type: "GeometryCollection", + geometries: [ + { type: "Point", coordinates: [0, 0] }, + { type: "LineString", coordinates: [[20, 20], [21, 21]] }, + ], + }; + const p: Geometry = { type: "Point", coordinates: [0.01, 0] }; + expectCloseRatio(ensureSymmetricDistance(p, collection, "nearest collection member"), 1_112, 0.02, "nearest collection member"); + }); + + it("visits isolated Point in GeometryCollection when A has no edges", () => { + const a: Geometry = { type: "Point", coordinates: [0, 80] }; + const b: Geometry = { + type: "GeometryCollection", + geometries: [ + { type: "Point", coordinates: [3, 80] }, // ~57.9 km away + { type: "LineString", coordinates: [[-5, 82], [5, 82]] }, // ~222 km away + ], + }; + const result = distance(a, b); + expectCloseRatio(ensureSymmetricDistance(a, b, "isolated Point in GC"), 57_919.97, 0.01, "isolated Point in GC"); + expectCloseRatio(result.distance, 57_919.97, 0.01, "isolated Point in GC"); + expect(result.point2).toEqual([3, 80]); + }); + }); + + describe("zero-distance edge cases", () => { + const cases: Array<{ name: string; g1: Geometry; g2: Geometry }> = [ + { + name: "point on polygon vertex", + g1: { type: "Point", coordinates: [0, 0] }, + g2: { type: "Polygon", coordinates: [[[0, 0], [10, 0], [10, 10], [0, 10], [0, 0]]] }, + }, + { + name: "point on polygon edge", + g1: { type: "Point", coordinates: [0.5, 0] }, + g2: { type: "Polygon", coordinates: [[[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]]] }, + }, + { + name: "two lines sharing an endpoint", + g1: { type: "LineString", coordinates: [[0, 0], [1, 1]] }, + g2: { type: "LineString", coordinates: [[1, 1], [2, 0]] }, + }, + { + name: "two polygons sharing an edge", + g1: { type: "Polygon", coordinates: [[[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]]] }, + g2: { type: "Polygon", coordinates: [[[1, 0], [2, 0], [2, 1], [1, 1], [1, 0]]] }, + }, + { + name: "two polygons touching at one vertex", + g1: { type: "Polygon", coordinates: [[[0, 0], [1, 0], [0, 1], [0, 0]]] }, + g2: { type: "Polygon", coordinates: [[[1, 0], [2, 0], [1, 1], [1, 0]]] }, + }, + { + name: "multipolygon touching multilinestring", + g1: { + type: "MultiPolygon", + coordinates: [ + [[[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]]], + [[[10, 10], [11, 10], [11, 11], [10, 11], [10, 10]]], + ], + }, + g2: { + type: "MultiLineString", + coordinates: [ + [[20, 20], [21, 21]], + [[10.5, 9], [10.5, 12]], + ], + }, + }, + { + name: "geometry collection member touching external point", + g1: { type: "Point", coordinates: [1, 0] }, + g2: { + type: "GeometryCollection", + geometries: [ + { type: "Point", coordinates: [50, 50] }, + { type: "Polygon", coordinates: [[[-1, -1], [1, -1], [1, 1], [-1, 1], [-1, -1]]] }, + ], + }, + }, + { + name: "point on multipolygon far-piece vertex", + g1: { type: "Point", coordinates: [101, 101] }, + g2: { + type: "MultiPolygon", + coordinates: [ + [[[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]]], + [[[100, 100], [101, 100], [101, 101], [100, 101], [100, 100]]], + ], + }, + }, + ]; + + it.each(cases)("returns zero for $name", ({ g1, g2, name }) => expectZeroBothWays(g1, g2, name)); + }); + + describe("line/segment regression scenarios", () => { + it("returns zero when line crosses polygon edge", () => { + const polygon: Polygon = { + type: "Polygon", + coordinates: [[[1, 1], [5, 1], [5, 5], [1, 5], [1, 1]]], + }; + const crossingLine: LineString = { + type: "LineString", + coordinates: [[0, 3], [6, 3]], + }; + const result = distance(crossingLine, polygon); + expectZeroBothWays(crossingLine, polygon, "line crosses polygon edge"); + expect(result.point1).toEqual(result.point2); + expect(result.point2[1]).toBeCloseTo(3, 2); + }); + + it("handles parallel non-intersecting lines", () => { + const line1: LineString = { type: "LineString", coordinates: [[1, 1], [11, 1]] }; + const line2: LineString = { type: "LineString", coordinates: [[1, 6], [11, 6]] }; + const result = distance(line1, line2); + expect(ensureSymmetricDistance(line1, line2, "parallel non-intersecting lines")).toBeGreaterThan(0); + expect(result.point1[0]).toBeGreaterThanOrEqual(1); + expect(result.point1[0]).toBeLessThanOrEqual(11); + expect(result.point2[0]).toBeGreaterThanOrEqual(1); + expect(result.point2[0]).toBeLessThanOrEqual(11); + }); + + it("handles multiple line segments with intersection in middle", () => { + const line1: LineString = { type: "LineString", coordinates: [[0, 0], [5, 5], [10, 10]] }; + const line2: LineString = { type: "LineString", coordinates: [[0, 10], [5, 5], [10, 0]] }; + const result = distance(line1, line2); + expectZeroBothWays(line1, line2, "multiple line segments with intersection in middle"); + expect(result.point1).toEqual(result.point2); + }); + + it("returns concrete non-dummy intersection coordinates", () => { + const geom1: LineString = { type: "LineString", coordinates: [[2.0, 48.0], [3.0, 49.0]] }; + const geom2: Polygon = { + type: "Polygon", + coordinates: [[[2.5, 48.0], [3.5, 48.0], [3.5, 49.0], [2.5, 49.0], [2.5, 48.0]]], + }; + const result = distance(geom1, geom2); + expectZeroBothWays(geom1, geom2, "concrete non-dummy intersection coordinates"); + expect(result.point1).toEqual(result.point2); + expect(result.point1).not.toEqual([0, 0]); + }); + + it("handles colinear overlapping lines", () => { + const line1: LineString = { type: "LineString", coordinates: [[0, 0], [10, 0]] }; + const line2: LineString = { type: "LineString", coordinates: [[5, 0], [15, 0]] }; + const result = distance(line1, line2); + expectZeroBothWays(line1, line2, "colinear overlapping lines"); + expect(result.point1[0]).toBeCloseTo(result.point2[0], 10); + expect(result.point1[1]).toBeCloseTo(result.point2[1], 10); + expect(result.point1[1]).toBeCloseTo(0, 10); + expect(result.point1[0]).toBeGreaterThanOrEqual(5); + expect(result.point1[0]).toBeLessThanOrEqual(10); + }); + + it("handles colinear non-overlapping lines", () => { + const line1: LineString = { type: "LineString", coordinates: [[0, 0], [5, 0]] }; + const line2: LineString = { type: "LineString", coordinates: [[10, 0], [15, 0]] }; + const result = distance(line1, line2); + expect(ensureSymmetricDistance(line1, line2, "colinear non-overlapping lines")).toBeGreaterThan(0); + expect(result.point1[1]).toBe(0); + expect(result.point2[1]).toBe(0); + const point1XInRange = (result.point1[0] >= 0 && result.point1[0] <= 5) || (result.point1[0] >= 10 && result.point1[0] <= 15); + const point2XInRange = (result.point2[0] >= 0 && result.point2[0] <= 5) || (result.point2[0] >= 10 && result.point2[0] <= 15); + expect(point1XInRange && point2XInRange).toBe(true); + }); + + it("handles line identical to polygon edge", () => { + const polygon: Polygon = { + type: "Polygon", + coordinates: [[[0, 0], [10, 0], [10, 10], [0, 10], [0, 0]]], + }; + const line: LineString = { type: "LineString", coordinates: [[0, 0], [10, 0]] }; + const result = distance(line, polygon); + expectZeroBothWays(line, polygon, "line identical to polygon edge"); + expect(result.point1).toEqual(result.point2); + }); + + it.each([ + { + name: "throws for a LineString that crosses the antimeridian", + g1: { + type: "Polygon", + coordinates: [[ + [23.135776547816135, -5.4896438887560635], + [27.977340911114993, -5.4896438887560635], + [27.977340911114993, 31.956263865426763], + [23.135776547816135, 31.956263865426763], + [23.135776547816135, -5.4896438887560635], + ]], + } as Geometry, + g2: { + type: "LineString", + coordinates: [ + [165.2812012336696, 6.549152703997294], + [-174.51797643760415, 6.852029146295976], + ], + } as Geometry, + }, + { + name: "throws for a LineString that crosses the antimeridian (former artifact regression)", + g1: { + type: "LineString", + coordinates: [ + [148.36866917715702, 46.84158577670763], + [92.27258167379557, 58.34804763832423], + ], + } as Geometry, + g2: { + type: "LineString", + coordinates: [ + [-157.2713800668627, -4.447939722741566], + [169.80694578648507, -14.439383953401865], + ], + } as Geometry, + }, + ])("$name", ({ g1, g2 }) => expectThrowBothWays(g1, g2, /RFC 7946/)); + + it("does not follow a naive s1/s2 result", () => { + const s1: Geometry = { + type: "LineString", + coordinates: [ + [-136.61127697878203, -67.24854737502018], + [3.6575704307094554, -65.92384272707795], + ], + }; + const s2: Geometry = { + type: "LineString", + coordinates: [ + [-69.05775139070103, -54.8979093430909], + [-86.16447630165361, -57.89895966720661], + ], + }; + + const result = distance(s1, s2); + const d = ensureSymmetricDistance(s1, s2, "s1/s2 distance"); + expectCloseRatio(d, 986_440, 0.005, "s1/s2 distance"); + expect(d).toBeLessThan(2_440_829.44); + expectCloseRatio(result.point1[0], -85.7, 0.02, "s1/s2 point1 lon"); + expectCloseRatio(result.point1[1], -66.77, 0.02, "s1/s2 point1 lat"); + expectCloseRatio(result.point2[0], -86.16447630165361, 0.001, "s1/s2 point2 lon"); + expectCloseRatio(result.point2[1], -57.89895966720661, 0.001, "s1/s2 point2 lat"); + }); + + it("finds nearest edge with Vincenty when lower bounds use haversine METERS_PER_DEGREE", () => { + const p: Geometry = { type: "Point", coordinates: [0, 0] }; + const b: Geometry = { + type: "MultiLineString", + coordinates: [ + [[0.9963, -1], [0.9963, 1]], // planar seed ~110907 m, lon-edge + [[-1, 1.0], [1, 1.0]], // true nearest ~110574 m, lat-edge + ], + }; + expectCloseRatio( + ensureSymmetricDistance(p, b, "Vincenty avoids over-pruning lat bounds", "vincenty"), + 110574.39, + 0.001, + "Vincenty avoids over-pruning lat bounds", + ); + }); + + it("regresses the concrete RBush false-pruning case with a poleward nearest edge", () => { + const point: Geometry = { type: "Point", coordinates: [-10, 72] }; + const ring = Array.from({ length: 66 }, (_, i) => { + const angle = (2 * Math.PI * i) / 65; + return [ + -2.5 + Math.cos(angle), + Math.max(-89, Math.min(89, 72.5 + 4 * Math.sin(angle))), + ] as [number, number]; + }); + const shape: Geometry = { + type: "MultiLineString", + coordinates: Array.from({ length: 65 }, (_, i) => [ring[i], ring[i + 1]]), + }; + + const result = distance(point, shape); + + expect(result.distance).toBe(223_111.88); + expect(result.point1).toEqual(point.coordinates); + expect(result.point2[0]).toBeLessThan(-3.49); + expect(result.point2[0]).toBeGreaterThan(-3.5); + expect(result.point2[1]).toBeGreaterThan(72.12); + expect(result.point2[1]).toBeLessThan(72.13); + }); + + // The previous test pins one witnessed failure. This one is broader: it fuzzes + // many high-latitude >64-edge shapes and checks the general invariant that the + // RBush-pruned result must not exceed the brute-force vertex minimum. It guards + // against nearby variants of the same under-sized search-box bug. + it("fuzzes the high-latitude RBush invariant for many >64-edge shapes", () => { + let lcg = 424242; + const rnd = () => (lcg = (lcg * 1103515245 + 12345) % 2147483648) / 2147483648; + const ring = (cLon: number, cLat: number, r: number, n: number, seed: number): Polygon => ({ + type: "Polygon", + coordinates: [Array.from({ length: n + 1 }, (_, i) => { + const a = (2 * Math.PI * i) / n; + const w = 1 + 0.3 * Math.sin(a * 5 + seed); + return [cLon + r * w * Math.cos(a), Math.max(-89, Math.min(89, cLat + r * w * Math.sin(a)))]; + })], + }); + + for (let k = 0; k < 40; k++) { + const latA = 40 + rnd() * 48, latB = 40 + rnd() * 48; + const lonA = -170 + rnd() * 40, lonB = 100 + rnd() * 70; + const a = ring(lonA, latA, 1 + rnd() * 8, 80, rnd() * 6); + const b = ring(lonB, latB, 1 + rnd() * 8, 80, rnd() * 6); + + let got: number; + try { got = distance(a, b).distance; } catch { continue; } + + let bruteForce = Infinity; + for (const va of a.coordinates[0]) for (const vb of b.coordinates[0]) bruteForce = Math.min(bruteForce, haversine(va, vb)); + + expect(got, `iteration ${k}: got=${got}, brute-force vertex minimum=${bruteForce}`).toBeLessThanOrEqual(bruteForce * 1.000001); + } + }); + }); + + describe("performance sanity", () => { + it("resolves two disjoint 500-vertex polygons in under a second", () => { + const circle = (centerLon: number, centerLat: number, radiusDeg: number, n: number): Geometry => ({ + type: "Polygon", + coordinates: [ + Array.from({ length: n + 1 }, (_, i) => { + const a = (2 * Math.PI * i) / n; + return [centerLon + radiusDeg * Math.cos(a), centerLat + radiusDeg * Math.sin(a)]; + }), + ], + }); + + const a = circle(2.35, 48.85, 0.3, 500); + const b = circle(2.35, 53.85, 0.3, 500); + + const start = performance.now(); + const result = distance(a, b); + const elapsedMs = performance.now() - start; + + expect(result.distance).toBeGreaterThan(0); + expect(elapsedMs).toBeLessThan(1_000); + }); + + it("resolves two nearby irregular ~5000-vertex polygons in under a second", () => { + const blob = (centerLon: number, centerLat: number, radiusDeg: number, n: number, seed: number): Geometry => ({ + type: "Polygon", + coordinates: [ + Array.from({ length: n + 1 }, (_, i) => { + const a = (2 * Math.PI * i) / n; + const wobble = 1 + 0.15 * Math.sin(a * 7 + seed) + 0.08 * Math.sin(a * 13 + seed * 2); + return [centerLon + radiusDeg * wobble * Math.cos(a), centerLat + radiusDeg * wobble * Math.sin(a)]; + }), + ], + }); + + const a = blob(2.0, 48.0, 0.5, 5000, 1); + const b = blob(4.0, 48.3, 0.5, 5000, 2); + + const start = performance.now(); + const result = distance(a, b); + const elapsedMs = performance.now() - start; + + expect(result.distance).toBeGreaterThan(0); + expect(elapsedMs).toBeLessThan(1_000); }); + }); }); diff --git a/test/wfs/response.test.ts b/test/wfs/response.test.ts index a503b4a1..2aac45fe 100644 --- a/test/wfs/response.test.ts +++ b/test/wfs/response.test.ts @@ -180,7 +180,7 @@ describe("wfs_engine/response", () => { }); const features = getFeatures(result); - expect(features[0].distance_to_filter_center as number).toBeCloseTo(1330.6551992128234, 6); + expect(features[0].distance_to_filter_center as number).toBeCloseTo(1330.65, 3); expect(features[0].intersection_area as number).toBeCloseTo(13010026.445506852, 3); }); @@ -209,7 +209,7 @@ describe("wfs_engine/response", () => { }); const features = getFeatures(result); - expect(features[0].distance_to_filter_center as number).toBeCloseTo(2340.9971606708805, 6); + expect(features[0].distance_to_filter_center as number).toBeCloseTo(2340.99, 3); }); it("should stay fast when computing centroid, area, and intersection_area for a large region against many polygons crossing its boundary", () => { diff --git a/tsconfig.test.json b/tsconfig.test.json index 8118f6bb..65098b24 100644 --- a/tsconfig.test.json +++ b/tsconfig.test.json @@ -12,6 +12,7 @@ }, "include": [ "./test/**/*", + "./src/types/**/*.d.ts", "./vitest*.mts" ], "exclude": [ From 7c15992d0cfab5474434a1584be5021372cd1656 Mon Sep 17 00:00:00 2001 From: Lionel Zoubritzky Date: Tue, 18 Aug 2026 19:35:02 +0200 Subject: [PATCH 02/11] feat: distance use azimuthal projections --- package-lock.json | 49 ++-- package.json | 2 +- src/helpers/distance.ts | 476 ++++++++++++++++++++-------------- test/helpers/distance.test.ts | 289 ++++++++++++++++++--- 4 files changed, 564 insertions(+), 252 deletions(-) diff --git a/package-lock.json b/package-lock.json index eb88b2ed..64eca969 100644 --- a/package-lock.json +++ b/package-lock.json @@ -19,12 +19,12 @@ "@turf/distance": "^7.4.0", "@turf/helpers": "^7.4.0", "@turf/length": "^7.4.0", + "@turf/midpoint": "^7.4.0", "earcut": "^3.2.4", "jsts": "^2.12.1", "mcp-framework": "^0.2.22", "node-fetch": "^3.3.2", "node-vincenty": "^0.0.6", - "rbush": "^4.0.1", "winston": "^3.19.0", "zod": "^3.25.76" }, @@ -1418,6 +1418,21 @@ "url": "https://opencollective.com/turf" } }, + "node_modules/@turf/bearing": { + "version": "7.4.0", + "resolved": "https://registry.npmjs.org/@turf/bearing/-/bearing-7.4.0.tgz", + "integrity": "sha512-utMyjTU5U3QfHWLRIKUAvJerCMmC6NLI7X6cmL6PA5+KYdzul/ymxL6CQrwJugrf9zt4Hnzvipwz0j5H1Yo87w==", + "license": "MIT", + "dependencies": { + "@turf/helpers": "7.4.0", + "@turf/invariant": "7.4.0", + "@types/geojson": "^7946.0.10", + "tslib": "^2.8.1" + }, + "funding": { + "url": "https://opencollective.com/turf" + } + }, "node_modules/@turf/centroid": { "version": "7.4.0", "resolved": "https://registry.npmjs.org/@turf/centroid/-/centroid-7.4.0.tgz", @@ -1552,6 +1567,23 @@ "url": "https://opencollective.com/turf" } }, + "node_modules/@turf/midpoint": { + "version": "7.4.0", + "resolved": "https://registry.npmjs.org/@turf/midpoint/-/midpoint-7.4.0.tgz", + "integrity": "sha512-+poTGThn9F2tveuUDjg1z5TE3QG8zp/IHVqpzkFwFLWX3xpeys0mev9H3LqlbzRU7waawZJ/lsqYlHg3mIWHJw==", + "license": "MIT", + "dependencies": { + "@turf/bearing": "7.4.0", + "@turf/destination": "7.4.0", + "@turf/distance": "7.4.0", + "@turf/helpers": "7.4.0", + "@types/geojson": "^7946.0.10", + "tslib": "^2.8.1" + }, + "funding": { + "url": "https://opencollective.com/turf" + } + }, "node_modules/@tybys/wasm-util": { "version": "0.10.4", "resolved": "https://registry.npmjs.org/@tybys/wasm-util/-/wasm-util-0.10.4.tgz", @@ -5205,12 +5237,6 @@ "dev": true, "license": "MIT" }, - "node_modules/quickselect": { - "version": "3.0.0", - "resolved": "https://registry.npmjs.org/quickselect/-/quickselect-3.0.0.tgz", - "integrity": "sha512-XdjUArbK4Bm5fLLvlm5KpTFOiOThgfWWI4axAZDWg4E/0mKdZyI9tNEfds27qCi1ze/vwTR16kvmmGhRra3c2g==", - "license": "ISC" - }, "node_modules/range-parser": { "version": "1.2.1", "resolved": "https://registry.npmjs.org/range-parser/-/range-parser-1.2.1.tgz", @@ -5235,15 +5261,6 @@ "node": ">= 0.10" } }, - "node_modules/rbush": { - "version": "4.0.1", - "resolved": "https://registry.npmjs.org/rbush/-/rbush-4.0.1.tgz", - "integrity": "sha512-IP0UpfeWQujYC8Jg162rMNc01Rf0gWMMAb2Uxus/Q0qOFw4lCcq6ZnQEZwUoJqWyUGJ9th7JjwI4yIWo+uvoAQ==", - "license": "MIT", - "dependencies": { - "quickselect": "^3.0.0" - } - }, "node_modules/react": { "version": "19.3.0", "resolved": "https://registry.npmjs.org/react/-/react-19.3.0.tgz", diff --git a/package.json b/package.json index a9dd89a4..7d8d2a35 100644 --- a/package.json +++ b/package.json @@ -66,11 +66,11 @@ "@turf/helpers": "^7.4.0", "@turf/length": "^7.4.0", "earcut": "^3.2.4", + "@turf/midpoint": "^7.4.0", "jsts": "^2.12.1", "mcp-framework": "^0.2.22", "node-fetch": "^3.3.2", "node-vincenty": "^0.0.6", - "rbush": "^4.0.1", "winston": "^3.19.0", "zod": "^3.25.76" }, diff --git a/src/helpers/distance.ts b/src/helpers/distance.ts index d0730c8f..c46f7447 100644 --- a/src/helpers/distance.ts +++ b/src/helpers/distance.ts @@ -1,8 +1,12 @@ -import RBush from "rbush"; import GeoJSONReader from "jsts/org/locationtech/jts/io/GeoJSONReader.js"; import DistanceOp from "jsts/org/locationtech/jts/operation/distance/DistanceOp.js"; +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 { GeometryFactory } from "jsts/org/locationtech/jts/geom.js"; +import GeometryLocation from "jsts/org/locationtech/jts/operation/distance/GeometryLocation.js"; +import { midpoint } from "@turf/midpoint"; /** * Minimal geodesic distance between two GeoJSON geometries, reported in meters @@ -11,15 +15,25 @@ import { distVincenty } from "node-vincenty"; * Antimeridian-crossing geometries (any Polygon ring or LineString segment with * |Δlon| > 180°) are rejected: per RFC 7946 they SHOULD be split before being * passed here (Polygon → MultiPolygon, LineString → MultiLineString). - * Antimeridian-adjacent geometry *pairs* (individually valid, but facing each - * other across ±180°) are rejected too: the planar nearest-point search below - * cannot resolve them correctly. + * + * Edges interpolate linearly in lon/lat, so a "point on the geometry" always + * means a point on that linear path. * * The implementation follows three steps: - * 1. ask JTS for an exact planar nearest-point seed and an early zero-distance check, - * 2. decompose both geometries into vertices and segments with an R-tree index, - * 3. refine the best geodesic answer by projecting vertices onto candidate segments, - * pruning most edges through cheap bounds and the R-tree index. + * 1. settle the zero-distance cases (touch, crossing, overlap, containment), + * skipping the test entirely when the two bounding boxes are apart, + * 2. seed the search with an exact planar nearest pair, computed in a plane + * whose two axes carry comparable ground distances, + * 3. refine that answer by re-running the planar search in an azimuthal + * equidistant projection centered on the current best pair. + * + * The search is monotone rather than convergent: step 3 only ever adopts a pair + * that beats the incumbent (with a 5mm tolerance), so the result is the best of + * the candidates seen. Every candidate is a genuine (point on A, point on B) + * pair, hence an upper bound on the true distance. + * + * Every planar search goes through an STR-tree of facets, so the cost stays + * near-linear in the vertex count (see `planarNearestLocations`). */ export interface DistanceResult { @@ -28,31 +42,39 @@ export interface DistanceResult { point2: Position; } +/** Point-to-point metric: spherical great-circle, or geodesic on WGS84. */ +export type Metric = "haversine" | "vincenty"; + +type JstsCoord = { x: number; y: number }; + +/** + * A jsts geometry. jsts ships its geometry classes untyped, so the members this + * file relies on are declared here rather than threaded as `any`. + */ +type JstsGeometry = { + getEnvelopeInternal(): { + distance(other: unknown): number; + getMinY(): number; + getMaxY(): number; + }; +}; + const EARTH_RADIUS_M = 6_371_000; const DEG_TO_RAD = Math.PI / 180; -// Smallest WGS84 meridional meters per latitude degree (near the equator). -// We use this conservative floor where bounds must never overestimate distance: -// pruning lower-bounds and latitude expansion in search boxes. -const MIN_LAT_METERS_PER_DEGREE = 110_574; -// Spherical meters per degree from EARTH_RADIUS_M. -// Used for spherical approximations (haversine and lon/lat scaling heuristics). -const SPHERE_METERS_PER_DEGREE = EARTH_RADIUS_M * DEG_TO_RAD; - -interface Edge { - a: Position; - b: Position; - minX: number; - maxX: number; - minY: number; - maxY: number; - lonMetersPerDegree: number; -} -interface Shape { - vertices: Position[]; - edges: Edge[]; - edgeIndex: RBush | null; -} +/** WGS84 ellipsoid parameters — the ellipsoid `distVincenty` 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); + +/** + * Refinement stops once a pass moves the answer by less than half a centimeter. + * Results are reported to the centimeter, so chasing more than that only buys + * extra planar searches. + */ +const TOLERANCE_M = 5e-3; + +const MAX_REFINEMENT_PASSES = 10; // 10 for security, but 3 is generally enough interface NearestPoint { distance: number; @@ -70,100 +92,39 @@ export function haversine(a: Position, b: Position): number { } /** - * Throws if any ring (Polygon/MultiPolygon) or segment (LineString/MultiLineString) - * has a longitude jump > 180°, indicating an antimeridian-crossing geometry that - * SHOULD be split per RFC 7946 before being used here. + * Throws if the edge has a longitude jump > 180°, indicating an + * antimeridian-crossing geometry that SHOULD be split per RFC 7946. */ -function assertNoAntimeridian(geom: Geometry): void { - function checkCoords(coords: Position[]): void { - for (let i = 0; i < coords.length - 1; i++) { - if (Math.abs(coords[i + 1][0] - coords[i][0]) > 180) - throw new Error( - "Antimeridian-crossing geometries SHOULD be split per RFC 7946 " + - "(Polygon → MultiPolygon, LineString → MultiLineString)", - ); - } - } - switch (geom.type) { - case "LineString": checkCoords(geom.coordinates); break; - case "MultiLineString": geom.coordinates.forEach(checkCoords); break; - case "Polygon": geom.coordinates.forEach(checkCoords); break; - case "MultiPolygon": geom.coordinates.forEach(poly => poly.forEach(checkCoords)); break; - case "GeometryCollection": geom.geometries.forEach(assertNoAntimeridian); break; - } -} - -/** Returns the [min, max] longitude spanned by a set of vertices. */ -function lonExtent(vertices: Position[]): [number, number] { - let min = Infinity, max = -Infinity; - for (const [lon] of vertices) { - if (lon < min) min = lon; - if (lon > max) max = lon; - } - return [min, max]; -} - -/** - * Throws if gA and gB's vertex extents are closer going the other way around - * ±180° than directly: JTS's planar nearest-point search and the R-tree pruning - * below assume a flat lon/lat plane and would then silently pick the wrong, - * far-side pairing (see the module docstring). - */ -function assertNoAntimeridianPair(verticesA: Position[], verticesB: Position[]): void { - const [aMin, aMax] = lonExtent(verticesA); - const [bMin, bMax] = lonExtent(verticesB); - const rawGap = Math.max(0, bMin - aMax, aMin - bMax); - const wrappedGap = 360 - Math.max(aMax, bMax) + Math.min(aMin, bMin); - if (wrappedGap < rawGap) +function checkEdge(c1: Position, c2: Position) { + if (Math.abs(c1[0] - c2[0]) > 180) { throw new Error( - "Antimeridian-adjacent geometry pair: the shortest path appears to cross ±180° " + - "longitude, which this planar-nearest-point implementation cannot resolve correctly", + "Antimeridian-crossing geometries SHOULD be split per RFC 7946 " + + "(Polygon → MultiPolygon, LineString → MultiLineString)", ); + } } -/** Precomputes one segment with its bbox and longitude scaling metadata. */ -function makeEdge(a: Position, b: Position): Edge { - const minX = Math.min(a[0], b[0]); - const maxX = Math.max(a[0], b[0]); - const minY = Math.min(a[1], b[1]); - const maxY = Math.max(a[1], b[1]); - const maxAbsLat = Math.max(Math.abs(minY), Math.abs(maxY)); - return { a, b, minX, maxX, minY, maxY, lonMetersPerDegree: SPHERE_METERS_PER_DEGREE * Math.cos(maxAbsLat*DEG_TO_RAD) }; -} - -/** - * Decomposes a geometry into its vertices and edges, then builds a search Shape. - * An R-tree index is created for shapes with more than 64 edges. - */ -function buildShape(geom: Geometry): Shape { - const vertices: Position[] = []; - const edges: Edge[] = []; - +/** Walks every edge of a geometry and rejects the ones crossing the antimeridian. */ +function assertNoAntimeridianCrossing(geom: Geometry) { function process(g: Geometry): void { switch (g.type) { - case "Point": vertices.push(g.coordinates); break; - case "MultiPoint": vertices.push(...g.coordinates); break; case "LineString": - vertices.push(...g.coordinates); - for (let i = 0; i < g.coordinates.length - 1; i++) edges.push(makeEdge(g.coordinates[i], g.coordinates[i + 1])); + for (let i = 0; i < g.coordinates.length - 1; i++) checkEdge(g.coordinates[i], g.coordinates[i + 1]); break; case "MultiLineString": for (const line of g.coordinates) { - vertices.push(...line); - for (let i = 0; i < line.length - 1; i++) edges.push(makeEdge(line[i], line[i + 1])); + for (let i = 0; i < line.length - 1; i++) checkEdge(line[i], line[i + 1]); } break; case "Polygon": for (const ring of g.coordinates) { - vertices.push(...ring); - for (let i = 0; i < ring.length - 1; i++) edges.push(makeEdge(ring[i], ring[i + 1])); + for (let i = 0; i < ring.length - 1; i++) checkEdge(ring[i], ring[i + 1]); } break; case "MultiPolygon": for (const poly of g.coordinates) for (const ring of poly) { - vertices.push(...ring); - for (let i = 0; i < ring.length - 1; i++) edges.push(makeEdge(ring[i], ring[i + 1])); + for (let i = 0; i < ring.length - 1; i++) checkEdge(ring[i], ring[i + 1]); } break; case "GeometryCollection": @@ -171,11 +132,7 @@ function buildShape(geom: Geometry): Shape { break; } } - process(geom); - const edgeIndex = edges.length > 64 // R-tree build cost is ~O(n log n); only worthwhile beyond this edge count - ? new RBush().load(edges) : null; - return { vertices, edges, edgeIndex }; } /** Interpolates a point at parameter t along a lon/lat segment. */ @@ -212,11 +169,10 @@ function nearestOnSegment(p: Position, a: Position, b: Position, pointDistance: let best = distanceAt(t); const segmentLengthMeters = pointDistance(a, b); const arcAngle = segmentLengthMeters / EARTH_RADIUS_M; - const requiredPrecisionMeters = Math.max(0.01, 0.005 * best); // stop when step < 0.5% of best distance, floored at 1 cm for ( let step = Math.min(0.5, 0.2 * arcAngle); // initial step: 0.2 of the arc angle (radians), capped at 0.5 t-units (half segment) - step * segmentLengthMeters > requiredPrecisionMeters; + step * segmentLengthMeters > TOLERANCE_M; // stop when step < 5 mm, since distance result is rounded to the cm step /= 2 ) { let improved = true; @@ -237,125 +193,251 @@ function nearestOnSegment(p: Position, a: Position, b: Position, pointDistance: return { distance: best, point: pointOnSegment(a, b, t) }; } -/** Computes a cheap metric lower bound from a point to an edge bbox. */ -function lowerBoundToEdge(p: Position, edge: Edge, pLonMetersPerDegree: number = SPHERE_METERS_PER_DEGREE*Math.cos(p[1]*DEG_TO_RAD)): number { - const latGap = Math.max(0, edge.minY - p[1], p[1] - edge.maxY); - const lonGap = Math.max(0, edge.minX - p[0], p[0] - edge.maxX); - const lonMetersPerDegree = Math.min(edge.lonMetersPerDegree, pLonMetersPerDegree); - return Math.max(latGap * MIN_LAT_METERS_PER_DEGREE, lonGap * lonMetersPerDegree); +/** + * Meridional and prime-vertical radii of curvature at `lat`, in meters. + * + * These give the projections below the same local scale as the metric the + * answer is finally reported in. That match is not cosmetic: candidate edges + * have to be *ranked* in the metric they are *measured* in, or a nearer edge + * can lose to a farther one. On WGS84 the meridional and parallel scales + * differ by ~0.67% at the equator, enough to flip the ranking of two edges + * 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 }; + const sinLat = Math.sin(lat * DEG_TO_RAD); + const w = 1 - WGS84_ECCENTRICITY_SQ * sinLat * sinLat; + return { + meridional: WGS84_SEMI_MAJOR_M * (1 - WGS84_ECCENTRICITY_SQ) / (w * Math.sqrt(w)), + primeVertical: WGS84_SEMI_MAJOR_M / Math.sqrt(w), + }; } /** - * Expands a metric radius around a point into a conservative lon/lat search box. - * - * latDegrees uses the same MIN_LAT_METERS_PER_DEGREE floor as lowerBoundToEdge, so - * it never underestimates the latitude span needed to reach `meters` away. - * - * lonDegrees solves the exact haversine equation for the worst-case latitude in - * [p.lat - latDegrees, p.lat + latDegrees] (the one closest to a pole, where a - * degree of longitude is shortest). The approximation can under-estimate - * the needed longitude span as soon as the box reaches into higher latitudes - * than p itself, letting the true nearest edge fall outside the search box. + * Builds the azimuthal equidistant projection centered at `center`. + * Returns `project` (lon/lat → [x,y] meters) and `unproject` ([x,y] meters → lon/lat). + * y points north, x points east from the center. */ -function safeSearchBox(p: Position, meters: number): Pick { - const latDegrees = meters / MIN_LAT_METERS_PER_DEGREE; - const maxAbsLatDeg = Math.min(89.999, Math.max(Math.abs(p[1] - latDegrees), Math.abs(p[1] + latDegrees))); - const cosProduct = Math.cos(p[1]*DEG_TO_RAD) * Math.cos(maxAbsLatDeg*DEG_TO_RAD); - - let lonDegrees = 180; - if (cosProduct > 1e-9) { - const halfAngle = meters / (2 * EARTH_RADIUS_M); - const sinSqHalfDeltaLambda = Math.sin(halfAngle) ** 2 / cosProduct; - if (sinSqHalfDeltaLambda >= 1) - lonDegrees = Math.min(180, (2 * Math.asin(Math.sqrt(sinSqHalfDeltaLambda))) / DEG_TO_RAD); +function makeAzimuthalEquidistant(center: Position, metric: Metric) { + const phi0 = center[1] * DEG_TO_RAD; + const lambda0 = center[0] * DEG_TO_RAD; + const sinPhi0 = Math.sin(phi0); + const cosPhi0 = Math.cos(phi0); + const { meridional, primeVertical } = curvatureRadii(center[1], metric); + + const reverseMap = new Map(); // shortcut to unproject for already `project`ed points + + /** + * 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. + */ + function radiusInDirection(east: number, north: number): number { + return 1 / (north * north / meridional + east * east / primeVertical); + } + + function project(p: Position): Position { + const phi = p[1] * DEG_TO_RAD; + const lambda = p[0] * DEG_TO_RAD; + const cosc = sinPhi0 * Math.sin(phi) + cosPhi0 * Math.cos(phi) * Math.cos(lambda - lambda0); + const c = Math.acos(Math.min(1, Math.max(-1, cosc))); + if (c < 1e-12) { + reverseMap.set("0,0", [p[0], p[1]]); + return [0, 0]; + } + const sinc = Math.sin(c); + const east = Math.cos(phi) * Math.sin(lambda - lambda0) / sinc; + const north = (cosPhi0 * Math.sin(phi) - sinPhi0 * Math.cos(phi) * Math.cos(lambda - lambda0)) / sinc; + const rho = radiusInDirection(east, north) * c; + const projected: Position = [rho * east, rho * north]; + reverseMap.set(`${projected[0]},${projected[1]}`, [p[0], p[1]]); + return projected; + } + + function unproject(p: Position): Position { + const [x, y] = p; + const rho = Math.sqrt(x * x + y * y); + if (rho < 1e-6) return center; + // (x, y) points the same way as (east, north), so the radius `project` + // scaled by is recoverable here and the round-trip stays exact. + const c = rho / radiusInDirection(x / rho, y / rho); + const sinC = Math.sin(c), cosC = Math.cos(c); + const phi = Math.asin(Math.min(1, Math.max(-1, cosC * sinPhi0 + y * sinC * cosPhi0 / rho))); + const lambda = lambda0 + Math.atan2(x * sinC, rho * cosPhi0 * cosC - y * sinPhi0 * sinC); + return [lambda / DEG_TO_RAD, phi / DEG_TO_RAD]; } + return { project, unproject, reverseMap }; +} + +/** + * Builds the equirectangular projection used for the initial planar seed: + * longitudes are scaled so that one x unit and one y unit span comparable + * ground distances around `lat0`. Raw lon/lat degrees are up to 1.4:1 + * anisotropic at French latitudes, which biases the planar nearest pair + * towards edges that are only nearer *in degrees* and costs the refinement + * loop a pass to recover. + * + * Both axes stay in degree-like units — only their ratio matters here, since + * every distance is remeasured with the real metric downstream. + */ +function makeEquirectangular(lat0: number, metric: Metric) { + const { meridional, primeVertical } = curvatureRadii(lat0, metric); + // Crushing longitudes towards a pole is the right answer, not a fallback: + // there a degree of longitude really does cover almost no ground. The floor + // only keeps `unproject` from dividing by zero exactly at ±90°. + const lonScale = Math.max(primeVertical * Math.cos(lat0 * DEG_TO_RAD) / meridional, 1e-6); return { - minX: p[0] - lonDegrees, - minY: p[1] - latDegrees, - maxX: p[0] + lonDegrees, - maxY: p[1] + latDegrees, + project: (p: Position): Position => [p[0] * lonScale, p[1]], + unproject: (p: Position): Position => [p[0] / lonScale, p[1]], }; } +/** + * Projects a GeoJSON Geometry through `project`, returning a new Geometry of + * the same type with projected coordinates. + */ +function projectGeometry(geom: Geometry, project: (_: Position) => Position): Geometry { + const projRing = (ring: Position[]) => ring.map(project); + switch (geom.type) { + case "Point": return { type: "Point", coordinates: project(geom.coordinates) }; + case "MultiPoint": return { type: "MultiPoint", coordinates: geom.coordinates.map(project) }; + case "LineString": return { type: "LineString", coordinates: geom.coordinates.map(project) }; + case "MultiLineString": return { type: "MultiLineString", coordinates: geom.coordinates.map(projRing) }; + case "Polygon": return { type: "Polygon", coordinates: geom.coordinates.map(projRing) }; + case "MultiPolygon": return { type: "MultiPolygon", coordinates: geom.coordinates.map(poly => poly.map(projRing)) }; + case "GeometryCollection": + return { type: "GeometryCollection", geometries: geom.geometries.map(g => projectGeometry(g, project)) }; + } +} -/** Finds the nearest point from p to a shape using vertex fallback and edge pruning. Returns null if nothing beats the upper bound. */ -function nearestOnShape(p: Position, shape: Shape, pointDistance: (a: Position, b: Position) => number, upperBound = Infinity): NearestPoint | null { - let best: NearestPoint | null = null; +/** + * Exact planar nearest pair between the *boundaries* of two JTS geometries, + * resolved through an STR-tree of facets. + * + * jsts' `DistanceOp` answers the same question with a nested loop over both + * segment lists, which is quadratic: for two 10 000-vertex polygons that is + * ~900 ms a pass against ~5 ms here. The trade is that facets carry no notion + * of interior, so this only finds the nearest pair once the geometries are + * already known not to intersect — see the `RelateOp.intersects` guard below. + */ +function planarNearestLocations(jstsA: JstsGeometry, jstsB: JstsGeometry): GeometryLocation[] { + return new IndexedFacetDistance(jstsA).nearestLocations(jstsB) as GeometryLocation[]; +} - if (shape.edges.length === 0) { - for (const v of shape.vertices) { - const d = pointDistance(p, v); - if (d < upperBound && (best === null || d < best.distance)) { - best = { distance: d, point: v }; - } +/** + * Turns a planar nearest pair into a geodesic one: the planar answer names the + * two facets involved, and the real distance is then minimized along them with + * `pointDistance`. + */ +function actualClosestOnGeometryLocation(locA: GeometryLocation, locB: GeometryLocation, pointDistance: (a: Position, b: Position) => number, unproject: (_: Position) => Position, reverseMap?: Map): DistanceResult { + function unproj(c: JstsCoord): Position { + return reverseMap?.get(`${c.x},${c.y}`) ?? unproject([c.x, c.y]); + } + + function getSegment(loc: GeometryLocation) { + const coords: JstsCoord[] = loc.getGeometryComponent().getCoordinates(); + if (coords.length > 1) { + const idx: number = loc.getSegmentIndex(); + return { start: unproj(coords[idx]), stop: unproj(coords[idx + 1]) }; } - return best; } - const pLonMetersPerDegree = SPHERE_METERS_PER_DEGREE * Math.cos(p[1]*DEG_TO_RAD); - const candidates = shape.edgeIndex - ? shape.edgeIndex.search(safeSearchBox(p, Math.min(upperBound, 20_000_000))) // 20 Mm ≈ half-Earth circumference: caps the bbox so safeSearchBox never overflows ±180°/±90° - : shape.edges; - - // Each edge first gets a cheap lower-bound test from its bbox in meters. - // Only the survivors pay for the more expensive point-to-segment refinement. - for (const edge of candidates) { - const bound = best === null ? upperBound : best.distance; - if (lowerBoundToEdge(p, edge, pLonMetersPerDegree) >= bound) continue; - const nearest = nearestOnSegment(p, edge.a, edge.b, pointDistance); - if (nearest.distance < (best === null ? upperBound : best.distance)) best = nearest; + // When the *query* side of an indexed search resolves to a bare Point facet, + // jsts tags that location with the base geometry's component and start index + // instead of its own (FacetSequence.nearestLocations, `isPointOther` branch), + // so both locations report the same component object. Only the coordinate is + // trustworthy there — which is all a vertex has anyway. + const mislabeledPointSide = locA.getGeometryComponent() === locB.getGeometryComponent(); + const segA = getSegment(locA); + const segB = mislabeledPointSide ? undefined : getSegment(locB); + const candidates = []; + if (segA) candidates.push({ point: unproj(locB.getCoordinate() as JstsCoord), seg: segA, rev: true }); + if (segB) candidates.push({ point: unproj(locA.getCoordinate() as JstsCoord), seg: segB, rev: false }); + if (candidates.length == 0) { // both sides are isolated points + const pA = unproj(locA.getCoordinate() as JstsCoord); + const pB = unproj(locB.getCoordinate() as JstsCoord); + return { distance: pointDistance(pA, pB), point1: pA, point2: pB }; } + let best: DistanceResult = { distance: Infinity, point1: [0, 0], point2: [0, 0] }; + for (const { point, seg, rev } of candidates) { + const nearest = nearestOnSegment(point, seg.start, seg.stop, pointDistance); + if (nearest.distance < best.distance) { + const [point1, point2] = rev ? [nearest.point, point] : [point, nearest.point]; + best = { distance: nearest.distance, point1, point2 }; + } + } return best; } /** * Compute geodesic distance in meters between gA and gB, with closest points. * - * JTS gives us a robust planar nearest-point seed and detects all exact touches, - * overlaps and containments. If the planar answer is not already zero, the final - * geodesic result is refined by scanning both directions: vertices of A against - * shape B, then vertices of B against shape A. - * * @param {object} gA GeoJSON Geometry * @param {object} gB GeoJSON Geometry * @param metric Point-to-point metric: "haversine" (default) or "vincenty". */ -export function distance(gA: Geometry, gB: Geometry, metric: "haversine" | "vincenty" = "haversine"): DistanceResult { +export function distance(gA: Geometry, gB: Geometry, metric: Metric = "haversine"): DistanceResult { const pointDistance = metric === "vincenty" ? distanceVincenty : haversine; - - assertNoAntimeridian(gA); - assertNoAntimeridian(gB); - - const shapeA = buildShape(gA); - const shapeB = buildShape(gB); - assertNoAntimeridianPair(shapeA.vertices, shapeB.vertices); - - const reader = new (GeoJSONReader as unknown as new () => { read(g: Geometry): unknown })(); - type JstsCoord = { x: number; y: number }; - const [n1, n2] = DistanceOp.nearestPoints(reader.read(gA), reader.read(gB)) as [JstsCoord, JstsCoord]; - const p1: Position = [n1.x, n1.y]; - const p2: Position = [n2.x, n2.y]; - - if (Math.hypot(n1.x - n2.x, n1.y - n2.y) < 1e-9) { // ~0.1 mm in degrees: below JTS floating-point noise → exact touch - return { distance: 0, point1: p1, point2: p1 }; + 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 }; + + assertNoAntimeridianCrossing(gA); + assertNoAntimeridianCrossing(gB); + + const reader = new GeoJSONReader(new GeometryFactory()); + const jstsA = reader.read(gA); + const jstsB = reader.read(gB); + const envA = jstsA.getEnvelopeInternal(); + const envB = jstsB.getEnvelopeInternal(); + + // Check whether the distance is 0. + // If an intersection is detected via the cheaper `RelateOp.intersects`, + // use the costlier DistanceOp to retrieve the intersection point. + if (envA.distance(envB) === 0 && RelateOp.intersects(jstsA, jstsB)) { + const [n1] = new DistanceOp(jstsA, jstsB).nearestPoints(); + return { distance: 0, point1: [n1.x, n1.y], point2: [n1.x, n1.y] }; } - let best: DistanceResult = { distance: pointDistance(p1, p2), point1: p1, point2: p2 }; - - for (const v of shapeA.vertices) { - const nearest = nearestOnShape(v, shapeB, pointDistance, best.distance); - if (nearest && nearest.distance < best.distance) { - best = { distance: nearest.distance, point1: v, point2: nearest.point }; + // Seed on the exact planar nearest pair. `actualClosestOnGeometryLocation` + // already minimizes over both facets of that pair, so the mirrored ordering + // yields the same distance and needs no second evaluation. + const seedLat = (Math.min(envA.getMinY(), envB.getMinY()) + Math.max(envA.getMaxY(), envB.getMaxY())) / 2; + const seed = makeEquirectangular(seedLat, metric); + const seedLocations = planarNearestLocations( + reader.read(projectGeometry(gA, seed.project)), + reader.read(projectGeometry(gB, seed.project)), + ); + let best = actualClosestOnGeometryLocation(seedLocations[0], seedLocations[1], pointDistance, seed.unproject); + + // Re-run the planar search in an equidistant azimuthal projection centered on + // the current best pair, keeping whichever pair measures shorter. + let converged = false; + for (let i = 0; i < MAX_REFINEMENT_PASSES && !converged; i++) { + const center = midpoint(best.point1, best.point2).geometry.coordinates; + const { project, unproject, reverseMap } = makeAzimuthalEquidistant(center, metric); + const projA = projectGeometry(gA, project); + const projB = projectGeometry(gB, project); + const [newLocA, newLocB] = planarNearestLocations(reader.read(projA), reader.read(projB)); + const newBest = actualClosestOnGeometryLocation(newLocA, newLocB, pointDistance, unproject, reverseMap); + if (newBest.distance > best.distance + TOLERANCE_M) { + // This projection has nothing better to offer, so stop and keep the + // incumbent. Not a fixed point, despite the flag: unlike the seed's + // equirectangular plane, which is affine in lon/lat and therefore maps + // edges to edges, a straight line in the azimuthal plane is a geodesic + // rather than the lon/lat-linear edge it stands for. On edges spanning + // thousands of kilometres the two curves separate far enough that the + // planar pair can name a point off the geometry, and this branch is what + // keeps such a pair from being adopted. + converged = true; + break; } + converged = best.distance - newBest.distance < TOLERANCE_M; + best = newBest; } - - for (const v of shapeB.vertices) { - const nearest = nearestOnShape(v, shapeA, pointDistance, best.distance); - if (nearest && nearest.distance < best.distance) { - best = { distance: nearest.distance, point1: nearest.point, point2: v }; - } + if (!converged) { + throw new Error("Convergence error in the distance algorithm: cannot compute the distance."); } return { ...best, distance: Math.round(best.distance * 100) / 100 }; diff --git a/test/helpers/distance.test.ts b/test/helpers/distance.test.ts index 124f9853..5e3e5e75 100644 --- a/test/helpers/distance.test.ts +++ b/test/helpers/distance.test.ts @@ -152,6 +152,65 @@ describe("distance helper", () => { expectZeroBothWays(gcA, gcB, "intersecting geometry collections"); }); + // Zero distance is settled in two steps: an indexed pass over the + // boundaries, then a point-in-area test per component for the strict + // containment the boundaries cannot show. These cases pin both outcomes of + // that second step. + it("returns zero for a polygon strictly inside another polygon", () => { + const outer: Polygon = { type: "Polygon", coordinates: [[[0, 0], [10, 0], [10, 10], [0, 10], [0, 0]]] }; + const inner: Polygon = { type: "Polygon", coordinates: [[[4, 4], [6, 4], [6, 6], [4, 6], [4, 4]]] }; + expectZeroBothWays(inner, outer, "polygon strictly inside a polygon"); + }); + + it("returns zero for a line strictly inside a polygon", () => { + const polygon: Polygon = { type: "Polygon", coordinates: [[[0, 0], [10, 0], [10, 10], [0, 10], [0, 0]]] }; + const line: LineString = { type: "LineString", coordinates: [[4, 4], [6, 6]] }; + expectZeroBothWays(line, polygon, "line strictly inside a polygon"); + }); + + it("returns zero when only one part of a MultiPolygon is inside", () => { + const outer: Polygon = { type: "Polygon", coordinates: [[[0, 0], [10, 0], [10, 10], [0, 10], [0, 0]]] }; + const parts: MultiPolygon = { + type: "MultiPolygon", + coordinates: [ + [[[4, 4], [6, 4], [6, 6], [4, 6], [4, 4]]], // inside + [[[40, 40], [41, 40], [41, 41], [40, 41], [40, 40]]], // far outside + ], + }; + expectZeroBothWays(parts, outer, "one MultiPolygon part inside"); + }); + + it("returns positive for a polygon sitting in another polygon's hole", () => { + const ringWithHole: Polygon = { + type: "Polygon", + coordinates: [ + [[0, 0], [10, 0], [10, 10], [0, 10], [0, 0]], + [[3, 3], [7, 3], [7, 7], [3, 7], [3, 3]], + ], + }; + const inHole: Polygon = { type: "Polygon", coordinates: [[[4, 4], [6, 4], [6, 6], [4, 6], [4, 4]]] }; + const d = ensureSymmetricDistance(inHole, ringWithHole, "polygon in a hole"); + // Nearest is the 1 degree gap from the square's edge to the hole's wall. + expect(d).toBeGreaterThan(0); + expectCloseRatio(d, haversine([6, 6], [7, 6]), 0.01, "polygon in a hole"); + }); + + it("returns positive for a square parked in a C-shape's notch", () => { + // The C-shape's bounding box swallows the square's whole bounding box, so + // the envelope shortcut cannot rule out an intersection and the + // containment test has to answer "not contained" on its own. At + // y in (1, 3) the C-shape's interior is only x in (0, 1). + const cShape: Polygon = { + type: "Polygon", + coordinates: [[[0, 0], [10, 0], [10, 1], [1, 1], [1, 3], [10, 3], [10, 4], [0, 4], [0, 0]]], + }; + const inNotch: Polygon = { type: "Polygon", coordinates: [[[2, 1.2], [9, 1.2], [9, 2.8], [2, 2.8], [2, 1.2]]] }; + const d = ensureSymmetricDistance(inNotch, cShape, "square in a C-shape notch"); + expect(d).toBeGreaterThan(0); + // Nearest is the 0.2 degree gap down to the C-shape's y = 1 edge. + expectCloseRatio(d, haversine([5, 1], [5, 1.2]), 0.01, "square in a C-shape notch"); + }); + it("does not misclassify outside point for polygon with duplicate vertex", () => { const polygonWithDuplicateVertex: Polygon = { type: "Polygon", @@ -166,11 +225,26 @@ describe("distance helper", () => { describe("antimeridian and poles", () => { // Antimeridian-crossing geometries are rejected by design (RFC 7946 SHOULD split them). - const capVertexCount = 36; - const ring = Array.from({ length: capVertexCount }, (_, i) => - [-180 + (i * 360) / capVertexCount, 80] as [number, number], - ); - ring.push(ring[0]); + // + // The rejection is load-bearing, not merely conservative. Two steps of the + // pipeline read raw lon/lat degrees as a plane: the zero-distance test + // (`RelateOp.intersects`) and the equirectangular seed. A ring with a + // |Δlon| > 180° edge is a different shape in that plane — it wraps the long + // way round the globe — so both steps answer about the wrong geometry: + // - containment goes undetected, because the planar interior is the + // complement of the intended one (a point inside a 170°..190° polygon + // is reported ~220 km away from it instead of at distance 0), + // - worse, the planar phantom of a crossing LineString sweeps the + // hemisphere it never visits, so `intersects` fires against unrelated + // geometry: a 165°..-174.5° line is reported as intersecting a box at + // lon 23°..28°, i.e. distance 0 instead of ~13 000 km. + // Only the azimuthal refinement is antimeridian-safe, since it projects each + // vertex spherically; it cannot rescue either step above. + const polarRing = (lat: number, sign = 1) => { + const ring = Array.from({ length: 36 }, (_, i) => [sign * (-180 + (i * 360) / 36), lat] as [number, number]); + ring.push(ring[0]); + return ring; + }; it.each([ { @@ -183,27 +257,98 @@ describe("distance helper", () => { [[175, 3], [175, 7], [178, 7], [178, 3], [175, 3]], ], } as Geometry, - error: /RFC 7946/, }, { name: "throws for a polar-cap polygon whose closing edge spans > 180° of longitude", g1: { type: "Point", coordinates: [0, 85] } as Geometry, - g2: { type: "Polygon", coordinates: [ring] } as Geometry, - error: /RFC 7946/, + g2: { type: "Polygon", coordinates: [polarRing(80)] } as Geometry, }, { - name: "throws for a pair of individually valid geometries facing each other across the antimeridian", - g1: { + name: "throws for a south polar cap", + g1: { type: "Point", coordinates: [0, -85] } as Geometry, + g2: { type: "Polygon", coordinates: [polarRing(-80)] } as Geometry, + }, + { + // Both rings encircle a pole, so both carry a > 180° closing edge. The + // guard has to walk holes, not just shells. + name: "throws for a polar annulus whose shell and hole both encircle the pole", + g1: { type: "Point", coordinates: [0, 80] } as Geometry, + g2: { type: "Polygon", coordinates: [polarRing(70), polarRing(85, -1)] } as Geometry, + }, + { + // Shell stays east of the antimeridian; only the hole crosses it. + name: "throws for a polygon whose hole alone crosses the antimeridian", + g1: { type: "Point", coordinates: [0, 0] } as Geometry, + g2: { type: "Polygon", - coordinates: [[[179.8, -0.1], [179.9, -0.1], [179.9, 0.1], [179.8, 0.1], [179.8, -0.1]]], + coordinates: [ + [[160, -10], [160, 20], [170, 20], [170, -10], [160, -10]], + [[175, 0], [175, 10], [-175, 10], [-175, 0], [175, 0]], + ], } as Geometry, + }, + { + name: "throws when only one MultiPolygon part crosses the antimeridian", + g1: { type: "Point", coordinates: [50, 50] } as Geometry, g2: { - type: "Polygon", - coordinates: [[[-179.9, -0.1], [-179.8, -0.1], [-179.8, 0.1], [-179.9, 0.1], [-179.9, -0.1]]], + type: "MultiPolygon", + coordinates: [ + [[[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]]], + [[[170, 0], [170, 10], [-170, 10], [-170, 0], [170, 0]]], + ], } as Geometry, - error: /Antimeridian-adjacent/, }, - ])("$name", ({ g1, g2, error }) => expectThrowBothWays(g1, g2, error)); + { + name: "throws when only one MultiLineString member crosses the antimeridian", + g1: { type: "Point", coordinates: [50, 50] } as Geometry, + g2: { type: "MultiLineString", coordinates: [[[0, 0], [1, 1]], [[170, 5], [-170, 5]]] } as Geometry, + }, + { + name: "throws when a GeometryCollection member crosses the antimeridian", + g1: { type: "Point", coordinates: [50, 50] } as Geometry, + g2: { + type: "GeometryCollection", + geometries: [ + { type: "Point", coordinates: [0, 0] }, + { type: "LineString", coordinates: [[170, 5], [-170, 5]] }, + ], + } as Geometry, + }, + ])("$name", ({ g1, g2 }) => expectThrowBothWays(g1, g2, /RFC 7946/)); + + // The guard triggers on |Δlon| > 180°, so a half-globe edge is still legal + // and must not be swept up with the crossing ones. + it("accepts an edge spanning exactly 180° of longitude", () => { + const line: Geometry = { type: "LineString", coordinates: [[-90, 10], [90, 10]] }; + const point: Geometry = { type: "Point", coordinates: [0, 0] }; + // Edges interpolate linearly in lon/lat, so the line passes through [0,10]. + expectCloseRatio( + ensureSymmetricDistance(point, line, "exactly 180 degrees of longitude"), + haversine([0, 0], [0, 10]), + 0.001, + "exactly 180 degrees of longitude", + ); + }); + + // The following case used to expect an error, but the new implementation + // supports antimeridian-adjacent geometries. Test it as a normal symmetric + // distance case instead. + it("accepts a pair of individually valid geometries facing each other across the antimeridian", () => { + const g1: Geometry = { + type: "Polygon", + coordinates: [[[179.8, -0.1], [179.9, -0.1], [179.9, 0.1], [179.8, 0.1], [179.8, -0.1]]], + } as Geometry; + const g2: Geometry = { + type: "Polygon", + coordinates: [[[-179.9, -0.1], [-179.8, -0.1], [-179.8, 0.1], [-179.9, 0.1], [-179.9, -0.1]]], + } as Geometry; + const d = ensureSymmetricDistance(g1, g2, "antimeridian-adjacent polygons"); + // Expected distance ~ haversine between 179.9 and -179.9 at equator (~22.3 km). + const expected = Math.round(haversine([179.9, 0], [-179.9, 0]) * 100) / 100; + // Check to the nearest decade of kilometers (±10 km) and sanity bounds. + expect(Math.abs(d - expected), "antimeridian distance close to expected").toBeLessThanOrEqual(100); + expect(Math.abs(22_300 - expected), "antimeridian distance close to expected").toBeLessThanOrEqual(100); + }); }); describe("antimeridian-safe Multi* geometries (RFC 7946 conforming)", () => { @@ -497,7 +642,10 @@ describe("distance helper", () => { }, ])("$name", ({ g1, g2 }) => expectThrowBothWays(g1, g2, /RFC 7946/)); - it("does not follow a naive s1/s2 result", () => { + // Two long sub-antarctic lines. Read as a plane of raw degrees they look + // ~2 440 km apart; the true nearest pair is ~986 km, on the far end of the + // first line. Pins that the search is not fooled by degree-space geometry. + it("picks the true nearest pair on two long sub-antarctic lines", () => { const s1: Geometry = { type: "LineString", coordinates: [ @@ -514,30 +662,36 @@ describe("distance helper", () => { }; const result = distance(s1, s2); - const d = ensureSymmetricDistance(s1, s2, "s1/s2 distance"); - expectCloseRatio(d, 986_440, 0.005, "s1/s2 distance"); + const d = ensureSymmetricDistance(s1, s2, "long sub-antarctic lines"); + expectCloseRatio(d, 986_440, 0.005, "long sub-antarctic lines"); expect(d).toBeLessThan(2_440_829.44); - expectCloseRatio(result.point1[0], -85.7, 0.02, "s1/s2 point1 lon"); - expectCloseRatio(result.point1[1], -66.77, 0.02, "s1/s2 point1 lat"); - expectCloseRatio(result.point2[0], -86.16447630165361, 0.001, "s1/s2 point2 lon"); - expectCloseRatio(result.point2[1], -57.89895966720661, 0.001, "s1/s2 point2 lat"); + expectCloseRatio(result.point1[0], -85.7, 0.02, "long sub-antarctic lines point1 lon"); + expectCloseRatio(result.point1[1], -66.77, 0.02, "long sub-antarctic lines point1 lat"); + expectCloseRatio(result.point2[0], -86.16447630165361, 0.001, "long sub-antarctic lines point2 lon"); + expectCloseRatio(result.point2[1], -57.89895966720661, 0.001, "long sub-antarctic lines point2 lat"); }); - it("finds nearest edge with Vincenty when lower bounds use haversine METERS_PER_DEGREE", () => { - const p: Geometry = { type: "Point", coordinates: [0, 0] }; + // Both edges span tens of degrees, so the linear lon/lat path they stand + // for and the geodesic the azimuthal plane refines along separate enough + // that the refinement's candidate pair drifts hundreds of km off the real + // edges (see the loop in `distance` for why that guard exists at all). + // Caught here as a thrown error rather than a silently wrong distance. + it("throws instead of returning a wrong distance for two very long, far-apart lines", () => { + const a: Geometry = { + type: "LineString", + coordinates: [ + [73.95631313323975, 38.032363414764404], + [41.048042762817566, -35.91032814939557], + ], + }; const b: Geometry = { - type: "MultiLineString", + type: "LineString", coordinates: [ - [[0.9963, -1], [0.9963, 1]], // planar seed ~110907 m, lon-edge - [[-1, 1.0], [1, 1.0]], // true nearest ~110574 m, lat-edge + [-33.47427845001221, 15.573780059814453], + [-38.53770555856234, -89], ], }; - expectCloseRatio( - ensureSymmetricDistance(p, b, "Vincenty avoids over-pruning lat bounds", "vincenty"), - 110574.39, - 0.001, - "Vincenty avoids over-pruning lat bounds", - ); + expectThrowBothWays(a, b, /Convergence error/); }); it("regresses the concrete RBush false-pruning case with a poleward nearest edge", () => { @@ -597,8 +751,67 @@ 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 + // 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. + describe("ranks candidate edges in the requested metric", () => { + // A meridian edge dLon away, and a parallel edge 1 degree north. Ranking + // them one way or the other is exactly the sphere/WGS84 disagreement: + // - a sphere prefers the meridian edge while dLon < 1 / cos(lat), + // - WGS84 prefers the parallel edge once dLon > M / (N cos(lat)), + // and M / N < 1 always, so every dLon in between is ranked differently by + // the two metrics. Each dLon here is the midpoint of that window. + const cases = [ + { lat: 0, dLon: 0.996653 }, // window (0.993306, 1.000000) + { lat: 45, dLon: 1.411839 }, // window (1.409464, 1.414214) + { lat: 60, dLon: 1.998318 }, // window (1.996636, 2.000000) + ]; + + it.each(cases)("prefers the parallel edge under Vincenty at lat $lat", ({ lat, dLon }) => { + const p: Geometry = { type: "Point", coordinates: [0, lat] }; + const edges: Geometry = { + type: "MultiLineString", + coordinates: [ + [[dLon, lat - 1], [dLon, lat + 1]], // meridian edge: nearer on a sphere + [[-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]); + 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"); + 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 }) => { + const p: Geometry = { type: "Point", coordinates: [0, lat] }; + const edges: Geometry = { + type: "MultiLineString", + coordinates: [ + [[dLon, lat - 1], [dLon, lat + 1]], + [[-dLon, lat + 1], [dLon, lat + 1]], + ], + }; + const parallelEdge = haversine([0, lat], [0, lat + 1]); + const meridianEdge = haversine([0, lat], [dLon, lat]); + expect(meridianEdge, `lat ${lat}: on a sphere the meridian edge must be the nearer one`).toBeLessThan(parallelEdge); + + // 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"); + 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}`); + }); + }); + describe("performance sanity", () => { - it("resolves two disjoint 500-vertex polygons in under a second", () => { + it("resolves two disjoint 500-vertex polygons in under 10 millisecond", () => { const circle = (centerLon: number, centerLat: number, radiusDeg: number, n: number): Geometry => ({ type: "Polygon", coordinates: [ @@ -617,10 +830,10 @@ describe("distance helper", () => { const elapsedMs = performance.now() - start; expect(result.distance).toBeGreaterThan(0); - expect(elapsedMs).toBeLessThan(1_000); + expect(elapsedMs).toBeLessThan(10); }); - it("resolves two nearby irregular ~5000-vertex polygons in under a second", () => { + it("resolves two nearby irregular ~10000-vertex polygons in under a second", () => { const blob = (centerLon: number, centerLat: number, radiusDeg: number, n: number, seed: number): Geometry => ({ type: "Polygon", coordinates: [ @@ -632,8 +845,8 @@ describe("distance helper", () => { ], }); - const a = blob(2.0, 48.0, 0.5, 5000, 1); - const b = blob(4.0, 48.3, 0.5, 5000, 2); + const a = blob(2.0, 48.0, 0.5, 10000, 1); + const b = blob(4.0, 48.3, 0.5, 10000, 2); const start = performance.now(); const result = distance(a, b); From c12fc3dffcff5def382b0949d1f89660e1cc3a89 Mon Sep 17 00:00:00 2001 From: Lionel Zoubritzky Date: Mon, 3 Aug 2026 13:26:50 +0200 Subject: [PATCH 03/11] feat: Add distance tool --- docs/mcp-tools.md | 130 ++++++++++++++++++++++++++++++++ src/tools/DistanceTool.ts | 93 +++++++++++++++++++++++ test/integration/samples.ts | 1 + test/tools/distance.test.ts | 52 +++++++++++++ test/tools/strict-input.test.ts | 9 +++ 5 files changed, 285 insertions(+) create mode 100644 src/tools/DistanceTool.ts create mode 100644 test/tools/distance.test.ts diff --git a/docs/mcp-tools.md b/docs/mcp-tools.md index f8c9cf24..16369915 100644 --- a/docs/mcp-tools.md +++ b/docs/mcp-tools.md @@ -52,6 +52,7 @@ Annotations MCP exposées dans la définition `tools/list` de chaque tool : - [`gpf_count_features`](#gpf_count_features) - [`gpf_get_feature_by_id`](#gpf_get_feature_by_id) - [`gpf_get_feature_by_id_layer`](#gpf_get_feature_by_id_layer) +- [`distance`](#distance) ## `geocode` @@ -2212,3 +2213,132 @@ Cet outil ne peut renvoyer qu'un unique objet (0 ou plusieurs résultats provoqu | --- | --- | --- | --- | | Succès | oui | oui | `content[0].text` est `JSON.stringify(structuredContent)`. | | Erreur | oui | non | `content[0].text` porte le message d'erreur ; aucun `structuredContent` n'est ajouté (réservé au `outputSchema` du cas de succès). | + +## `distance` + +Code Source : [src/tools/DistanceTool.ts](../src/tools/DistanceTool.ts) + +### Titre + +Distance entre deux points + +### Description du tool + +``` +Renvoie la distance (en mètres) entre deux points à partir de leur longitude et latitude. +``` + +### Schéma d’entrée + +| Champ | Type | Requis | Description | +| --- | --- | --- | --- | +| `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. | + +

+Schéma d’entrée brut + +```json +{ + "type": "object", + "properties": { + "departure": { + "type": "object", + "properties": { + "lon": { + "type": "number", + "minimum": -180, + "maximum": 180, + "description": "La longitude du point de départ." + }, + "lat": { + "type": "number", + "minimum": -90, + "maximum": 90, + "description": "La latitude du point de départ." + } + }, + "required": [ + "lon", + "lat" + ], + "additionalProperties": false, + "description": "Le point de départ" + }, + "arrival": { + "type": "object", + "properties": { + "lon": { + "type": "number", + "minimum": -180, + "maximum": 180, + "description": "La longitude du point d'arrivée." + }, + "lat": { + "type": "number", + "minimum": -90, + "maximum": 90, + "description": "La latitude du point d'arrivée." + } + }, + "required": [ + "lon", + "lat" + ], + "additionalProperties": false, + "description": "Le point d'arrivée" + }, + "profile": { + "type": "string", + "enum": [ + "direct", + "vincenty" + ], + "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`." + } + }, + "required": [ + "departure", + "arrival" + ], + "additionalProperties": false, + "$schema": "http://json-schema.org/draft-07/schema#" +} +``` + +
+ +### Schéma de sortie + +| Champ | Type | Requis | Description | +| --- | --- | --- | --- | +| `distance` | number | oui | La distance entre les deux points, en mètres. | + +
+Schéma de sortie brut + +```json +{ + "type": "object", + "properties": { + "distance": { + "type": "number", + "description": "La distance entre les deux points, en mètres." + } + }, + "required": [ + "distance" + ] +} +``` + +
+ +### Réponse MCP + +| Cas | `content` | `structuredContent` | Relation entre `content` et `structuredContent` | +| --- | --- | --- | --- | +| Succès | oui | oui | `content[0].text` est `JSON.stringify(structuredContent)`. | +| Erreur | oui | non | `content[0].text` porte le message d'erreur ; aucun `structuredContent` n'est ajouté (réservé au `outputSchema` du cas de succès). | diff --git a/src/tools/DistanceTool.ts b/src/tools/DistanceTool.ts new file mode 100644 index 00000000..f8e9e209 --- /dev/null +++ b/src/tools/DistanceTool.ts @@ -0,0 +1,93 @@ +/** + * MCP tool exposing distance lookup for a single geographic position. + */ + +import BaseTool from "./BaseTool.js"; +import { z } from "zod"; + +import { READ_ONLY_OPEN_WORLD_TOOL_ANNOTATIONS } from "../helpers/toolAnnotations.js"; +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"; + +// --- Schemas --- + +const distanceInputSchema = z.object({ + departure: z.object({ + lon: lonSchema.describe("La longitude du point de départ."), + lat: latSchema.describe("La latitude du point de départ."), + }).describe("Le point de départ"), + arrival: z.object({ + lon: lonSchema.describe("La longitude du point d'arrivée."), + lat: latSchema.describe("La latitude du point d'arrivée."), + }).describe("Le point d'arrivée"), + profile: z + .enum(["direct", "vincenty"]) + .default("direct") + .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`." + ].join("")), +}).strict(); + +const distanceOutputSchema = z.object({ + distance: z.number().describe("La distance entre les deux points, en mètres."), +}); + +// --- Types --- + +type DistanceInput = z.infer; + +// --- Tool --- + +const DISTANCE_TOOL_DESCRIPTION = `Renvoie la distance (en mètres) entre deux points à partir de leur longitude et latitude.`; + +class DistanceTool extends BaseTool { + name = "distance"; + title = "Distance entre deux points"; + annotations = READ_ONLY_OPEN_WORLD_TOOL_ANNOTATIONS; + description = DISTANCE_TOOL_DESCRIPTION; + protected outputSchemaShape = distanceOutputSchema; + + schema = distanceInputSchema; + + // The framework requires a plain Zod object here to publish a compatible + // input schema, but its `.default()` fields are (incorrectly) marked required. + get inputSchema() { + return generatePublishedInputSchema(distanceInputSchema); + } + + /** + * Resolves the distance query. + * + * @param input Normalized tool input. + * @returns The distance. + */ + async execute(input: DistanceInput) { + logger.info(`[tool] execute ${this.name} ...`, { + input: input + }); + + switch (input.profile) { + case "direct": + case "vincenty": { + const pointDistance = input.profile == "direct" ? haversine : distanceVincenty; + const raw = pointDistance( + [input.departure.lon, input.departure.lat], + [input.arrival.lon, input.arrival.lat] + ); + return { + distance: Math.round(raw * 100) / 100 + }; + } + default: { + const profile: never = input.profile; + throw new Error(`Impossible profile ${profile}`); + } + } + } +} + +export default DistanceTool; diff --git a/test/integration/samples.ts b/test/integration/samples.ts index 2107b45c..e85668c5 100644 --- a/test/integration/samples.ts +++ b/test/integration/samples.ts @@ -19,6 +19,7 @@ export const EXPECTED_TOOL_NAMES = [ "cadastre", "urbanisme", "assiette_sup", + "distance", "gpf_search_types", "gpf_describe_type", "gpf_get_features", diff --git a/test/tools/distance.test.ts b/test/tools/distance.test.ts new file mode 100644 index 00000000..fa749949 --- /dev/null +++ b/test/tools/distance.test.ts @@ -0,0 +1,52 @@ +import { describe, it, expect } from "vitest"; + +import DistanceTool from "../../src/tools/DistanceTool"; +import { validateStructuredContentAgainstOutputSchema } from "./helpers/outputSchema"; +import { expectErrorText } from "./helpers/errorAssertions"; + +describe("Test DistanceTool", () => { + const departure = { lon: 2.3522, lat: 48.8566 }; + const arrival = { lon: 2.2945, lat: 48.8584 }; + + it("should publish an optional profile and a distance output schema", () => { + const tool = new 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", + }); + expect(tool.toolDefinition.outputSchema).toBeDefined(); + }); + + it.each([undefined, "direct", "vincenty"])("should return a structured distance for profile %s", async (profile) => { + const tool = new DistanceTool(); + const response = await tool.toolCall({ + params: { + name: "distance", + arguments: { departure, arrival, ...(profile && { profile }) }, + }, + }); + + expect(response.isError).toBeUndefined(); + expect(response.structuredContent).toMatchObject({ distance: expect.any(Number) }); + expect((response.structuredContent as { distance: number }).distance).toBeGreaterThan(0); + expect(response.content[0]).toMatchObject({ type: "text" }); + expect(validateStructuredContentAgainstOutputSchema( + tool.toolDefinition.outputSchema, + response.structuredContent, + )).toBeNull(); + }); + + it("should reject invalid coordinates at the tool boundary", async () => { + const tool = new DistanceTool(); + const response = await tool.toolCall({ + params: { + name: "distance", + arguments: { departure: { lon: 600, lat: 48.8566 }, arrival }, + }, + }); + + expect(expectErrorText(response)).toContain("departure.lon: La valeur doit être au plus 180."); + }); +}); diff --git a/test/tools/strict-input.test.ts b/test/tools/strict-input.test.ts index f9398bff..284fb450 100644 --- a/test/tools/strict-input.test.ts +++ b/test/tools/strict-input.test.ts @@ -4,6 +4,7 @@ import AdminexpressTool from "../../src/tools/AdminexpressTool"; import AltitudeTool from "../../src/tools/AltitudeTool"; import AssietteSupTool from "../../src/tools/AssietteSupTool"; import CadastreTool from "../../src/tools/CadastreTool"; +import DistanceTool from "../../src/tools/DistanceTool"; import GeocodeTool from "../../src/tools/GeocodeTool"; import GpfCountFeaturesTool from "../../src/tools/GpfCountFeaturesTool"; import GpfDescribeTypeTool from "../../src/tools/GpfDescribeTypeTool"; @@ -33,6 +34,14 @@ const strictInputCases = [ tool: new CadastreTool(), validArguments: { lon: 2.3522, lat: 48.8566 }, }, + { + label: "DistanceTool", + tool: new DistanceTool(), + validArguments: { + departure: { lon: 2.3522, lat: 48.8566 }, + arrival: { lon: 2.2945, lat: 48.8584 }, + }, + }, { label: "GeocodeTool", tool: new GeocodeTool(), From d3382a2ee16aed80b64f3720694bbf4f6ad6ac28 Mon Sep 17 00:00:00 2001 From: esgn <5435148+esgn@users.noreply.github.com> Date: Thu, 1 Oct 2026 17:43:39 +0200 Subject: [PATCH 04/11] fix(distance): replace node-vincenty with geographiclib-geodesic node-vincenty (0.0.6, unmaintained, untyped) returns an undefined distance for antipodal or nearly antipodal points, which the wrapper let through. Karney's algorithm always converges; over 40,000 random pairs it differs from Vincenty by at most 0.56 mm. The hand-written node-vincenty types and their tsconfig.test.json entry go away. --- package-lock.json | 14 +++++++------- package.json | 4 ++-- src/helpers/distance.ts | 18 +++++++++--------- src/types/node-vincenty.d.ts | 15 --------------- test/helpers/distance.test.ts | 6 ++++++ tsconfig.test.json | 1 - 6 files changed, 24 insertions(+), 34 deletions(-) delete mode 100644 src/types/node-vincenty.d.ts diff --git a/package-lock.json b/package-lock.json index 64eca969..cb6bc530 100644 --- a/package-lock.json +++ b/package-lock.json @@ -21,10 +21,10 @@ "@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 +3414,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 +4723,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..dc26f79a 100644 --- a/package.json +++ b/package.json @@ -65,12 +65,12 @@ "@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..20a9ebfe 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"; @@ -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 `distanceVincenty` 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); @@ -443,14 +443,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. */ +/** + * 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 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; + // 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/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..b8ec79cd 100644 --- a/test/helpers/distance.test.ts +++ b/test/helpers/distance.test.ts @@ -53,6 +53,12 @@ describe("distance helper", () => { 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(distanceVincenty(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"); 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": [ From 2931922e778e7221f988d64c17959e18b1fffd1f Mon Sep 17 00:00:00 2001 From: esgn <5435148+esgn@users.noreply.github.com> Date: Fri, 2 Oct 2026 10:32:22 +0200 Subject: [PATCH 05/11] refactor(distance): name distances after the Earth model, not Vincenty The distance tool profiles `direct` and `vincenty` become `spherical` (still the default) and `ellipsoidal`, and the helper's `distanceVincenty` and its `"haversine"`/`"vincenty"` metrics become `ellipsoidalDistance` and `"spherical"`/`"ellipsoidal"`: since the switch to geographiclib, nothing runs Vincenty's formula anymore. --- docs/mcp-tools.md | 10 +++++----- src/helpers/distance.ts | 16 ++++++++-------- src/tools/DistanceTool.ts | 18 +++++++++--------- test/helpers/distance.test.ts | 30 +++++++++++++++--------------- test/tools/distance.test.ts | 6 +++--- 5 files changed, 40 insertions(+), 40 deletions(-) diff --git a/docs/mcp-tools.md b/docs/mcp-tools.md index 16369915..94ff8dd9 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 et coûteuse, précision à 1mm). Par défaut : `spherical`. 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 et coûteuse, précision à 1mm). Par défaut : `spherical`." } }, "required": [ diff --git a/src/helpers/distance.ts b/src/helpers/distance.ts index 20a9ebfe..90068e0b 100644 --- a/src/helpers/distance.ts +++ b/src/helpers/distance.ts @@ -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 `distanceVincenty` 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); @@ -376,10 +376,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 }; @@ -448,7 +448,7 @@ export function distance(gA: Geometry, gB: Geometry, metric: Metric = "haversine * Karney's algorithm (geographiclib). Unlike Vincenty's iteration, it is * defined for every pair, antipodal points included. */ -export function distanceVincenty(a: Position, b: Position) { +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!; } diff --git a/src/tools/DistanceTool.ts b/src/tools/DistanceTool.ts index f8e9e209..1c5af9cf 100644 --- a/src/tools/DistanceTool.ts +++ b/src/tools/DistanceTool.ts @@ -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,12 @@ 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 et coûteuse, précision à 1mm)", + ". Par défaut : `spherical`." ].join("")), }).strict(); @@ -71,9 +71,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/test/helpers/distance.test.ts b/test/helpers/distance.test.ts index b8ec79cd..98ecce4d 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,13 +41,13 @@ 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); @@ -56,7 +56,7 @@ describe("distance helper", () => { 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(distanceVincenty(paris.coordinates, antipode)).toBeCloseTo(20003931.46, 2); + expect(ellipsoidalDistance(paris.coordinates, antipode)).toBeCloseTo(20003931.46, 2); }); it("computes Paris->Marseille", () => { @@ -760,7 +760,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. @@ -777,7 +777,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", @@ -786,15 +786,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", @@ -810,7 +810,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: { From ffdbb72fbde74276180976e5d6c6faa5be4c8975 Mon Sep 17 00:00:00 2001 From: esgn <5435148+esgn@users.noreply.github.com> Date: Fri, 2 Oct 2026 10:43:12 +0200 Subject: [PATCH 06/11] fix(distance): update description to clarify distance calculation between two geographic positions --- src/tools/DistanceTool.ts | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tools/DistanceTool.ts b/src/tools/DistanceTool.ts index 1c5af9cf..588b55f5 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"; From c99b6334fe2d21b3272901f5e61d6ee3b363f2d4 Mon Sep 17 00:00:00 2001 From: esgn <5435148+esgn@users.noreply.github.com> Date: Fri, 2 Oct 2026 11:03:15 +0200 Subject: [PATCH 07/11] fix(distance): handle a nearest location on the last vertex of a line 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. getSegment then read a vertex past the end and threw a TypeError, which failed urbanisme and assiette_sup calls whose nearest object is a line reached at its end, and silently nulled distance_to_filter_center. Clamp the index to the last segment. --- src/helpers/distance.ts | 5 ++++- test/helpers/distance.test.ts | 13 +++++++++++++ 2 files changed, 17 insertions(+), 1 deletion(-) diff --git a/src/helpers/distance.ts b/src/helpers/distance.ts index 90068e0b..5031108b 100644 --- a/src/helpers/distance.ts +++ b/src/helpers/distance.ts @@ -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]) }; } } diff --git a/test/helpers/distance.test.ts b/test/helpers/distance.test.ts index 98ecce4d..6227045d 100644 --- a/test/helpers/distance.test.ts +++ b/test/helpers/distance.test.ts @@ -527,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", From 0009b0b410c0ecede494f4b4d67af64bf71603ac Mon Sep 17 00:00:00 2001 From: esgn <5435148+esgn@users.noreply.github.com> Date: Fri, 2 Oct 2026 11:05:44 +0200 Subject: [PATCH 08/11] fix(dependencies): remove unused @turf/distance package --- package-lock.json | 1 - package.json | 1 - 2 files changed, 2 deletions(-) diff --git a/package-lock.json b/package-lock.json index cb6bc530..37da64d0 100644 --- a/package-lock.json +++ b/package-lock.json @@ -16,7 +16,6 @@ "@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", diff --git a/package.json b/package.json index dc26f79a..59d39955 100644 --- a/package.json +++ b/package.json @@ -62,7 +62,6 @@ "@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", From 7a1313a5b4c5309f0aead3de440e51d489271c38 Mon Sep 17 00:00:00 2001 From: esgn <5435148+esgn@users.noreply.github.com> Date: Fri, 2 Oct 2026 11:30:42 +0200 Subject: [PATCH 09/11] fix(distance): state the ellipsoidal precision the tool actually returns MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The result is rounded to the centimeter, so 0.5 cm, not 1 mm. Drop "coûteuse", the cost being a few microseconds, and the default already published by the schema. --- docs/mcp-tools.md | 4 ++-- src/tools/DistanceTool.ts | 3 +-- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/docs/mcp-tools.md b/docs/mcp-tools.md index 94ff8dd9..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 : `spherical` distance à vol d'oiseau (Terre ronde, précision à 0.5%), `ellipsoidal` distance à vol d'oiseau (Terre ellipsoïde, plus précise et coûteuse, précision à 1mm). Par défaut : `spherical`. Valeurs : spherical, ellipsoidal. Valeur par défaut : spherical. | +| `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 @@ -2296,7 +2296,7 @@ Renvoie la distance (en mètres) entre deux points à partir de leur longitude e "ellipsoidal" ], "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 et coûteuse, précision à 1mm). Par défaut : `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/src/tools/DistanceTool.ts b/src/tools/DistanceTool.ts index 588b55f5..b4a660a5 100644 --- a/src/tools/DistanceTool.ts +++ b/src/tools/DistanceTool.ts @@ -27,8 +27,7 @@ const distanceInputSchema = z.object({ .default("spherical") .describe(["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 et coûteuse, précision à 1mm)", - ". Par défaut : `spherical`." + " `ellipsoidal` distance à vol d'oiseau (Terre ellipsoïde, plus précise, précision à 0.5cm).", ].join("")), }).strict(); From fd6e2fed41436fb1c50c49be282ed54334b22848 Mon Sep 17 00:00:00 2001 From: esgn <5435148+esgn@users.noreply.github.com> Date: Fri, 2 Oct 2026 13:10:10 +0200 Subject: [PATCH 10/11] test: run the timed tests on their own Wall-clock assertions failed whenever other test files ran in parallel. The four timed tests move to *.perf.test.ts files, which test:perf (vitest.perf.config.mts) runs one file at a time and verify:fast chains after test:unit. The 500-vertex distance test now times the median of warm runs, under 15 ms. --- docs/dev.md | 35 +++++--- package.json | 3 +- test/helpers/distance.perf.test.ts | 57 ++++++++++++ test/helpers/distance.test.ts | 47 ---------- test/wfs/response.perf.test.ts | 135 +++++++++++++++++++++++++++++ test/wfs/response.test.ts | 117 ------------------------- vitest.config.mts | 3 +- vitest.perf.config.mts | 23 +++++ 8 files changed, 241 insertions(+), 179 deletions(-) create mode 100644 test/helpers/distance.perf.test.ts create mode 100644 test/wfs/response.perf.test.ts create mode 100644 vitest.perf.config.mts diff --git a/docs/dev.md b/docs/dev.md index 88aa6502..3c80d322 100644 --- a/docs/dev.md +++ b/docs/dev.md @@ -122,18 +122,19 @@ Les niveaux 1 et 2 nécessitent un build à jour (`npm run build`) et un accès ### Vue d'ensemble des commandes -| Commande | Rôle | -| --------------------------- | --------------------------------------------------------------- | -| `npm run typecheck` | Type-check de l'application (`tsconfig.json`) | -| `npm run typecheck:test` | Type-check des fichiers de test (`tsconfig.test.json`) | -| `npm test` / `test:unit` | Tests unitaires | -| `npm run test:integration` | Tests d'intégration niveau 1 | -| `npm run test:e2e` | Tests E2E agent niveau 2 | -| `npm run test:coverage` | Tests unitaires avec couverture | -| `npm run bench` | Benchmark du calcul de `intersection_area` | -| `npm run verify:fast` | `typecheck` + `typecheck:test` + `build` + `test:unit` | -| `npm run verify` | `verify:fast` + `test:integration` | -| `npm run verify:full` | `verify` + `test:e2e` | +| Commande | Rôle | +| -------------------------- | -------------------------------------------------------------------- | +| `npm run typecheck` | Type-check de l'application (`tsconfig.json`) | +| `npm run typecheck:test` | Type-check des fichiers de test (`tsconfig.test.json`) | +| `npm test` / `test:unit` | Tests unitaires | +| `npm run test:perf` | Tests chronométrés (`*.perf.test.ts`), un fichier à la fois | +| `npm run test:integration` | Tests d'intégration niveau 1 | +| `npm run test:e2e` | Tests E2E agent niveau 2 | +| `npm run test:coverage` | Tests unitaires avec couverture | +| `npm run bench` | Benchmark du calcul de `intersection_area` | +| `npm run verify:fast` | `typecheck` + `typecheck:test` + `build` + `test:unit` + `test:perf` | +| `npm run verify` | `verify:fast` + `test:integration` | +| `npm run verify:full` | `verify` + `test:e2e` | ### Tests unitaires @@ -143,6 +144,14 @@ npm run test:unit npm test ``` +### Tests chronométrés + +Les tests qui mesurent un temps (`*.perf.test.ts`) tournent à part, un fichier à la fois, pour que les autres tests ne faussent pas la mesure. `verify:fast` les lance après les tests unitaires. + +```bash +npm run test:perf +``` + ### Tests d'intégration (niveau 1) ```bash @@ -166,7 +175,7 @@ npm run test:coverage ### Vérifications combinées ```bash -npm run verify:fast # typecheck + build + tests unitaires +npm run verify:fast # typecheck + build + tests unitaires et chronométrés npm run verify # verify:fast + tests d'intégration niveau 1 npm run verify:full # verify + tests E2E niveau 2 ``` diff --git a/package.json b/package.json index 59d39955..90bd66a7 100644 --- a/package.json +++ b/package.json @@ -41,11 +41,12 @@ "start": "node --use-env-proxy dist/index.js", "test": "npm run test:unit", "test:unit": "vitest run", + "test:perf": "vitest run --config vitest.perf.config.mts", "test:integration": "vitest run --config vitest.integration.config.mts", "test:e2e": "vitest run --config vitest.e2e.config.mts", "test:coverage": "vitest run --coverage", "bench": "vitest bench --run --reporter=verbose", - "verify:fast": "npm run typecheck && npm run typecheck:test && npm run build && npm run test:unit", + "verify:fast": "npm run typecheck && npm run typecheck:test && npm run build && npm run test:unit && npm run test:perf", "verify": "npm run verify:fast && npm run test:integration", "verify:full": "npm run verify && npm run test:e2e", "fresh": "npm run reset:local && npm ci && npm run verify", diff --git a/test/helpers/distance.perf.test.ts b/test/helpers/distance.perf.test.ts new file mode 100644 index 00000000..106a4126 --- /dev/null +++ b/test/helpers/distance.perf.test.ts @@ -0,0 +1,57 @@ +import { describe, expect, it } from "vitest"; +import type { Geometry } from "geojson"; + +import distance from "../../src/helpers/distance.js"; + +describe("distance helper", () => { + describe("performance sanity", () => { + it("resolves two disjoint 500-vertex polygons in under 15 milliseconds once warm", () => { + const circle = (centerLon: number, centerLat: number, radiusDeg: number, n: number): Geometry => ({ + type: "Polygon", + coordinates: [ + Array.from({ length: n + 1 }, (_, i) => { + const a = (2 * Math.PI * i) / n; + return [centerLon + radiusDeg * Math.cos(a), centerLat + radiusDeg * Math.sin(a)]; + }), + ], + }); + + const a = circle(2.35, 48.85, 0.3, 500); + const b = circle(2.35, 53.85, 0.3, 500); + + expect(distance(a, b).distance).toBeGreaterThan(0); + + // Time warm runs only: the first calls run before V8 has optimized the code. + for (let i = 0; i < 5; i++) distance(a, b); + const samples = Array.from({ length: 15 }, () => { + const start = performance.now(); + distance(a, b); + return performance.now() - start; + }).sort((x, y) => x - y); + expect(samples[7]).toBeLessThan(15); // median of the warm runs + }); + + it("resolves two nearby irregular ~10000-vertex polygons in under a second", () => { + const blob = (centerLon: number, centerLat: number, radiusDeg: number, n: number, seed: number): Geometry => ({ + type: "Polygon", + coordinates: [ + Array.from({ length: n + 1 }, (_, i) => { + const a = (2 * Math.PI * i) / n; + const wobble = 1 + 0.15 * Math.sin(a * 7 + seed) + 0.08 * Math.sin(a * 13 + seed * 2); + return [centerLon + radiusDeg * wobble * Math.cos(a), centerLat + radiusDeg * wobble * Math.sin(a)]; + }), + ], + }); + + const a = blob(2.0, 48.0, 0.5, 10000, 1); + const b = blob(4.0, 48.3, 0.5, 10000, 2); + + const start = performance.now(); + const result = distance(a, b); + const elapsedMs = performance.now() - start; + + expect(result.distance).toBeGreaterThan(0); + expect(elapsedMs).toBeLessThan(1_000); + }); + }); +}); diff --git a/test/helpers/distance.test.ts b/test/helpers/distance.test.ts index 6227045d..0f793a26 100644 --- a/test/helpers/distance.test.ts +++ b/test/helpers/distance.test.ts @@ -828,51 +828,4 @@ describe("distance helper", () => { expectCloseRatio(d, meridianEdge, 0.001, `haversine edge ranking at lat ${lat}`); }); }); - - describe("performance sanity", () => { - it("resolves two disjoint 500-vertex polygons in under 10 millisecond", () => { - const circle = (centerLon: number, centerLat: number, radiusDeg: number, n: number): Geometry => ({ - type: "Polygon", - coordinates: [ - Array.from({ length: n + 1 }, (_, i) => { - const a = (2 * Math.PI * i) / n; - return [centerLon + radiusDeg * Math.cos(a), centerLat + radiusDeg * Math.sin(a)]; - }), - ], - }); - - const a = circle(2.35, 48.85, 0.3, 500); - const b = circle(2.35, 53.85, 0.3, 500); - - const start = performance.now(); - const result = distance(a, b); - const elapsedMs = performance.now() - start; - - expect(result.distance).toBeGreaterThan(0); - expect(elapsedMs).toBeLessThan(10); - }); - - it("resolves two nearby irregular ~10000-vertex polygons in under a second", () => { - const blob = (centerLon: number, centerLat: number, radiusDeg: number, n: number, seed: number): Geometry => ({ - type: "Polygon", - coordinates: [ - Array.from({ length: n + 1 }, (_, i) => { - const a = (2 * Math.PI * i) / n; - const wobble = 1 + 0.15 * Math.sin(a * 7 + seed) + 0.08 * Math.sin(a * 13 + seed * 2); - return [centerLon + radiusDeg * wobble * Math.cos(a), centerLat + radiusDeg * wobble * Math.sin(a)]; - }), - ], - }); - - const a = blob(2.0, 48.0, 0.5, 10000, 1); - const b = blob(4.0, 48.3, 0.5, 10000, 2); - - const start = performance.now(); - const result = distance(a, b); - const elapsedMs = performance.now() - start; - - expect(result.distance).toBeGreaterThan(0); - expect(elapsedMs).toBeLessThan(1_000); - }); - }); }); diff --git a/test/wfs/response.perf.test.ts b/test/wfs/response.perf.test.ts new file mode 100644 index 00000000..b312b0f6 --- /dev/null +++ b/test/wfs/response.perf.test.ts @@ -0,0 +1,135 @@ +import { performance } from "node:perf_hooks"; +import { describe, expect, it } from "vitest"; +import { intersect } from "@turf/intersect"; +import { feature as turfFeature, featureCollection } from "@turf/helpers"; +import type { Polygon } from "geojson"; +import area from "../../src/helpers/area"; + +import { transformFeatureCollectionResponse } from "../../src/wfs/response"; +import { type SpatialExtraOptions } from "../../src/wfs/schema"; + +describe("wfs_engine/response", () => { + function getFeatures( + result: ReturnType, + ): NonNullable["features"]> { + expect(result.features).toBeDefined(); + return result.features!; + } + + describe("transformFeatureCollectionResponse", () => { + it("should stay fast when computing centroid, area, and intersection_area for a large region against many polygons crossing its boundary", () => { + const regionCenterLon = 2.35; + const regionCenterLat = 48.85; + const regionVertices = 3000; + const polygonSides = 20; + const featureCount = 5000; + + /** Close a ring by repeating its first position, as GeoJSON requires. */ + function closeRing(positions: T[]) { + return [...positions, positions[0]]; + } + + function regularRing(centerLon: number, centerLat: number, radius: number, sides: number) { + return closeRing(Array.from({ length: sides }, (_, index) => { + const angle = (Math.PI * 2 * index) / sides - Math.PI / 2; + return [centerLon + radius * Math.cos(angle), centerLat + radius * Math.sin(angle)] as [number, number]; + })); + } + + /** Position on the region boundary at the given angle from its center. */ + function regionBoundary(angle: number) { + const radius = 0.22 + 0.03 * Math.sin(7 * angle); + return [ + regionCenterLon + radius * Math.cos(angle), + regionCenterLat + radius * Math.sin(angle), + ]; + } + + const regionGeometry = { + type: "Polygon" as const, + coordinates: [ + closeRing(Array.from({ length: regionVertices }, (_, index) => regionBoundary((Math.PI * 2 * index) / regionVertices))), + ], + }; + + // Every polygon is centered on the region boundary, so that each one needs + // an actual intersection to be computed. + const featureCollection = { + type: "FeatureCollection", + features: Array.from({ length: featureCount }, (_, index) => { + const [centerLon, centerLat] = regionBoundary((Math.PI * 2 * index) / featureCount); + const radius = 0.002 + (index % 7) * 0.00025; + + return { + id: `poly.${index}`, + geometry: { + type: "Polygon", + coordinates: [regularRing(centerLon, centerLat, radius, polygonSides)], + }, + properties: { name: `poly-${index}` }, + }; + }), + }; + + const input = { + typename: "TEST:type", + spatial_extras: ["centroid", "area", "intersection_area"] as SpatialExtraOptions[], + intersects_feature_filter: { + typename: "TEST:region", + feature_id: "region.1", + }, + }; + + const start = performance.now(); + const result = transformFeatureCollectionResponse(featureCollection, input, regionGeometry); + const elapsedMs = performance.now() - start; + + const features = getFeatures(result); + expect(features).toHaveLength(featureCount); + expect(features[0].centroid).toMatchObject({ lon: expect.any(Number), lat: expect.any(Number) }); + const notPartiallyIntersecting = features.filter((feature) => !( + typeof feature.intersection_area === "number" + && feature.intersection_area > 0 + && feature.intersection_area < (feature.area as number) + )); + expect(notPartiallyIntersecting).toEqual([]); + expect(elapsedMs).toBeLessThan(1000); + }); + + it("should stay fast when computing intersection_area for a large geometry with holes crossing a detailed reference boundary", () => { + // Like a forest with clearings across the boundary of an isochrone: both + // geometries have many positions in the same place, which used to cost + // their product (over a second here). + function wavyRing(centerLon: number, centerLat: number, radius: number, vertices: number, waves: number) { + const positions = Array.from({ length: vertices }, (_, index) => { + const angle = (Math.PI * 2 * index) / vertices; + const wavyRadius = radius * (1 + 0.02 * Math.sin(waves * angle)); + return [centerLon + wavyRadius * Math.cos(angle), centerLat + wavyRadius * Math.sin(angle)]; + }); + return [...positions, positions[0]]; + } + + const reference: Polygon = { type: "Polygon", coordinates: [wavyRing(2.35, 48.85, 0.3, 20000, 200)] }; + // Centered on the eastern boundary of the reference, with clearings on a grid. + const clearings = Array.from({ length: 100 }, (_, index) => [2.65 + 0.014 * (index % 10 - 4.5), 48.85 + 0.014 * (Math.floor(index / 10) - 4.5)]) + .filter(([lon, lat]) => Math.hypot(lon - 2.65, lat - 48.85) < 0.08) + .map(([lon, lat]) => wavyRing(lon, lat, 0.003, 16, 3).reverse()); + const forest: Polygon = { type: "Polygon", coordinates: [wavyRing(2.65, 48.85, 0.1, 20000, 300), ...clearings] }; + + const start = performance.now(); + const result = transformFeatureCollectionResponse({ + type: "FeatureCollection", + features: [{ id: "forest.1", geometry: forest, properties: {} }], + }, { + typename: "TEST:type", + spatial_extras: ["intersection_area"], + intersects_feature_filter: { typename: "TEST:region", feature_id: "region.1" }, + }, reference); + const elapsedMs = performance.now() - start; + + const inter = intersect(featureCollection([turfFeature(forest), turfFeature(reference)])); + expect(getFeatures(result)[0].intersection_area as number).toBeCloseTo(area(inter!.geometry), 3); + expect(elapsedMs).toBeLessThan(300); + }); + }); +}); diff --git a/test/wfs/response.test.ts b/test/wfs/response.test.ts index 2aac45fe..46c0a207 100644 --- a/test/wfs/response.test.ts +++ b/test/wfs/response.test.ts @@ -1,4 +1,3 @@ -import { performance } from "node:perf_hooks"; import { describe, expect, it } from "vitest"; import { intersect } from "@turf/intersect"; import { feature as turfFeature, featureCollection } from "@turf/helpers"; @@ -11,7 +10,6 @@ import { transformFeatureCollectionResponse, postProcessFeatureCollection, } from "../../src/wfs/response"; -import { type SpatialExtraOptions } from "../../src/wfs/schema"; describe("wfs_engine/response", () => { function getFeatures( @@ -212,121 +210,6 @@ describe("wfs_engine/response", () => { expect(features[0].distance_to_filter_center as number).toBeCloseTo(2340.99, 3); }); - it("should stay fast when computing centroid, area, and intersection_area for a large region against many polygons crossing its boundary", () => { - const regionCenterLon = 2.35; - const regionCenterLat = 48.85; - const regionVertices = 3000; - const polygonSides = 20; - const featureCount = 5000; - - /** Close a ring by repeating its first position, as GeoJSON requires. */ - function closeRing(positions: T[]) { - return [...positions, positions[0]]; - } - - function regularRing(centerLon: number, centerLat: number, radius: number, sides: number) { - return closeRing(Array.from({ length: sides }, (_, index) => { - const angle = (Math.PI * 2 * index) / sides - Math.PI / 2; - return [centerLon + radius * Math.cos(angle), centerLat + radius * Math.sin(angle)] as [number, number]; - })); - } - - /** Position on the region boundary at the given angle from its center. */ - function regionBoundary(angle: number) { - const radius = 0.22 + 0.03 * Math.sin(7 * angle); - return [ - regionCenterLon + radius * Math.cos(angle), - regionCenterLat + radius * Math.sin(angle), - ]; - } - - const regionGeometry = { - type: "Polygon" as const, - coordinates: [ - closeRing(Array.from({ length: regionVertices }, (_, index) => regionBoundary((Math.PI * 2 * index) / regionVertices))), - ], - }; - - // Every polygon is centered on the region boundary, so that each one needs - // an actual intersection to be computed. - const featureCollection = { - type: "FeatureCollection", - features: Array.from({ length: featureCount }, (_, index) => { - const [centerLon, centerLat] = regionBoundary((Math.PI * 2 * index) / featureCount); - const radius = 0.002 + (index % 7) * 0.00025; - - return { - id: `poly.${index}`, - geometry: { - type: "Polygon", - coordinates: [regularRing(centerLon, centerLat, radius, polygonSides)], - }, - properties: { name: `poly-${index}` }, - }; - }), - }; - - const input = { - typename: "TEST:type", - spatial_extras: ["centroid", "area", "intersection_area"] as SpatialExtraOptions[], - intersects_feature_filter: { - typename: "TEST:region", - feature_id: "region.1", - }, - }; - - const start = performance.now(); - const result = transformFeatureCollectionResponse(featureCollection, input, regionGeometry); - const elapsedMs = performance.now() - start; - - const features = getFeatures(result); - expect(features).toHaveLength(featureCount); - expect(features[0].centroid).toMatchObject({ lon: expect.any(Number), lat: expect.any(Number) }); - const notPartiallyIntersecting = features.filter((feature) => !( - typeof feature.intersection_area === "number" - && feature.intersection_area > 0 - && feature.intersection_area < (feature.area as number) - )); - expect(notPartiallyIntersecting).toEqual([]); - expect(elapsedMs).toBeLessThan(1000); - }); - - it("should stay fast when computing intersection_area for a large geometry with holes crossing a detailed reference boundary", () => { - // Like a forest with clearings across the boundary of an isochrone: both - // geometries have many positions in the same place, which used to cost - // their product (over a second here). - function wavyRing(centerLon: number, centerLat: number, radius: number, vertices: number, waves: number) { - const positions = Array.from({ length: vertices }, (_, index) => { - const angle = (Math.PI * 2 * index) / vertices; - const wavyRadius = radius * (1 + 0.02 * Math.sin(waves * angle)); - return [centerLon + wavyRadius * Math.cos(angle), centerLat + wavyRadius * Math.sin(angle)]; - }); - return [...positions, positions[0]]; - } - - const reference: Polygon = { type: "Polygon", coordinates: [wavyRing(2.35, 48.85, 0.3, 20000, 200)] }; - // Centered on the eastern boundary of the reference, with clearings on a grid. - const clearings = Array.from({ length: 100 }, (_, index) => [2.65 + 0.014 * (index % 10 - 4.5), 48.85 + 0.014 * (Math.floor(index / 10) - 4.5)]) - .filter(([lon, lat]) => Math.hypot(lon - 2.65, lat - 48.85) < 0.08) - .map(([lon, lat]) => wavyRing(lon, lat, 0.003, 16, 3).reverse()); - const forest: Polygon = { type: "Polygon", coordinates: [wavyRing(2.65, 48.85, 0.1, 20000, 300), ...clearings] }; - - const start = performance.now(); - const result = transformFeatureCollectionResponse({ - type: "FeatureCollection", - features: [{ id: "forest.1", geometry: forest, properties: {} }], - }, { - typename: "TEST:type", - spatial_extras: ["intersection_area"], - intersects_feature_filter: { typename: "TEST:region", feature_id: "region.1" }, - }, reference); - const elapsedMs = performance.now() - start; - - const inter = intersect(featureCollection([turfFeature(forest), turfFeature(reference)])); - expect(getFeatures(result)[0].intersection_area as number).toBeCloseTo(area(inter!.geometry), 3); - expect(elapsedMs).toBeLessThan(300); - }); - // Regression: `dwithin_point` used to short-circuit `intersection_area` and // return the feature's own geometry, on the false premise that a DWITHIN // match implies containment. DWITHIN matches as soon as ANY part of the diff --git a/vitest.config.mts b/vitest.config.mts index 7e36da44..d8c367f2 100644 --- a/vitest.config.mts +++ b/vitest.config.mts @@ -15,7 +15,8 @@ export default defineConfig({ globals: false, environment: "node", include: ["test/**/*.test.ts"], - exclude: ["test/integration/**/*"], + // Timed tests run on their own (`test:perf`, vitest.perf.config.mts). + exclude: ["test/integration/**/*", "test/**/*.perf.test.ts"], testTimeout: 60 * MILLISECONDS, // Default value, set explicitly because the suite relies on per-file isolation // (it also stops Vitest from suggesting `isolate: false`). diff --git a/vitest.perf.config.mts b/vitest.perf.config.mts new file mode 100644 index 00000000..e98a7e5c --- /dev/null +++ b/vitest.perf.config.mts @@ -0,0 +1,23 @@ +import { defineConfig } from "vitest/config"; + +const MILLISECONDS = 1000; + +export default defineConfig({ + resolve: { + alias: [ + { + find: /^(\.{1,2}\/.*)\.js$/, + replacement: "$1", + }, + ], + }, + test: { + globals: false, + environment: "node", + include: ["test/**/*.perf.test.ts"], + testTimeout: 60 * MILLISECONDS, + // Run one file at a time: these tests measure wall-clock time, which test + // files running in parallel would inflate. + fileParallelism: false, + }, +}); From e252e3ed73d5263479e2599376c5b84026d24812 Mon Sep 17 00:00:00 2001 From: Lionel Zoubritzky Date: Fri, 2 Oct 2026 13:57:09 +0200 Subject: [PATCH 11/11] chore(dependencies): remove unused `@types/rbush` package --- package-lock.json | 8 -------- package.json | 1 - 2 files changed, 9 deletions(-) diff --git a/package-lock.json b/package-lock.json index 37da64d0..aeb74494 100644 --- a/package-lock.json +++ b/package-lock.json @@ -42,7 +42,6 @@ "@turf/intersect": "^7.4.0", "@types/geojson": "^7946.0.16", "@types/node": "^26.6.3", - "@types/rbush": "^4.0.0", "@types/supertest": "^7.2.1", "@vitest/coverage-v8": "^5.0.3", "ajv": "^8.20.0", @@ -1681,13 +1680,6 @@ "kleur": "^3.0.3" } }, - "node_modules/@types/rbush": { - "version": "4.0.0", - "resolved": "https://registry.npmjs.org/@types/rbush/-/rbush-4.0.0.tgz", - "integrity": "sha512-+N+2H39P8X+Hy1I5mC6awlTX54k3FhiUmvt7HWzGJZvF+syUAAxP/stwppS8JE84YHqFgRMv6fCy31202CMFxQ==", - "dev": true, - "license": "MIT" - }, "node_modules/@types/superagent": { "version": "8.1.11", "resolved": "https://registry.npmjs.org/@types/superagent/-/superagent-8.1.11.tgz", diff --git a/package.json b/package.json index 90bd66a7..02bc3865 100644 --- a/package.json +++ b/package.json @@ -85,7 +85,6 @@ "@turf/intersect": "^7.4.0", "@types/geojson": "^7946.0.16", "@types/node": "^26.6.3", - "@types/rbush": "^4.0.0", "@types/supertest": "^7.2.1", "@vitest/coverage-v8": "^5.0.3", "ajv": "^8.20.0",