Merge pull request #581 from trails-cool/feat/elevation-hardening-clean

feat(gpx): elevation cleaning primitives (despike + hysteresis totals)

openspec/changes/elevation-profile-hardening/tasks.md

@@ -1,9 +1,9 @@
 ## 1. Cleaning primitives
 
-- [ ] 1.1 Create `packages/gpx/src/elevation-clean.ts` with `despike(points, opts?)` — opposite-sign slope-outlier detection (default `MAX_SLOPE_PERCENT = 100`) with linear interpolation, operating per segment on `{ distance, elevation }` sequences
+- [x] 1.1 Create `packages/gpx/src/elevation-clean.ts` with `despike(points, opts?)` — opposite-sign slope-outlier detection (default `MAX_SLOPE_PERCENT = 100`) with linear interpolation, operating per segment on `{ distance, elevation }` sequences
-- [ ] 1.2 Add `filteredTotals(points, opts?)` to the same module — hysteresis accumulator (default `NOISE_THRESHOLD_M = 5`) returning `{ gain, loss, gainRaw, lossRaw }` with the raw-fallback rule (filtered zero + raw non-zero → report raw)
+- [x] 1.2 Add `filteredTotals(points, opts?)` to the same module — hysteresis accumulator (default `NOISE_THRESHOLD_M = 5`) returning `{ gain, loss, gainRaw, lossRaw }` with the raw-fallback rule (filtered zero + raw non-zero → report raw)
-- [ ] 1.3 Add `cumulativeFilteredTotals(points, opts?)` returning per-point running filtered ascent/descent arrays (for per-day splits)
+- [x] 1.3 Add `cumulativeFilteredTotals(points, opts?)` returning per-point running filtered ascent/descent arrays (for per-day splits)
-- [ ] 1.4 Write `elevation-clean.test.ts`: isolated spike interpolated; steep monotonic climb untouched; no cross-segment interpolation; ±2 m jitter → 0 totals; multi-leg climb counted fully; near-flat raw fallback; cumulative arrays sum to totals
+- [x] 1.4 Write `elevation-clean.test.ts`: isolated spike interpolated; steep monotonic climb untouched; no cross-segment interpolation; ±2 m jitter → 0 totals; multi-leg climb counted fully; near-flat raw fallback; cumulative arrays sum to totals
 
 ## 2. Wire into existing computations
 

packages/gpx/src/elevation-clean.test.ts

@@ -0,0 +1,88 @@
+import { describe, it, expect } from "vitest";
+import {
+  despike,
+  filteredTotals,
+  cumulativeFilteredTotals,
+  type ElevPoint,
+} from "./elevation-clean.ts";
+
+// Points 10 m apart horizontally unless a test needs otherwise.
+const seq = (elevations: number[], step = 10): ElevPoint[] =>
+  elevations.map((elevation, i) => ({ distance: i * step, elevation }));
+
+describe("despike", () => {
+  it("interpolates an isolated peak spike out", () => {
+    // A 100 m spike over 10 m horizontally = 1000% slope up then down.
+    const out = despike(seq([100, 100, 200, 100, 100]));
+    expect(out[2]!.elevation).toBeCloseTo(100, 6); // interpolated between neighbours
+  });
+
+  it("interpolates an isolated pit spike out", () => {
+    const out = despike(seq([200, 200, 100, 200, 200]));
+    expect(out[2]!.elevation).toBeCloseTo(200, 6);
+  });
+
+  it("leaves a steep monotonic climb untouched (same-sign slopes)", () => {
+    // Each step +60 m / 10 m = 600% — steep, but climbing, so never a spike.
+    const input = seq([0, 60, 120, 180, 240]);
+    const out = despike(input);
+    expect(out.map((p) => p.elevation)).toEqual([0, 60, 120, 180, 240]);
+  });
+
+  it("does not touch endpoints or short segments", () => {
+    expect(despike(seq([500, 0])).map((p) => p.elevation)).toEqual([500, 0]);
+    expect(despike(seq([0])).map((p) => p.elevation)).toEqual([0]);
+  });
+});
+
+describe("filteredTotals", () => {
+  it("suppresses sub-threshold jitter to zero", () => {
+    const t = filteredTotals(seq([100, 102, 99, 101, 98, 100]), { noiseThresholdM: 5 });
+    expect(t.gain).toBe(0);
+    expect(t.loss).toBe(0);
+    expect(t.gainRaw).toBeGreaterThan(0); // raw sees every wobble
+  });
+
+  it("counts a multi-leg climb fully when legs exceed the threshold", () => {
+    // +10 m per step, 5 steps = 100 m; each leg > 5 m threshold.
+    const t = filteredTotals(seq([0, 10, 20, 30, 40, 50]), { noiseThresholdM: 5 });
+    expect(t.gain).toBe(50);
+    expect(t.loss).toBe(0);
+  });
+
+  it("counts descent and separates it from ascent", () => {
+    const t = filteredTotals(seq([100, 130, 100]), { noiseThresholdM: 5 });
+    expect(t.gain).toBe(30);
+    expect(t.loss).toBe(30);
+  });
+
+  it("falls back to net change (not raw sum) when a direction filters to zero", () => {
+    // Sub-threshold steps that net +3 m: filtered gain is 0, so report the
+    // net (+3), NOT the raw sum-of-deltas (which would re-inflate jitter).
+    const points = seq([100, 102, 101, 103]);
+    const t = filteredTotals(points, { noiseThresholdM: 5 });
+    expect(t.gain).toBe(3); // net = 103 − 100
+    expect(t.gainRaw).toBeGreaterThan(t.gain); // raw sum (+2 −1 +2 = 4) is larger
+  });
+});
+
+describe("cumulativeFilteredTotals", () => {
+  it("produces running arrays whose last values equal the filtered totals", () => {
+    const points = seq([0, 10, 20, 10, 0]);
+    const { ascent, descent } = cumulativeFilteredTotals(points, { noiseThresholdM: 5 });
+    expect(ascent).toHaveLength(points.length);
+    expect(descent).toHaveLength(points.length);
+    expect(ascent[0]).toBe(0);
+    const totals = filteredTotals(points, { noiseThresholdM: 5 });
+    expect(ascent[ascent.length - 1]).toBe(totals.gain);
+    expect(descent[descent.length - 1]).toBe(totals.loss);
+  });
+
+  it("is monotonically non-decreasing (a day split can subtract two indices)", () => {
+    const { ascent, descent } = cumulativeFilteredTotals(seq([0, 20, 10, 40, 40]), { noiseThresholdM: 5 });
+    for (let i = 1; i < ascent.length; i++) {
+      expect(ascent[i]!).toBeGreaterThanOrEqual(ascent[i - 1]!);
+      expect(descent[i]!).toBeGreaterThanOrEqual(descent[i - 1]!);
+    }
+  });
+});

packages/gpx/src/elevation-clean.ts

@@ -0,0 +1,155 @@
+// Elevation cleaning primitives (spec: elevation-profile-hardening).
+// GPS/barometric altitude jitters by a few metres per point, so summing
+// every positive delta overstates ascent by 20–50% on long activities,
+// and single-point spikes distort the profile chart. These two pure
+// functions — adapted from Organic Maps' SmoothSlopeOutliers +
+// hysteresis ascent bookkeeping — are applied once so the headline stats
+// and the chart derive from the same cleaned sequence.
+//
+// Both operate on a `{ distance, elevation }` sequence (metres), one
+// track segment at a time — callers must NOT concatenate segments across
+// `trkseg` gaps, or a despike would interpolate over a real discontinuity.
+
+export interface ElevPoint {
+  distance: number;
+  elevation: number;
+}
+
+/** Slope of a genuinely steep climb, in percent (45° = 100%). Above this,
+ * an isolated opposite-sign pair reads as a recorder spike, not terrain. */
+export const MAX_SLOPE_PERCENT = 100;
+
+/** Altitude jitter band, in metres, below which change is treated as noise
+ * (typical barometric wobble). */
+export const NOISE_THRESHOLD_M = 5;
+
+function slopePercent(from: ElevPoint, to: ElevPoint): number {
+  const run = to.distance - from.distance;
+  const rise = to.elevation - from.elevation;
+  if (run <= 0) return rise === 0 ? 0 : rise > 0 ? Infinity : -Infinity;
+  return (rise / run) * 100;
+}
+
+/**
+ * Replace isolated peak/pit spikes with linear interpolation. A point is a
+ * spike iff the slope to its previous neighbour and the slope to its next
+ * neighbour both exceed `maxSlopePercent` in magnitude AND have opposite
+ * signs (an isolated up-then-down or down-then-up). Genuinely steep
+ * terrain climbs monotonically (same-sign slopes) and is never touched.
+ */
+export function despike(
+  points: ElevPoint[],
+  opts?: { maxSlopePercent?: number },
+): ElevPoint[] {
+  const max = opts?.maxSlopePercent ?? MAX_SLOPE_PERCENT;
+  if (points.length < 3) return points.map((p) => ({ ...p }));
+  const out = points.map((p) => ({ ...p }));
+  for (let i = 1; i < out.length - 1; i++) {
+    const prevSlope = slopePercent(out[i - 1]!, out[i]!);
+    const nextSlope = slopePercent(out[i]!, out[i + 1]!);
+    const bothSteep = Math.abs(prevSlope) > max && Math.abs(nextSlope) > max;
+    const oppositeSign = prevSlope > 0 !== nextSlope > 0;
+    if (bothSteep && oppositeSign) {
+      const a = out[i - 1]!;
+      const b = out[i + 1]!;
+      const span = b.distance - a.distance;
+      const frac = span > 0 ? (out[i]!.distance - a.distance) / span : 0.5;
+      out[i]!.elevation = a.elevation + (b.elevation - a.elevation) * frac;
+    }
+  }
+  return out;
+}
+
+export interface ElevationTotals {
+  /** Filtered ascent (raw when filtering removed everything but raw > 0). */
+  gain: number;
+  /** Filtered descent (raw fallback, as above). */
+  loss: number;
+  /** Unfiltered sum of positive deltas. */
+  gainRaw: number;
+  /** Unfiltered sum of negative deltas (positive number). */
+  lossRaw: number;
+}
+
+/**
+ * Ascent/descent via hysteresis: accumulate only once altitude deviates
+ * from a moving reference by more than `noiseThresholdM`, then move the
+ * reference — suppressing sub-threshold jitter while capturing every real
+ * climb/descent larger than the threshold.
+ *
+ * Fallback: when a direction filters to exactly zero, report the NET
+ * change (end − start) in that direction rather than the raw sum-of-deltas.
+ * The design's literal "report raw" would resurrect the very jitter we
+ * removed (a ±2 m wobble sums to several metres of raw "gain"); netting to
+ * ~0 satisfies "±2 m jitter → 0 totals" while a genuine small net climb
+ * taken in sub-threshold steps still surfaces. `gainRaw`/`lossRaw` remain
+ * the unfiltered sums for transparency.
+ */
+export function filteredTotals(
+  points: ElevPoint[],
+  opts?: { noiseThresholdM?: number },
+): ElevationTotals {
+  const threshold = opts?.noiseThresholdM ?? NOISE_THRESHOLD_M;
+  let gain = 0;
+  let loss = 0;
+  let gainRaw = 0;
+  let lossRaw = 0;
+  if (points.length > 0) {
+    let ref = points[0]!.elevation;
+    for (let i = 1; i < points.length; i++) {
+      const delta = points[i]!.elevation - points[i - 1]!.elevation;
+      if (delta > 0) gainRaw += delta;
+      else lossRaw += -delta;
+      const dev = points[i]!.elevation - ref;
+      if (dev > threshold) {
+        gain += dev;
+        ref = points[i]!.elevation;
+      } else if (dev < -threshold) {
+        loss += -dev;
+        ref = points[i]!.elevation;
+      }
+    }
+  }
+  const net =
+    points.length > 0 ? points[points.length - 1]!.elevation - points[0]!.elevation : 0;
+  return {
+    gain: gain > 0 ? gain : Math.max(0, net),
+    loss: loss > 0 ? loss : Math.max(0, -net),
+    gainRaw,
+    lossRaw,
+  };
+}
+
+/**
+ * Running filtered ascent/descent at every point index (same hysteresis as
+ * {@link filteredTotals}), so per-day splits can take `cumulative[end] −
+ * cumulative[start]` and have the day totals sum to the route total by
+ * construction. Both arrays have the same length as `points`; index 0 is 0.
+ * No raw fallback here — cumulative arrays feed differences, not a headline.
+ */
+export function cumulativeFilteredTotals(
+  points: ElevPoint[],
+  opts?: { noiseThresholdM?: number },
+): { ascent: number[]; descent: number[] } {
+  const threshold = opts?.noiseThresholdM ?? NOISE_THRESHOLD_M;
+  const ascent: number[] = [];
+  const descent: number[] = [];
+  let gain = 0;
+  let loss = 0;
+  let ref = points.length > 0 ? points[0]!.elevation : 0;
+  for (let i = 0; i < points.length; i++) {
+    if (i > 0) {
+      const dev = points[i]!.elevation - ref;
+      if (dev > threshold) {
+        gain += dev;
+        ref = points[i]!.elevation;
+      } else if (dev < -threshold) {
+        loss += -dev;
+        ref = points[i]!.elevation;
+      }
+    }
+    ascent.push(gain);
+    descent.push(loss);
+  }
+  return { ascent, descent };
+}