refactoring

website/src/lib/geo/closest-point.test.ts

@@ -0,0 +1,50 @@
+import { describe, it } from 'node:test';
+import assert from 'node:assert/strict';
+import { closestPointIndex, closestPointIndexIn } from '$lib/geo/closest-point';
+
+describe('closestPointIndex', () => {
+    const line = [0, 0, 1, 0, 2, 0, 2, 1];
+
+    it('has no point on an empty polyline', () => {
+        assert.equal(closestPointIndex([], { lng: 0, lat: 0 }), undefined);
+    });
+
+    it('is the only point of a polyline of one point', () => {
+        assert.equal(closestPointIndex([3, 4], { lng: 0, lat: 0 }), 0);
+    });
+
+    it('is the nearest end of the nearest segment', () => {
+        assert.equal(closestPointIndex(line, { lng: 0.2, lat: 0.1 }), 0);
+        assert.equal(closestPointIndex(line, { lng: 0.8, lat: 0.1 }), 1);
+        assert.equal(closestPointIndex(line, { lng: 1.6, lat: -0.1 }), 2);
+        assert.equal(closestPointIndex(line, { lng: 2.1, lat: 0.9 }), 3);
+    });
+
+    it('uses the distance to the segment, not to its ends', () => {
+        // far from every point but right next to the long segment
+        const coordinates = [0, 0, 10, 0, 10, 5];
+        assert.equal(closestPointIndex(coordinates, { lng: 7, lat: 0.01 }), 1);
+        assert.equal(closestPointIndex(coordinates, { lng: 3, lat: 0.01 }), 0);
+    });
+});
+
+describe('closestPointIndexIn', () => {
+    const lng = [0, 1, 2, 2, 5, 6];
+    const lat = [0, 0, 0, 1, 5, 5];
+
+    it('gives the index in the arrays, looking at the given range only', () => {
+        assert.equal(closestPointIndexIn(lng, lat, 0, 4, { lng: 0.8, lat: 0.1 }), 1);
+        // the range starts after the closest points
+        assert.equal(closestPointIndexIn(lng, lat, 4, 6, { lng: 0.8, lat: 0.1 }), 4);
+        assert.equal(closestPointIndexIn(lng, lat, 4, 6, { lng: 5.9, lat: 5.1 }), 5);
+    });
+
+    it('has no point in an empty range', () => {
+        assert.equal(closestPointIndexIn(lng, lat, 3, 3, { lng: 0, lat: 0 }), undefined);
+        assert.equal(closestPointIndexIn(lng, lat, 4, 2, { lng: 0, lat: 0 }), undefined);
+    });
+
+    it('is the point itself for a range of one point', () => {
+        assert.equal(closestPointIndexIn(lng, lat, 3, 4, { lng: 9, lat: 9 }), 3);
+    });
+});

website/src/lib/geo/closest-point.ts

@@ -0,0 +1,71 @@
+const metersPerDegree = 111320;
+
+/**
+ * Index of the point of a polyline (flat `[lng, lat, ...]` coordinates) that is the closest to
+ * `point`, among the two ends of the segment of the polyline that passes the closest to it.
+ * `undefined` for an empty polyline. Planar approximation, like the one used to simplify tracks.
+ */
+export function closestPointIndex(
+    coordinates: ArrayLike<number>,
+    point: { lng: number; lat: number }
+): number | undefined {
+    return closestIndex(
+        Math.floor(coordinates.length / 2),
+        (i) => coordinates[2 * i],
+        (i) => coordinates[2 * i + 1],
+        point
+    );
+}
+
+/**
+ * Like `closestPointIndex`, for the points `from` (included) to `to` (excluded) of a polyline
+ * given as the arrays of its longitudes and latitudes. The index is the one in the arrays.
+ */
+export function closestPointIndexIn(
+    lng: ArrayLike<number>,
+    lat: ArrayLike<number>,
+    from: number,
+    to: number,
+    point: { lng: number; lat: number }
+): number | undefined {
+    const closest = closestIndex(
+        Math.max(0, to - from),
+        (i) => lng[from + i],
+        (i) => lat[from + i],
+        point
+    );
+    return closest === undefined ? undefined : from + closest;
+}
+
+function closestIndex(
+    count: number,
+    lngAt: (index: number) => number,
+    latAt: (index: number) => number,
+    point: { lng: number; lat: number }
+): number | undefined {
+    if (count === 0) {
+        return undefined;
+    }
+    const kx = Math.cos((point.lat * Math.PI) / 180) * metersPerDegree;
+    const ky = metersPerDegree;
+    const px = point.lng * kx;
+    const py = point.lat * ky;
+    const x = (i: number) => lngAt(i) * kx - px;
+    const y = (i: number) => latAt(i) * ky - py;
+
+    let best = 0;
+    let bestDistance = Number.MAX_VALUE;
+    for (let i = 0; i < count - 1; i++) {
+        const [ax, ay, bx, by] = [x(i), y(i), x(i + 1), y(i + 1)];
+        const [dx, dy] = [bx - ax, by - ay];
+        const length = dx * dx + dy * dy;
+        // the point is at the origin: the closest position on the segment, as a ratio
+        const t = length === 0 ? 0 : Math.max(0, Math.min(1, -(ax * dx + ay * dy) / length));
+        const distance = Math.hypot(ax + t * dx, ay + t * dy);
+        if (distance < bestDistance) {
+            bestDistance = distance;
+            best = t < 0.5 ? i : i + 1;
+        }
+    }
+    return best;
+}

website/src/lib/geo/coordinates.test.ts

@@ -0,0 +1,31 @@
+import assert from 'node:assert/strict';
+import { describe, it } from 'node:test';
+import { distance } from '$lib/geo/coordinates';
+
+describe('distance', () => {
+    it('is zero between a point and itself', () => {
+        assert.equal(distance({ lng: 4, lat: 50 }, { lng: 4, lat: 50 }), 0);
+    });
+
+    it('measures a degree of latitude and longitude', () => {
+        const meridian = distance({ lng: 0, lat: 0 }, { lng: 0, lat: 1 });
+        assert.ok(Math.abs(meridian - 111195) < 10, `${meridian}`);
+        const equator = distance({ lng: 0, lat: 0 }, { lng: 1, lat: 0 });
+        assert.ok(Math.abs(equator - 111195) < 10, `${equator}`);
+        // degrees of longitude get shorter away from the equator
+        const north = distance({ lng: 0, lat: 60 }, { lng: 1, lat: 60 });
+        assert.ok(Math.abs(north - 111195 / 2) < 100, `${north}`);
+    });
+
+    it('is symmetric and handles the antipodes and the antimeridian', () => {
+        const a = { lng: 4.35, lat: 50.85 };
+        const b = { lng: -74, lat: 40.7 };
+        assert.equal(distance(a, b), distance(b, a));
+        assert.ok(Math.abs(distance(a, b) - 5900_000) < 100_000);
+        const half = Math.PI * 6371008.8;
+        assert.ok(Math.abs(distance({ lng: 0, lat: 0 }, { lng: 180, lat: 0 }) - half) < 1);
+        assert.ok(Math.abs(distance({ lng: 0, lat: 90 }, { lng: 0, lat: -90 }) - half) < 1);
+        const across = distance({ lng: 179.5, lat: 0 }, { lng: -179.5, lat: 0 });
+        assert.ok(Math.abs(across - 111195) < 10, `${across}`);
+    });
+});

website/src/lib/geo/coordinates.ts

@@ -0,0 +1,11 @@
+import maplibregl from 'maplibre-gl';
+
+export type Coordinates = {
+    lat: number;
+    lng: number;
+};
+
+/** The distance in meters between two coordinates (haversine formula). */
+export function distance(a: Coordinates, b: Coordinates): number {
+    return new maplibregl.LngLat(a.lng, a.lat).distanceTo(new maplibregl.LngLat(b.lng, b.lat));
+}