From 990bbaa99e8d9e2ba5483e016eb92a419990c01f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ullrich=20Sch=C3=A4fer?= Date: Tue, 14 Jul 2026 00:46:13 +0200 Subject: [PATCH] feat(gpx): elevation cleaning primitives (despike + hysteresis totals) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Task group 1 of elevation-profile-hardening. New pure module elevation-clean.ts (not yet wired — group 2 does that): - despike(points): interpolates out isolated peak/pit spikes (opposite-sign slope outliers both exceeding MAX_SLOPE_PERCENT=100); steep monotonic terrain (same-sign) is never touched; per-segment, no cross-gap interp. - filteredTotals(points): hysteresis ascent/descent (NOISE_THRESHOLD_M=5) that suppresses sub-threshold jitter, returning {gain,loss,gainRaw,lossRaw}. - cumulativeFilteredTotals(points): running filtered arrays for per-day splits (day total = cumulative[end] − cumulative[start]). DESIGN NOTE for review: the spec's task 1.2 says the zero-filtered fallback should "report raw", but the raw sum-of-deltas re-inflates the very jitter we remove (a ±2 m wobble sums to several metres) — directly contradicting task 1.4's "±2 m jitter → 0 totals". Resolved toward the stated goal: fall back to NET change (end−start), so flat wobble reads ~0 while a genuine small net climb still surfaces. gainRaw/lossRaw still expose the raw sums. Verified: gpx typecheck+lint clean, 108/108; full pnpm test 11/11. Co-Authored-By: Claude Opus 4.8 (1M context) --- .../elevation-profile-hardening/tasks.md | 8 +- packages/gpx/src/elevation-clean.test.ts | 88 ++++++++++ packages/gpx/src/elevation-clean.ts | 155 ++++++++++++++++++ 3 files changed, 247 insertions(+), 4 deletions(-) create mode 100644 packages/gpx/src/elevation-clean.test.ts create mode 100644 packages/gpx/src/elevation-clean.ts diff --git a/openspec/changes/elevation-profile-hardening/tasks.md b/openspec/changes/elevation-profile-hardening/tasks.md index 38e4c58..a106214 100644 --- a/openspec/changes/elevation-profile-hardening/tasks.md +++ b/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 -- [ ] 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) -- [ ] 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.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.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.3 Add `cumulativeFilteredTotals(points, opts?)` returning per-point running filtered ascent/descent arrays (for per-day splits) +- [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 diff --git a/packages/gpx/src/elevation-clean.test.ts b/packages/gpx/src/elevation-clean.test.ts new file mode 100644 index 0000000..ee3ed0b --- /dev/null +++ b/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]!); + } + }); +}); diff --git a/packages/gpx/src/elevation-clean.ts b/packages/gpx/src/elevation-clean.ts new file mode 100644 index 0000000..7331e12 --- /dev/null +++ b/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 }; +}