https://bikeatthebeach.com/js/geo/track-math.js

1// js/geo/track-math.js
2//
3// Pure measurement helpers for GPS tracks. No map, no DOM, no state.
4//
5// Written for the trip-import flow rather than reached for from Turf on
6// purpose. That flow is used on a phone in the field, and Turf arrives
7// through the page's import map from a CDN; every dependency in the
8// upload path is one more thing that can fail on a weak connection at a
9// put-in. Nothing here reaches for Turf, and every function is small
10// enough to read in one sitting.
11//
12// The trip-import flow uses the first five. The rest were added on
13// August 22, 2026 for the route corridor and day-stage work, which needs
14// the same measurements against a route rather than a recorded track.
15//
16// Coordinates are GeoJSON order throughout: [lng, lat] or [lng, lat, ele].
17
18const EARTH_RADIUS_M = 6371008.8;   // IUGG mean radius
19const M_PER_DEG_LAT  = 111320;      // close enough anywhere on the globe
20const toRad = deg => (deg * Math.PI) / 180;
21
22/**
23 * Great-circle distance in meters between two [lng, lat] positions.
24 * Haversine is accurate to well under a meter at trip scale, which is
25 * far finer than consumer GPS.
26 */
27export function haversineMeters(a, b) {
28  const [lng1, lat1] = a;
29  const [lng2, lat2] = b;
30  const dLat = toRad(lat2 - lat1);
31  const dLng = toRad(lng2 - lng1);
32  const s =
33    Math.sin(dLat / 2) ** 2 +
34    Math.cos(toRad(lat1)) * Math.cos(toRad(lat2)) * Math.sin(dLng / 2) ** 2;
35  return 2 * EARTH_RADIUS_M * Math.asin(Math.min(1, Math.sqrt(s)));
36}
37
38/** Total length of a coordinate ring in meters. */
39export function lineLengthMeters(coords) {
40  if (!Array.isArray(coords) || coords.length < 2) return 0;
41  let total = 0;
42  for (let i = 1; i < coords.length; i++) {
43    total += haversineMeters(coords[i - 1], coords[i]);
44  }
45  return total;
46}
47
48/**
49 * Remove `meters` of track from the start and the same from the end.
50 *
51 * This is the privacy trim. Someone who paddles from their own dock or
52 * rides from their driveway has their home address sitting in the first
53 * few hundred meters of the file, and blurring the whole track would
54 * ruin the observation it exists to support. Cutting the ends keeps the
55 * middle exact.
56 *
57 * The cut lands on an interpolated position rather than the nearest
58 * recorded point, so "500 m" means 500 m and not "somewhere between 480
59 * and 690 depending on how often the watch sampled".
60 *
61 * Returns null when the track is too short to survive the trim — the
62 * caller must treat that as "do not share this track" rather than
63 * falling back to the untrimmed original.
64 */
65export function trimLineEnds(coords, meters) {
66  if (!Array.isArray(coords) || coords.length < 2) return null;
67  if (!meters) return coords.slice();
68
69  const total = lineLengthMeters(coords);
70  // Need something meaningful left in the middle, not two points a
71  // meter apart that still pinpoint where the trim happened.
72  if (total <= meters * 2 + 50) return null;
73
74  const start = positionAtDistance(coords, meters);
75  const end = positionAtDistance(coords, total - meters);
76  if (!start || !end) return null;
77
78  const kept = coords.filter(
79    (_, i) => i > start.index && i <= end.index
80  );
81
82  return [start.position, ...kept, end.position];
83}
84
85/**
86 * Walk `meters` along the line and return { position, index } where
87 * `index` is the last vertex at or before the returned position.
88 *
89 * Also used for "drop a pin anywhere along my track" — the participant
90 * scrubs a slider and this resolves it to a real coordinate.
91 */
92export function positionAtDistance(coords, meters) {
93  if (!Array.isArray(coords) || coords.length < 2) return null;
94  if (meters <= 0) return { position: coords[0].slice(0, 2), index: 0 };
95
96  let travelled = 0;
97  for (let i = 1; i < coords.length; i++) {
98    const seg = haversineMeters(coords[i - 1], coords[i]);
99    if (travelled + seg >= meters) {
100      // Linear interpolation across the segment. Over a GPS sample
101      // interval the great-circle and the straight line differ by far
102      // less than the fix accuracy, so this is not worth doing properly.
103      const t = seg === 0 ? 0 : (meters - travelled) / seg;
104      const [x1, y1] = coords[i - 1];
105      const [x2, y2] = coords[i];
106      return {
107        position: [x1 + (x2 - x1) * t, y1 + (y2 - y1) * t],
108        index: i - 1
109      };
110    }
111    travelled += seg;
112  }
113
114  const last = coords[coords.length - 1];
115  return { position: last.slice(0, 2), index: coords.length - 1 };
116}
117
118/**
119 * Thin a dense track by dropping points closer than `minMeters` to the
120 * last point kept. First and last are always kept.
121 *
122 * A watch recording at one-second intervals produces thousands of points
123 * for a short ride, most of them within a meter of each other. This is
124 * deliberately not Ramer-Douglas-Peucker: RDP is shape-aware and better,
125 * but it needs the whole line in memory at once and is markedly slower on
126 * a phone. Distance thinning is one pass and good enough for a track that
127 * is about to be drawn as a route.
128 */
129export function thinByDistance(coords, minMeters = 10) {
130  if (!Array.isArray(coords) || coords.length < 3) return coords?.slice() || [];
131
132  const out = [coords[0]];
133  for (let i = 1; i < coords.length - 1; i++) {
134    if (haversineMeters(out[out.length - 1], coords[i]) >= minMeters) {
135      out.push(coords[i]);
136    }
137  }
138  out.push(coords[coords.length - 1]);
139  return out;
140}
141
142/**
143 * Flatten a LineString or MultiLineString geometry, or a Feature
144 * wrapping one, to a single coordinate array.
145 *
146 * MultiLineString parts are joined end to end. That is a real
147 * simplification: it invents a segment across each gap. Callers that
148 * measure a route can live with it, because an imported GPX splits into
149 * parts where the recording paused rather than where the rider teleported.
150 */
151export function lineCoords(input) {
152  const geom = input?.type === 'Feature' ? input.geometry : input;
153  if (!geom) return [];
154  if (geom.type === 'LineString') return geom.coordinates || [];
155  if (geom.type === 'MultiLineString') return (geom.coordinates || []).flat();
156  return [];
157}
158
159/**
160 * Project a point onto a segment. Returns the distance to it in meters
161 * and `t`, the fraction along the segment where the foot of the
162 * perpendicular lands, clamped to the segment ends.
163 *
164 * Projected flat around the point before the math. Over the distances
165 * this is asked about, a corridor width or the offset of a campground
166 * from a road, the error from ignoring curvature is centimeters, and
167 * doing it this way keeps the inner loop to arithmetic.
168 */
169export function projectPointToSegment(p, a, b) {
170  const kx = M_PER_DEG_LAT * Math.cos((p[1] * Math.PI) / 180);
171  const ky = M_PER_DEG_LAT;
172
173  const ax = (a[0] - p[0]) * kx, ay = (a[1] - p[1]) * ky;
174  const bx = (b[0] - p[0]) * kx, by = (b[1] - p[1]) * ky;
175
176  const dx = bx - ax, dy = by - ay;
177  const lenSq = dx * dx + dy * dy;
178
179  // A zero-length segment has no direction to project onto, so the foot
180  // is the segment itself.
181  let t = lenSq === 0 ? 0 : ((0 - ax) * dx + (0 - ay) * dy) / lenSq;
182  t = Math.max(0, Math.min(1, t));
183
184  return { meters: Math.hypot(ax + t * dx, ay + t * dy), t };
185}
186
187/**
188 * The piece of a line between two distances along it, in meters.
189 *
190 * Both ends are interpolated rather than snapped to the nearest vertex.
191 * That is the whole point where this is used: cutting a tour into days
192 * at the overnight stops. An imported route can carry a kilometer
193 * between vertices, and snapping would move a day's end by that much.
194 *
195 * Returns null rather than a degenerate line if the two distances do not
196 * enclose anything.
197 */
198export function sliceByDistance(coords, startMeters, endMeters) {
199  if (!Array.isArray(coords) || coords.length < 2) return null;
200
201  const total = lineLengthMeters(coords);
202  const from = Math.max(0, Math.min(startMeters, total));
203  const to   = Math.max(0, Math.min(endMeters, total));
204  if (to - from <= 0) return null;
205
206  const start = positionAtDistance(coords, from);
207  const end   = positionAtDistance(coords, to);
208  if (!start || !end) return null;
209
210  const out = [start.position];
211  for (let i = start.index + 1; i <= end.index; i++) {
212    // A cut that lands exactly on a vertex would otherwise repeat it.
213    if (haversineMeters(out[out.length - 1], coords[i]) > 0.01) {
214      out.push(coords[i].slice(0, 2));
215    }
216  }
217  if (haversineMeters(out[out.length - 1], end.position) > 0.01) {
218    out.push(end.position);
219  }
220
221  return out.length >= 2 ? out : null;
222}
223
224export const metersToMiles = m => m / 1609.344;
225export const metersToFeet = m => m * 3.28084;