fix: DEM elevation overcounting and add hysteresis-only recalculation button

- dem.py: apply 45s median filter before hysteresis to suppress SRTM
  tile-boundary steps that were accumulating through the 5m threshold;
  raise DEM hysteresis threshold from 5m to 10m
- dem.py: back up elevation_m as elevation_m_original in timeseries
  before the first DEM overwrite, so original sensor data is preserved
- dem.py: add recalculate_elevation_hysteresis() — recomputes gain/loss
  from original recorded elevation (reads elevation_m_original if a DEM
  run already replaced elevation_m) using source-aware thresholds
  (5m barometric, 10m GPS/unknown); does not touch the elevation array
- edit/server.py, serve/server.py: split /recalculate-elevation into
  two endpoints: /recalculate-elevation/dem and
  /recalculate-elevation/hysteresis
- EditDrawer.svelte: replace single DEM button with two side-by-side
  buttons — "Recalculate (hysteresis)" (fast, offline) and
  "Recalculate (DEM)" (SRTM lookup)
This commit is contained in:
Davide Scaini
2026-04-20 21:41:23 +02:00
parent 2b7a37ed41
commit ebac3f50f4
4 changed files with 215 additions and 42 deletions
+19 -2
View File
@@ -428,8 +428,8 @@ async def save_activity(activity_id: str, payload: dict[str, Any]) -> JSONRespon
return JSONResponse({"ok": True, "sidecar": str(sidecar_path)}) return JSONResponse({"ok": True, "sidecar": str(sidecar_path)})
@app.post("/api/activity/{activity_id}/recalculate-elevation") @app.post("/api/activity/{activity_id}/recalculate-elevation/dem")
async def recalculate_elevation_endpoint(activity_id: str) -> JSONResponse: async def recalculate_elevation_dem_endpoint(activity_id: str) -> JSONResponse:
"""Replace GPS altitude with DEM terrain elevation and recompute gain/loss. """Replace GPS altitude with DEM terrain elevation and recompute gain/loss.
Requires --dem-url to be set when starting bincio edit. Requires --dem-url to be set when starting bincio edit.
@@ -450,6 +450,23 @@ async def recalculate_elevation_endpoint(activity_id: str) -> JSONResponse:
raise HTTPException(422, str(e)) raise HTTPException(422, str(e))
@app.post("/api/activity/{activity_id}/recalculate-elevation/hysteresis")
async def recalculate_elevation_hysteresis_endpoint(activity_id: str) -> JSONResponse:
"""Recompute gain/loss from original recorded elevation using source-aware hysteresis."""
dd = _get_data_dir()
_check_id(activity_id)
try:
from bincio.extract.dem import recalculate_elevation_hysteresis
from bincio.render.merge import merge_one
result = recalculate_elevation_hysteresis(dd, activity_id)
merge_one(dd, activity_id)
return JSONResponse(result)
except FileNotFoundError as e:
raise HTTPException(404, str(e))
except ValueError as e:
raise HTTPException(422, str(e))
@app.post("/api/activity/{activity_id}/images") @app.post("/api/activity/{activity_id}/images")
async def upload_image(activity_id: str, file: UploadFile = File(...)) -> JSONResponse: async def upload_image(activity_id: str, file: UploadFile = File(...)) -> JSONResponse:
dd = _get_data_dir() dd = _get_data_dir()
+144 -24
View File
@@ -13,6 +13,7 @@ Pass the base URL (without path) to bincio serve/edit via --dem-url.
from __future__ import annotations from __future__ import annotations
import json import json
import statistics
import urllib.error import urllib.error
import urllib.request import urllib.request
from pathlib import Path from pathlib import Path
@@ -26,9 +27,31 @@ _DEFAULT_SAMPLE_INTERVAL_S = 10
# Maximum locations per API request. Open-Elevation supports ~512 per POST. # Maximum locations per API request. Open-Elevation supports ~512 per POST.
_DEFAULT_BATCH_SIZE = 512 _DEFAULT_BATCH_SIZE = 512
# Hysteresis threshold after DEM correction. DEM data is already smooth # Hysteresis threshold after DEM correction.
# (30 m horizontal resolution) so 5 m is a generous dead-band. # SRTM30 at 1 Hz produces tile-boundary steps of ~13 m every few seconds.
_DEM_HYSTERESIS_M = 5.0 # A 5 m threshold barely filters them; 10 m suppresses them reliably.
_DEM_HYSTERESIS_M = 10.0
# Median filter window (seconds / samples at 1 Hz) applied to DEM-interpolated
# series before hysteresis. 45 s smooths SRTM tile steps while keeping real
# climbs (typical cycling ramp > 100 m over > 2 min).
_MEDIAN_WINDOW_S = 45
def _median_filter(values: list[float], window: int) -> list[float]:
"""Apply a sliding-window median filter to *values*.
The window is centred on each sample; edges are handled by shrinking the
window symmetrically (same as scipy's 'reflect' / 'nearest' default).
"""
half = window // 2
n = len(values)
out: list[float] = []
for i in range(n):
lo = max(0, i - half)
hi = min(n, i + half + 1)
out.append(statistics.median(values[lo:hi]))
return out
def lookup_elevations( def lookup_elevations(
@@ -77,6 +100,28 @@ def lookup_elevations(
return results return results
def _hysteresis_gain_loss(
elevations: list[float], threshold_m: float
) -> tuple[float, float]:
"""Compute elevation gain and loss using a hysteresis dead-band.
Only commits a new elevation level when it differs from the last committed
value by at least *threshold_m*. Returns (gain, loss) in metres, both
as positive numbers.
"""
gain = loss = 0.0
committed = elevations[0]
for e in elevations[1:]:
diff = e - committed
if abs(diff) >= threshold_m:
if diff > 0:
gain += diff
else:
loss += abs(diff)
committed = e
return gain, loss
def recalculate_elevation( def recalculate_elevation(
user_dir: Path, user_dir: Path,
activity_id: str, activity_id: str,
@@ -91,10 +136,13 @@ def recalculate_elevation(
2. Subsample GPS points every *sample_interval_s* seconds. 2. Subsample GPS points every *sample_interval_s* seconds.
3. Query the DEM API for those points (batched). 3. Query the DEM API for those points (batched).
4. Linearly interpolate DEM elevation back to every GPS-valid second. 4. Linearly interpolate DEM elevation back to every GPS-valid second.
5. Apply :data:`_DEM_HYSTERESIS_M` hysteresis to compute gain / loss. 5. Apply a 45 s median filter to smooth SRTM tile-boundary noise.
6. Write the updated ``elevation_m`` array to the timeseries JSON. 6. Apply :data:`_DEM_HYSTERESIS_M` (10 m) hysteresis to compute gain/loss.
7. Patch ``elevation_gain_m`` / ``elevation_loss_m`` in the detail JSON. 7. Preserve the original elevation as ``elevation_m_original`` in the
8. Patch ``elevation_gain_m`` in ``index.json`` (summary entry for the feed). timeseries (only on the first DEM run — never overwrites a prior backup).
8. Write the smoothed DEM elevation array as ``elevation_m``.
9. Patch ``elevation_gain_m`` / ``elevation_loss_m`` in the detail JSON.
10. Patch ``elevation_gain_m`` in ``index.json`` (summary entry for feed).
Returns Returns
------- -------
@@ -165,7 +213,6 @@ def recalculate_elevation(
j = 0 j = 0
for i in gps_idx: for i in gps_idx:
t = t_arr[i] t = t_arr[i]
# Advance j so valid_pairs[j] is the left anchor for t
while j + 1 < len(valid_pairs) - 1 and valid_pairs[j + 1][0] <= t: while j + 1 < len(valid_pairs) - 1 and valid_pairs[j + 1][0] <= t:
j += 1 j += 1
t0, e0 = valid_pairs[j] t0, e0 = valid_pairs[j]
@@ -176,33 +223,37 @@ def recalculate_elevation(
else: else:
new_ele[i] = round(e0, 1) new_ele[i] = round(e0, 1)
# ── 4. Hysteresis accumulation ──────────────────────────────────────────── # ── 4. Median filter to suppress SRTM tile-boundary steps ────────────────
valid_eles = [e for e in new_ele if e is not None] valid_indices = [i for i, e in enumerate(new_ele) if e is not None]
gain = loss = 0.0 valid_eles_raw = [new_ele[i] for i in valid_indices] # type: ignore[misc]
committed = valid_eles[0] smoothed = _median_filter(valid_eles_raw, _MEDIAN_WINDOW_S) # type: ignore[arg-type]
for e in valid_eles[1:]:
diff = e - committed # Write smoothed values back into new_ele
if abs(diff) >= _DEM_HYSTERESIS_M: for idx, e in zip(valid_indices, smoothed):
if diff > 0: new_ele[idx] = round(e, 1)
gain += diff
else: # ── 5. Hysteresis accumulation on smoothed series ─────────────────────────
loss += diff smoothed_valid = [e for e in new_ele if e is not None]
committed = e gain, loss = _hysteresis_gain_loss(smoothed_valid, _DEM_HYSTERESIS_M) # type: ignore[arg-type]
gain_r = round(gain, 1) gain_r = round(gain, 1)
loss_r = round(abs(loss), 1) loss_r = round(loss, 1)
# ── 5. Write timeseries ─────────────────────────────────────────────────── # ── 6. Preserve original elevation (only if not already backed up) ────────
if "elevation_m_original" not in ts:
ts["elevation_m_original"] = ts.get("elevation_m")
# ── 7. Write timeseries ───────────────────────────────────────────────────
ts["elevation_m"] = new_ele ts["elevation_m"] = new_ele
ts_path.write_text(json.dumps(ts, indent=2, ensure_ascii=False), encoding="utf-8") ts_path.write_text(json.dumps(ts, indent=2, ensure_ascii=False), encoding="utf-8")
# ── 6. Patch activity detail JSON ───────────────────────────────────────── # ── 8. Patch activity detail JSON ─────────────────────────────────────────
detail = json.loads(json_path.read_text(encoding="utf-8")) detail = json.loads(json_path.read_text(encoding="utf-8"))
detail["elevation_gain_m"] = gain_r detail["elevation_gain_m"] = gain_r
detail["elevation_loss_m"] = loss_r detail["elevation_loss_m"] = loss_r
json_path.write_text(json.dumps(detail, indent=2, ensure_ascii=False), encoding="utf-8") json_path.write_text(json.dumps(detail, indent=2, ensure_ascii=False), encoding="utf-8")
# ── 7. Patch index.json summary (so merge_all picks up the new value) ───── # ── 9. Patch index.json summary ───────────────────────────────────────────
index_path = user_dir / "index.json" index_path = user_dir / "index.json"
if index_path.exists(): if index_path.exists():
index = json.loads(index_path.read_text(encoding="utf-8")) index = json.loads(index_path.read_text(encoding="utf-8"))
@@ -219,3 +270,72 @@ def recalculate_elevation(
"elevation_loss_m": loss_r, "elevation_loss_m": loss_r,
"points_queried": n_queried, "points_queried": n_queried,
} }
def recalculate_elevation_hysteresis(user_dir: Path, activity_id: str) -> dict:
"""Recompute elevation gain/loss from the original recorded elevation data.
Uses the same source-aware hysteresis thresholds as the extract pipeline:
- 5 m for barometric altimeters (FIT files with ``enhanced_altitude``)
- 10 m for GPS-derived altitude (GPX, TCX, FIT without barometric)
The elevation array in the timeseries is **not** modified. If a DEM
correction was previously applied, the backup in ``elevation_m_original``
is used as the source so the original sensor data is recovered.
Returns
-------
dict with keys ``elevation_gain_m``, ``elevation_loss_m``.
"""
acts_dir = user_dir / "activities"
json_path = acts_dir / f"{activity_id}.json"
ts_path = acts_dir / f"{activity_id}.timeseries.json"
if not json_path.exists():
raise FileNotFoundError(f"Activity not found: {activity_id}")
if not ts_path.exists():
raise ValueError("Activity has no timeseries data")
ts = json.loads(ts_path.read_text(encoding="utf-8"))
# Use original elevation if a DEM backup exists, otherwise use current
ele_arr: list[Optional[float]] = (
ts.get("elevation_m_original") or ts.get("elevation_m") or []
)
elevations = [e for e in ele_arr if e is not None]
if len(elevations) < 2:
raise ValueError("Not enough elevation data to compute gain/loss")
# Determine threshold from altitude_source stored in detail JSON
detail = json.loads(json_path.read_text(encoding="utf-8"))
altitude_source = detail.get("altitude_source", "unknown")
threshold = 5.0 if altitude_source == "barometric" else 10.0
gain, loss = _hysteresis_gain_loss(elevations, threshold)
gain_r = round(gain, 1)
loss_r = round(loss, 1)
# Patch detail JSON
detail["elevation_gain_m"] = gain_r
detail["elevation_loss_m"] = loss_r
json_path.write_text(json.dumps(detail, indent=2, ensure_ascii=False), encoding="utf-8")
# Patch index.json summary
index_path = user_dir / "index.json"
if index_path.exists():
index = json.loads(index_path.read_text(encoding="utf-8"))
for s in index.get("activities", []):
if s.get("id") == activity_id:
s["elevation_gain_m"] = gain_r
break
index_path.write_text(
json.dumps(index, indent=2, ensure_ascii=False), encoding="utf-8"
)
return {
"elevation_gain_m": gain_r,
"elevation_loss_m": loss_r,
"threshold_m": threshold,
"source": altitude_source,
}
+26 -2
View File
@@ -1267,8 +1267,8 @@ async def post_activity(
return JSONResponse({"ok": True}) return JSONResponse({"ok": True})
@app.post("/api/activity/{activity_id}/recalculate-elevation") @app.post("/api/activity/{activity_id}/recalculate-elevation/dem")
async def recalculate_elevation_endpoint( async def recalculate_elevation_dem_endpoint(
activity_id: str, activity_id: str,
bincio_session: Optional[str] = Cookie(default=None), bincio_session: Optional[str] = Cookie(default=None),
) -> JSONResponse: ) -> JSONResponse:
@@ -1296,6 +1296,30 @@ async def recalculate_elevation_endpoint(
raise HTTPException(422, str(e)) raise HTTPException(422, str(e))
@app.post("/api/activity/{activity_id}/recalculate-elevation/hysteresis")
async def recalculate_elevation_hysteresis_endpoint(
activity_id: str,
bincio_session: Optional[str] = Cookie(default=None),
) -> JSONResponse:
"""Recompute gain/loss from original recorded elevation using source-aware hysteresis."""
user = _require_user(bincio_session)
_check_id(activity_id)
dd = _get_data_dir() / user.handle
if not (dd / "activities" / f"{activity_id}.json").exists():
raise HTTPException(404, "Activity not found")
try:
from bincio.extract.dem import recalculate_elevation_hysteresis
from bincio.render.merge import merge_one
result = recalculate_elevation_hysteresis(dd, activity_id)
merge_one(dd, activity_id)
_trigger_rebuild(user.handle)
return JSONResponse(result)
except FileNotFoundError as e:
raise HTTPException(404, str(e))
except ValueError as e:
raise HTTPException(422, str(e))
@app.delete("/api/activity/{activity_id}", response_model=GenericResponse) @app.delete("/api/activity/{activity_id}", response_model=GenericResponse)
async def delete_activity( async def delete_activity(
activity_id: str, activity_id: str,
+22 -10
View File
@@ -24,8 +24,8 @@
let confirmDelete = false; let confirmDelete = false;
let deleting = false; let deleting = false;
// Elevation recalculation from DEM // Elevation recalculation
let recalculating = false; let recalculating: '' | 'dem' | 'hysteresis' = '';
let recalcStatus = ''; let recalcStatus = '';
let recalcOk = false; let recalcOk = false;
@@ -124,12 +124,12 @@
: [...hideStats, key]; : [...hideStats, key];
} }
async function recalculateElevation() { async function recalculateElevation(method: 'dem' | 'hysteresis') {
recalculating = true; recalculating = method;
recalcStatus = ''; recalcStatus = '';
recalcOk = false; recalcOk = false;
try { try {
const res = await fetch(`${api}/recalculate-elevation`, { method: 'POST' }); const res = await fetch(`${api}/recalculate-elevation/${method}`, { method: 'POST' });
const d = await res.json(); const d = await res.json();
if (!res.ok) throw new Error(d.detail ?? await res.text()); if (!res.ok) throw new Error(d.detail ?? await res.text());
recalcOk = true; recalcOk = true;
@@ -140,7 +140,7 @@
recalcStatus = e.message; recalcStatus = e.message;
recalcOk = false; recalcOk = false;
} finally { } finally {
recalculating = false; recalculating = '';
} }
} }
@@ -292,14 +292,26 @@
<!-- Elevation recalculation --> <!-- Elevation recalculation -->
<div class="mb-4"> <div class="mb-4">
<p class="text-xs text-zinc-500 mb-2">Elevation</p> <p class="text-xs text-zinc-500 mb-2">Elevation</p>
<div class="flex gap-2">
<button <button
type="button" type="button"
class="w-full flex items-center justify-center gap-2 px-3 py-2 rounded-lg border border-zinc-700 text-xs text-zinc-400 hover:border-zinc-500 hover:text-zinc-200 transition-colors disabled:opacity-40" class="flex-1 flex items-center justify-center gap-1 px-3 py-2 rounded-lg border border-zinc-700 text-xs text-zinc-400 hover:border-zinc-500 hover:text-zinc-200 transition-colors disabled:opacity-40"
disabled={recalculating} disabled={recalculating !== ''}
on:click={recalculateElevation} on:click={() => recalculateElevation('hysteresis')}
title="Recompute from the original recorded elevation using noise-filtering (fast, no network)"
> >
{recalculating ? 'Querying terrain data…' : ' Recalculate from terrain map (DEM)'} {recalculating === 'hysteresis' ? 'Computing…' : '📐 Recalculate (hysteresis)'}
</button> </button>
<button
type="button"
class="flex-1 flex items-center justify-center gap-1 px-3 py-2 rounded-lg border border-zinc-700 text-xs text-zinc-400 hover:border-zinc-500 hover:text-zinc-200 transition-colors disabled:opacity-40"
disabled={recalculating !== ''}
on:click={() => recalculateElevation('dem')}
title="Replace elevation with SRTM terrain data from the internet (slower, most accurate for GPS-only devices)"
>
{recalculating === 'dem' ? 'Querying terrain…' : '⛰ Recalculate (DEM)'}
</button>
</div>
{#if recalcStatus} {#if recalcStatus}
<p class="text-xs mt-1.5 text-center" class:text-green-400={recalcOk} class:text-red-400={!recalcOk}> <p class="text-xs mt-1.5 text-center" class:text-green-400={recalcOk} class:text-red-400={!recalcOk}>
{recalcStatus} {recalcStatus}