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/docs/mcp-tools.md b/docs/mcp-tools.md
index f8c9cf24..d1c2045c 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 : `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
+
+```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": [
+ "spherical",
+ "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, précision à 0.5cm)."
+ }
+ },
+ "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/package-lock.json b/package-lock.json
index 3b03e050..aeb74494 100644
--- a/package-lock.json
+++ b/package-lock.json
@@ -16,10 +16,11 @@
"@turf/bbox-polygon": "^7.4.0",
"@turf/centroid": "^7.4.0",
"@turf/circle": "^7.4.0",
- "@turf/distance": "^7.4.0",
"@turf/helpers": "^7.4.0",
"@turf/length": "^7.4.0",
+ "@turf/midpoint": "^7.4.0",
"earcut": "^3.2.4",
+ "geographiclib-geodesic": "^2.2.0",
"jsts": "^2.12.1",
"mcp-framework": "^0.2.22",
"node-fetch": "^3.3.2",
@@ -992,22 +993,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": {
@@ -1412,6 +1416,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",
@@ -1546,10 +1565,27 @@
"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.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,
@@ -3369,6 +3405,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",
diff --git a/package.json b/package.json
index df08b4d9..02bc3865 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",
@@ -62,10 +63,11 @@
"@turf/bbox-polygon": "^7.4.0",
"@turf/centroid": "^7.4.0",
"@turf/circle": "^7.4.0",
- "@turf/distance": "^7.4.0",
"@turf/helpers": "^7.4.0",
"@turf/length": "^7.4.0",
+ "@turf/midpoint": "^7.4.0",
"earcut": "^3.2.4",
+ "geographiclib-geodesic": "^2.2.0",
"jsts": "^2.12.1",
"mcp-framework": "^0.2.22",
"node-fetch": "^3.3.2",
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).
+ *
+ * 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. 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 {
+ distance: number;
+ point1: Position;
+ point2: Position;
+}
+
+/** Point-to-point metric: spherical great-circle, or geodesic on WGS84. */
+export type Metric = "spherical" | "ellipsoidal";
+
+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;
+
+/** 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);
+
+/**
+ * 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;
+ 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));
+}
/**
- * Compute approximative distance in meters between gA and gB.
- *
- * TODO: replace the lon/lat planar nearest-point step with a geodesic geometry distance.
+ * Throws if the edge has a longitude jump > 180°, indicating an
+ * antimeridian-crossing geometry that SHOULD be split per RFC 7946.
+ */
+function checkEdge(c1: Position, c2: Position) {
+ if (Math.abs(c1[0] - c2[0]) > 180) {
+ throw new Error(
+ "Antimeridian-crossing geometries SHOULD be split per RFC 7946 " +
+ "(Polygon → MultiPolygon, LineString → MultiLineString)",
+ );
+ }
+}
+
+/** 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 "LineString":
+ 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) {
+ for (let i = 0; i < line.length - 1; i++) checkEdge(line[i], line[i + 1]);
+ }
+ break;
+ case "Polygon":
+ for (const ring of g.coordinates) {
+ 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) {
+ for (let i = 0; i < ring.length - 1; i++) checkEdge(ring[i], ring[i + 1]);
+ }
+ break;
+ case "GeometryCollection":
+ for (const child of g.geometries) process(child);
+ break;
+ }
+ }
+ process(geom);
+}
+
+/** 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;
+
+ 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 > TOLERANCE_M; // stop when step < 5 mm, since distance result is rounded to the cm
+ 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) };
+}
+
+/**
+ * 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 === "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 {
+ meridional: WGS84_SEMI_MAJOR_M * (1 - WGS84_ECCENTRICITY_SQ) / (w * Math.sqrt(w)),
+ primeVertical: WGS84_SEMI_MAJOR_M / Math.sqrt(w),
+ };
+}
+
+/**
+ * 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 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 "spherical", 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 {
+ 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)) };
+ }
+}
+
+/**
+ * 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[];
+}
+
+/**
+ * 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) {
+ // 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]) };
+ }
+ }
+
+ // 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.
*
* @param {object} gA GeoJSON Geometry
* @param {object} gB GeoJSON Geometry
+ * @param metric Point-to-point metric: "spherical" (default, haversine) or "ellipsoidal" (geographiclib).
*/
-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: 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 };
+
+ 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] };
+ }
+
+ // 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;
+ }
+ if (!converged) {
+ throw new Error("Convergence error in the distance algorithm: cannot compute the distance.");
+ }
- /*
- * 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' }
- );
+ return { ...best, distance: Math.round(best.distance * 100) / 100 };
}
+
+/**
+ * Computes the geodesic distance on WGS84 between two lon/lat points, with
+ * Karney's algorithm (geographiclib). Unlike Vincenty's iteration, it is
+ * defined for every pair, antipodal points included.
+ */
+export function ellipsoidalDistance(a: Position, b: Position) {
+ // Inverse takes lat1, lon1, lat2, lon2; asking for DISTANCE always sets s12.
+ return geodesic.Geodesic.WGS84.Inverse(a[1], a[0], b[1], b[0], geodesic.Geodesic.DISTANCE).s12!;
+}
+
+export default distance;
diff --git a/src/tools/DistanceTool.ts b/src/tools/DistanceTool.ts
new file mode 100644
index 00000000..b4a660a5
--- /dev/null
+++ b/src/tools/DistanceTool.ts
@@ -0,0 +1,92 @@
+/**
+ * MCP tool exposing the distance between two geographic positions.
+ */
+
+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 { ellipsoidalDistance, 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(["spherical", "ellipsoidal"])
+ .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, précision à 0.5cm).",
+ ].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 "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]
+ );
+ 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/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.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 09f9586f..0f793a26 100644
--- a/test/helpers/distance.test.ts
+++ b/test/helpers/distance.test.ts
@@ -1,39 +1,831 @@
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, { ellipsoidalDistance, 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?: "spherical" | "ellipsoidal",
+) {
+ 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 an ellipsoidal point metric", () => {
+ const defaultDistance = distance(paris, marseille).distance;
+ 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(ellipsoidal);
+ expect(result.distance).not.toBe(defaultDistance);
+ expect(result.point1).toEqual(paris.coordinates);
+ expect(result.point2).toEqual(marseille.coordinates);
+ });
+
+ it("computes the ellipsoidal distance between antipodal points", () => {
+ // Vincenty's iteration has no answer here; Karney's algorithm goes over the pole.
+ const antipode = [paris.coordinates[0] - 180, -paris.coordinates[1]];
+ expect(ellipsoidalDistance(paris.coordinates, antipode)).toBeCloseTo(20003931.46, 2);
+ });
+
+ it("computes Paris->Marseille", () => {
+ const result = distance(paris, marseille);
+ expectCloseRatio(result.distance, 662_488.38, 0.0005, "Paris-Marseille distance");
+ 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");
+ });
+
+ // 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",
+ 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).
+ //
+ // 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([
+ {
+ 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,
+ },
+ {
+ 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: [polarRing(80)] } as Geometry,
+ },
+ {
+ 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: [
+ [[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: "MultiPolygon",
+ coordinates: [
+ [[[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]]],
+ [[[170, 0], [170, 10], [-170, 10], [-170, 0], [170, 0]]],
+ ],
+ } as Geometry,
+ },
+ {
+ 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)", () => {
+ 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("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",
+ 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/));
+
+ // 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: [
+ [-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, "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, "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");
+ });
+
+ // 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: "LineString",
+ coordinates: [
+ [-33.47427845001221, 15.573780059814453],
+ [-38.53770555856234, -89],
+ ],
+ };
+ expectThrowBothWays(a, b, /Convergence error/);
+ });
+
+ 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);
+ }
+ });
+ });
+
+ // 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 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.
+ 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 on the ellipsoid 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 = 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, `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 on the sphere 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, `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/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..1f175d18
--- /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: ["spherical", "ellipsoidal"],
+ default: "spherical",
+ });
+ expect(tool.toolDefinition.outputSchema).toBeDefined();
+ });
+
+ 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: {
+ 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(),
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 a503b4a1..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(
@@ -180,7 +178,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,122 +207,7 @@ describe("wfs_engine/response", () => {
});
const features = getFeatures(result);
- expect(features[0].distance_to_filter_center as number).toBeCloseTo(2340.9971606708805, 6);
- });
-
- 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);
+ expect(features[0].distance_to_filter_center as number).toBeCloseTo(2340.99, 3);
});
// Regression: `dwithin_point` used to short-circuit `intersection_area` and
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,
+ },
+});