Skip to content

Commit d5f69eb

Browse files
committed
fix(tides): robust extrema detection (despike + hysteresis) for gauge data
Observation-derived tide events were computed with a ±N-sample local-extremum scan using non-strict comparisons, so raw gauge noise produced hundreds of spurious highs/lows and flat plateaus registered as both a high and a low at the same time. A single sentinel/spike (e.g. IOC's -1.0 m glitch) also inflated the range and suppressed the real extremum. Add a shared tideExtrema util: despikeSeries (rolling-median outlier removal) + findTideExtrema (hysteresis/zigzag detection with a min-prominence delta, strictly alternating H/L, edge-aware so a recent peak near "now" is still reported). Wire it into the IOC and Pegelonline providers, replacing their duplicated naive logic.
1 parent 61a5ad9 commit d5f69eb

5 files changed

Lines changed: 267 additions & 69 deletions

File tree

integrations/knowledge-tides-ioc/index.ts

Lines changed: 25 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,10 @@
1-
import { createPlace, type Place, USER_AGENT } from "@openmapx/core";
1+
import {
2+
createPlace,
3+
despikeSeries,
4+
findTideExtrema,
5+
type Place,
6+
USER_AGENT,
7+
} from "@openmapx/core";
28
import type { IntegrationContext } from "@openmapx/integration-framework";
39
import { registerPlaceResolver } from "@openmapx/place-ids";
410

@@ -168,38 +174,17 @@ async function fetchObservations(stationCode: string): Promise<IocDataPoint[]> {
168174
return best.sort((a, b) => a.stime.localeCompare(b.stime));
169175
}
170176

171-
/** Detect local extrema in a sorted observation series to derive past H/L events. */
177+
/**
178+
* Derive past H/L events from the observation series (metres). Uses the shared
179+
* hysteresis detector so sensor noise / flat plateaus don't produce spurious
180+
* extrema. Threshold: 3 cm or 12 % of the observed range, whichever is larger.
181+
*/
172182
function deriveExtrema(curve: Array<{ time: string; value: number }>): TideEvent[] {
173-
if (curve.length < 5) return [];
174-
const events: TideEvent[] = [];
175-
// Use a small window to skip noise; events must be a strict local max/min
176-
// over ±2 samples (~10–60 min depending on station's XMtInt).
177-
const W = 2;
178-
for (let i = W; i < curve.length - W; i++) {
179-
const v = curve[i].value;
180-
let isHigh = true;
181-
let isLow = true;
182-
for (let j = i - W; j <= i + W; j++) {
183-
if (j === i) continue;
184-
if (curve[j].value > v) isHigh = false;
185-
if (curve[j].value < v) isLow = false;
186-
if (!isHigh && !isLow) break;
187-
}
188-
if (isHigh) {
189-
events.push({
190-
time: curve[i].time,
191-
type: "H",
192-
valueFt: Math.round(v * M_TO_FT * 100) / 100,
193-
});
194-
} else if (isLow) {
195-
events.push({
196-
time: curve[i].time,
197-
type: "L",
198-
valueFt: Math.round(v * M_TO_FT * 100) / 100,
199-
});
200-
}
201-
}
202-
return events;
183+
return findTideExtrema(curve, { minDelta: 0.03, relativeDelta: 0.12 }).map((e) => ({
184+
time: e.time,
185+
type: e.type,
186+
valueFt: Math.round(e.value * M_TO_FT * 100) / 100,
187+
}));
203188
}
204189

205190
async function buildTidesResponse(
@@ -215,10 +200,14 @@ async function buildTidesResponse(
215200
const obs = await fetchObservations(station.code);
216201
if (obs.length === 0) return null;
217202

218-
const curveRaw = obs.map((p) => ({
219-
time: reformatIocTime(p.stime),
220-
value: p.slevel,
221-
}));
203+
// IOC raw gauge data carries occasional spikes / sentinels (e.g. an exact
204+
// -1.0 m glitch); drop them before plotting or deriving extrema, since one
205+
// outlier corrupts the range, threshold and event list. 0.25 m is far beyond
206+
// any real tide change over the ~15 s sampling interval.
207+
const curveRaw = despikeSeries(
208+
obs.map((p) => ({ time: reformatIocTime(p.stime), value: p.slevel })),
209+
0.25,
210+
);
222211
const curve = curveRaw.map((p) => ({
223212
time: p.time,
224213
valueFt: Math.round(p.value * M_TO_FT * 100) / 100,

integrations/knowledge-tides-pegelonline/index.ts

Lines changed: 24 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,10 @@
1-
import { createPlace, type Place, USER_AGENT } from "@openmapx/core";
1+
import {
2+
createPlace,
3+
despikeSeries,
4+
findTideExtrema,
5+
type Place,
6+
USER_AGENT,
7+
} from "@openmapx/core";
28
import type { IntegrationContext } from "@openmapx/integration-framework";
39
import { registerPlaceResolver } from "@openmapx/place-ids";
410

@@ -147,35 +153,18 @@ async function fetchMeasurements(uuid: string): Promise<PegelMeasurement[]> {
147153
return raw.filter((_, i) => i % 15 === 0);
148154
}
149155

156+
/**
157+
* Derive H/L events from the water-level series (cm). Uses the shared
158+
* hysteresis detector so sensor noise / flat plateaus don't produce spurious
159+
* extrema. Threshold: 5 cm or 12 % of the observed range, whichever is larger.
160+
*/
150161
function deriveExtrema(curve: Array<{ time: string; valueCm: number }>): TideEvent[] {
151-
if (curve.length < 5) return [];
152-
const events: TideEvent[] = [];
153-
const W = 3;
154-
for (let i = W; i < curve.length - W; i++) {
155-
const v = curve[i].valueCm;
156-
let isHigh = true;
157-
let isLow = true;
158-
for (let j = i - W; j <= i + W; j++) {
159-
if (j === i) continue;
160-
if (curve[j].valueCm > v) isHigh = false;
161-
if (curve[j].valueCm < v) isLow = false;
162-
if (!isHigh && !isLow) break;
163-
}
164-
if (isHigh) {
165-
events.push({
166-
time: curve[i].time,
167-
type: "H",
168-
valueFt: Math.round(v * CM_TO_FT * 100) / 100,
169-
});
170-
} else if (isLow) {
171-
events.push({
172-
time: curve[i].time,
173-
type: "L",
174-
valueFt: Math.round(v * CM_TO_FT * 100) / 100,
175-
});
176-
}
177-
}
178-
return events;
162+
const samples = curve.map((p) => ({ time: p.time, value: p.valueCm }));
163+
return findTideExtrema(samples, { minDelta: 5, relativeDelta: 0.12 }).map((e) => ({
164+
time: e.time,
165+
type: e.type,
166+
valueFt: Math.round(e.value * CM_TO_FT * 100) / 100,
167+
}));
179168
}
180169

181170
async function buildTidesResponse(
@@ -191,10 +180,12 @@ async function buildTidesResponse(
191180
const obs = await fetchMeasurements(station.uuid);
192181
if (obs.length === 0) return null;
193182

194-
const curveRaw = obs.map((p) => ({
195-
time: reformatPegelTime(p.timestamp),
196-
valueCm: p.value,
197-
}));
183+
// Drop spike/sentinel outliers (>25 cm from the local median) before plotting
184+
// or deriving extrema — one bad reading corrupts the range, threshold and events.
185+
const curveRaw = despikeSeries(
186+
obs.map((p) => ({ time: reformatPegelTime(p.timestamp), value: p.value })),
187+
25,
188+
).map((p) => ({ time: p.time, valueCm: p.value }));
198189
const curve = curveRaw.map((p) => ({
199190
time: p.time,
200191
valueFt: Math.round(p.valueCm * CM_TO_FT * 100) / 100,

packages/core/src/index.ts

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -415,6 +415,8 @@ export { pointInIsochroneGeometry } from "./utils/pointInPolygon";
415415
export { decodePolyline, encodePolyline } from "./utils/polyline";
416416
export { searchResultToCategoryPlace } from "./utils/searchResultToCategoryPlace";
417417
export { sectionSlug } from "./utils/sectionSlug";
418+
export type { TideExtremaOptions, TideExtreme, TideSample } from "./utils/tideExtrema";
419+
export { despikeSeries, findTideExtrema } from "./utils/tideExtrema";
418420
export {
419421
USER_AGENT,
420422
USER_AGENT_ADMIN,
Lines changed: 87 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,87 @@
1+
import { describe, expect, it } from "vitest";
2+
import { despikeSeries, findTideExtrema } from "./tideExtrema";
3+
4+
// Synthesize a sampled sine with optional noise. `t` in hours.
5+
function series(
6+
hours: number,
7+
stepMin: number,
8+
amplitude: number,
9+
periodH: number,
10+
noise = 0,
11+
): { time: string; value: number }[] {
12+
const out: { time: string; value: number }[] = [];
13+
for (let m = 0; m <= hours * 60; m += stepMin) {
14+
const t = m / 60;
15+
const wobble = noise ? (((m * 9301 + 49297) % 233280) / 233280 - 0.5) * 2 * noise : 0;
16+
out.push({
17+
time: new Date(Date.UTC(2026, 0, 1, 0, m)).toISOString(),
18+
value: amplitude * Math.sin((2 * Math.PI * t) / periodH) + wobble,
19+
});
20+
}
21+
return out;
22+
}
23+
24+
describe("findTideExtrema", () => {
25+
it("finds the expected number of alternating extrema on a clean semidiurnal curve", () => {
26+
// 24h, 12.42h period → ~2 highs + ~2 lows; amplitude 1m.
27+
const extrema = findTideExtrema(series(24, 6, 1, 12.42), {
28+
minDelta: 0.03,
29+
relativeDelta: 0.1,
30+
});
31+
// ~2 cycles over 24h → a handful of extrema, not dozens.
32+
expect(extrema.length).toBeGreaterThanOrEqual(2);
33+
expect(extrema.length).toBeLessThanOrEqual(6);
34+
// strictly alternating
35+
for (let i = 1; i < extrema.length; i++) {
36+
expect(extrema[i].type).not.toBe(extrema[i - 1].type);
37+
}
38+
});
39+
40+
it("rejects sensor noise instead of reporting dozens of micro-extrema", () => {
41+
// Same tide with 2cm noise sampled every minute (the bug's conditions).
42+
const noisy = series(24, 1, 1, 12.42, 0.02);
43+
const extrema = findTideExtrema(noisy, { minDelta: 0.03, relativeDelta: 0.1 });
44+
expect(extrema.length).toBeLessThanOrEqual(5);
45+
for (let i = 1; i < extrema.length; i++) {
46+
expect(extrema[i].type).not.toBe(extrema[i - 1].type);
47+
}
48+
});
49+
50+
it("emits a recent peak that sits near the end of the window (trailing extreme)", () => {
51+
// Descend to a trough, rise to a peak, then dip slightly — like a station's
52+
// past-24h ending just after a midday high (the Trieste/Poreč case).
53+
const pts: { time: string; value: number }[] = [];
54+
let m = 0;
55+
const push = (value: number) => {
56+
pts.push({ time: new Date(Date.UTC(2026, 0, 1, 0, m)).toISOString(), value });
57+
m += 10;
58+
};
59+
for (let v = 0.5; v >= -1; v -= 0.1) push(Number(v.toFixed(2))); // down to trough
60+
for (let v = -0.9; v <= 1; v += 0.1) push(Number(v.toFixed(2))); // up to peak
61+
push(0.95);
62+
push(0.9); // small post-peak dip at the very end
63+
const extrema = findTideExtrema(pts, { minDelta: 0.03, relativeDelta: 0.12 });
64+
expect(extrema.map((e) => e.type)).toEqual(["L", "H"]);
65+
});
66+
67+
it("despikeSeries removes a sentinel outlier so it doesn't spawn spurious extrema", () => {
68+
const clean = series(24, 6, 1, 12.42);
69+
const withSpike = clean.map((s, i) => (i === 20 ? { ...s, value: -1 } : s));
70+
const before = findTideExtrema(withSpike, { minDelta: 0.03, relativeDelta: 0.12 });
71+
const after = findTideExtrema(despikeSeries(withSpike, 0.25), {
72+
minDelta: 0.03,
73+
relativeDelta: 0.12,
74+
});
75+
// The spike injects spurious H/L around index 20; despiking removes them.
76+
expect(after.length).toBeLessThan(before.length);
77+
expect(despikeSeries(withSpike, 0.25)).toHaveLength(withSpike.length - 1);
78+
});
79+
80+
it("never emits a high and a low at the same flat plateau", () => {
81+
const flat = Array.from({ length: 60 }, (_, m) => ({
82+
time: new Date(Date.UTC(2026, 0, 1, 0, m)).toISOString(),
83+
value: -0.7,
84+
}));
85+
expect(findTideExtrema(flat, { minDelta: 0.03, relativeDelta: 0.1 })).toEqual([]);
86+
});
87+
});
Lines changed: 129 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,129 @@
1+
export interface TideSample {
2+
time: string;
3+
value: number;
4+
}
5+
6+
export interface TideExtreme {
7+
time: string;
8+
value: number;
9+
type: "H" | "L";
10+
}
11+
12+
export interface TideExtremaOptions {
13+
/** Minimum swing (in the samples' own units) below which a reversal is noise. */
14+
minDelta: number;
15+
/** Extra threshold as a fraction of the series range, added via `max`. Default 0. */
16+
relativeDelta?: number;
17+
}
18+
19+
/**
20+
* Derive tide highs/lows from a noisy observation series using hysteresis
21+
* (a "zigzag" peak detector). A reversal is only registered once the series
22+
* has moved at least `delta = max(minDelta, range * relativeDelta)` away from
23+
* the running peak/trough.
24+
*
25+
* This replaces a naive ±N-sample local-extremum scan, which on raw gauge data
26+
* reports dozens of spurious highs/lows — and, because it compared
27+
* non-strictly, flagged flat plateaus as both a high AND a low at the same
28+
* time. The hysteresis output is noise-free and strictly alternates H/L.
29+
*/
30+
function median(values: number[]): number {
31+
const sorted = [...values].sort((a, b) => a - b);
32+
const mid = sorted.length >> 1;
33+
return sorted.length % 2 ? sorted[mid] : (sorted[mid - 1] + sorted[mid]) / 2;
34+
}
35+
36+
/**
37+
* Drop spike/sentinel outliers from a gauge series: a sample is removed when it
38+
* deviates from the median of a centred ±`window` neighbourhood by more than
39+
* `maxDeviation` (in the samples' own units). Tide levels change slowly, so a
40+
* large jump over a few seconds is a sensor glitch (e.g. IOC's `-1.0` sentinel)
41+
* — and one such outlier wrecks both the range/threshold and the extrema. Run
42+
* this before {@link findTideExtrema} and before plotting the curve.
43+
*/
44+
export function despikeSeries(
45+
samples: TideSample[],
46+
maxDeviation: number,
47+
window = 5,
48+
): TideSample[] {
49+
if (samples.length <= window * 2) return samples;
50+
const values = samples.map((s) => s.value);
51+
const kept: TideSample[] = [];
52+
for (let i = 0; i < samples.length; i++) {
53+
const lo = Math.max(0, i - window);
54+
const hi = Math.min(values.length, i + window + 1);
55+
if (Math.abs(values[i] - median(values.slice(lo, hi))) <= maxDeviation) kept.push(samples[i]);
56+
}
57+
return kept;
58+
}
59+
60+
export function findTideExtrema(samples: TideSample[], opts: TideExtremaOptions): TideExtreme[] {
61+
if (samples.length < 5) return [];
62+
63+
let lo = Number.POSITIVE_INFINITY;
64+
let hi = Number.NEGATIVE_INFINITY;
65+
for (const s of samples) {
66+
if (s.value < lo) lo = s.value;
67+
if (s.value > hi) hi = s.value;
68+
}
69+
const delta = Math.max(opts.minDelta, (hi - lo) * (opts.relativeDelta ?? 0));
70+
if (!(delta > 0)) return [];
71+
72+
const extrema: TideExtreme[] = [];
73+
let mn = Number.POSITIVE_INFINITY;
74+
let mx = Number.NEGATIVE_INFINITY;
75+
let mnIdx = 0;
76+
let mxIdx = 0;
77+
// null until the first significant move decides whether to open on a high or low.
78+
let lookForMax: boolean | null = null;
79+
// Value of the last confirmed extreme, to gauge whether the trailing swing is real.
80+
let pivot: number | null = null;
81+
82+
for (let i = 0; i < samples.length; i++) {
83+
const v = samples[i].value;
84+
if (v > mx) {
85+
mx = v;
86+
mxIdx = i;
87+
}
88+
if (v < mn) {
89+
mn = v;
90+
mnIdx = i;
91+
}
92+
93+
if (lookForMax !== false && v < mx - delta) {
94+
// Skip index 0: a drop from the first sample means the window opened past
95+
// a peak — a window edge, not a real turning point. Still transition.
96+
if (mxIdx > 0) {
97+
extrema.push({ time: samples[mxIdx].time, value: samples[mxIdx].value, type: "H" });
98+
}
99+
pivot = samples[mxIdx].value;
100+
mn = v;
101+
mnIdx = i;
102+
lookForMax = false;
103+
} else if (lookForMax !== true && v > mn + delta) {
104+
if (mnIdx > 0) {
105+
extrema.push({ time: samples[mnIdx].time, value: samples[mnIdx].value, type: "L" });
106+
}
107+
pivot = samples[mnIdx].value;
108+
mx = v;
109+
mxIdx = i;
110+
lookForMax = true;
111+
}
112+
}
113+
114+
// Emit the final turning point. Hysteresis only confirms an extreme once the
115+
// series reverses past `delta`, so a peak/trough sitting at the very end of
116+
// the window (e.g. a recent high near "now") would otherwise be dropped.
117+
// Require it to be an actual turn (not the last sample) and ≥ delta from the
118+
// last pivot, so a still-rising/falling tail isn't reported as an extreme.
119+
const lastIdx = samples.length - 1;
120+
if (pivot !== null) {
121+
if (lookForMax && mxIdx < lastIdx && samples[mxIdx].value - pivot >= delta) {
122+
extrema.push({ time: samples[mxIdx].time, value: samples[mxIdx].value, type: "H" });
123+
} else if (lookForMax === false && mnIdx < lastIdx && pivot - samples[mnIdx].value >= delta) {
124+
extrema.push({ time: samples[mnIdx].time, value: samples[mnIdx].value, type: "L" });
125+
}
126+
}
127+
128+
return extrema;
129+
}

0 commit comments

Comments
 (0)