movement

Runnable sheet for fetching OwnTracks locations directly from the SQLite database on gravlax.

import datetime as _dt
import json
import shlex
import subprocess

GRAVLAX = 'john@gravlax.critchley.biz'
DB_PATH = '/var/lib/owntracks/locations.db'

def _remote_sqlite_json(sql, params=None):
    remote_script = r'''
import json
import sqlite3
import sys
payload = json.load(sys.stdin)
conn = sqlite3.connect(payload['db_path'])
conn.row_factory = sqlite3.Row
rows = conn.execute(payload['sql'], payload.get('params') or []).fetchall()
print(json.dumps([dict(row) for row in rows]))
'''
    payload = json.dumps({'db_path': DB_PATH, 'sql': sql, 'params': params or []})
    remote_cmd = 'sudo python3 -c ' + shlex.quote(remote_script)
    proc = subprocess.run(
        ['ssh', GRAVLAX, remote_cmd],
        input=payload,
        text=True,
        capture_output=True,
        timeout=30,
    )
    if proc.returncode != 0:
        raise RuntimeError(proc.stderr.strip() or f'ssh returned {proc.returncode}')
    return json.loads(proc.stdout)

def fetch_locations(limit=100, start=None, end=None):
    where = []
    params = []
    if start is not None:
        where.append('tst >= ?')
        params.append(int(start))
    if end is not None:
        where.append('tst <= ?')
        params.append(int(end))
    clause = (' where ' + ' and '.join(where)) if where else ''
    sql = f'''
        select received_at, tst, lat, lon, acc, batt, topic
        from locations
        {clause}
        order by tst desc
        limit ?
    '''
    params.append(int(limit))
    rows = _remote_sqlite_json(sql, params)
    rows.reverse()
    return rows

locations = fetch_locations(limit=100)
print(f'fetched {len(locations)} locations from {GRAVLAX}:{DB_PATH}')
if locations:
    first = _dt.datetime.fromtimestamp(locations[0]['tst'], _dt.UTC).isoformat()
    last = _dt.datetime.fromtimestamp(locations[-1]['tst'], _dt.UTC).isoformat()
    print('range:', first, 'to', last)
    print('latest:', locations[-1])

Summarise the fetched fixes and print the latest rows.

import datetime as _dt
import math

if 'locations' not in globals():
    raise RuntimeError('Run fetch_locations first')

def haversine_m(a, b):
    r = 6371000.0
    lat1, lon1 = math.radians(a['lat']), math.radians(a['lon'])
    lat2, lon2 = math.radians(b['lat']), math.radians(b['lon'])
    dlat = lat2 - lat1
    dlon = lon2 - lon1
    x = math.sin(dlat / 2) ** 2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon / 2) ** 2
    return 2 * r * math.asin(math.sqrt(x))

total_m = 0.0
moving_pairs = 0
for prev, cur in zip(locations, locations[1:]):
    d = haversine_m(prev, cur)
    total_m += d
    if d >= 5:
        moving_pairs += 1

print(f'points: {len(locations)}')
print(f'path distance over fetched points: {total_m / 1000:.2f} km')
print(f'pairs moving at least 5m: {moving_pairs}')
print('\nlast 10 fixes:')
for row in locations[-10:]:
    ts = _dt.datetime.fromtimestamp(row['tst'], _dt.UTC).strftime('%Y-%m-%d %H:%M:%S UTC')
    acc = row.get('acc')
    batt = row.get('batt')
    print(f"{ts}  {row['lat']:.6f}, {row['lon']:.6f}  acc={acc}m batt={batt}%")

Optional bounded query template. Edit start and end Unix timestamps, then run this cell.

# Example: fetch a bounded time range by Unix timestamps.
# Set start/end, then run this cell.
start = None  # e.g. 1783680000
end = None    # e.g. 1783683600
limit = 1000

range_locations = fetch_locations(limit=limit, start=start, end=end)
print(f'fetched {len(range_locations)} rows for range start={start!r} end={end!r}')

Seaborn/matplotlib plot of yesterday’s speed after combining nearby/noisy readings. Direction colour: east red, west blue, north/south green, diagonals blended.

import datetime as _dt
import math

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection
from matplotlib.colors import to_rgba
import pandas as pd
import seaborn as sns

if 'fetch_locations' not in globals():
    raise RuntimeError('Run fetch_locations first so the remote SQLite helper is defined')

MAX_PLAUSIBLE_KPH = 260.0
BIN_SECONDS = 5 * 60
MAX_ACCURACY_M = 100.0

# Default to yesterday's journey. UTC day boundaries are fine here because
# there was no travel around midnight.
target_day = _dt.date.today() - _dt.timedelta(days=1)
start_of_day_dt = _dt.datetime(target_day.year, target_day.month, target_day.day, tzinfo=_dt.UTC)
end_of_day_dt = start_of_day_dt + _dt.timedelta(days=1)
start_of_day = int(start_of_day_dt.timestamp())
end_of_day = int(end_of_day_dt.timestamp())
raw = fetch_locations(limit=10000, start=start_of_day, end=end_of_day)
raw = [r for r in raw if r.get('lat') is not None and r.get('lon') is not None]
raw = [r for r in raw if float(r.get('acc') or 9999) <= MAX_ACCURACY_M]
raw.sort(key=lambda r: int(r['tst']))

if len(raw) < 2:
    raise RuntimeError(f'Need at least two usable fixes for today; got {len(raw)}')

def haversine_m(a, b):
    r = 6371000.0
    lat1, lon1 = math.radians(float(a['lat'])), math.radians(float(a['lon']))
    lat2, lon2 = math.radians(float(b['lat'])), math.radians(float(b['lon']))
    dlat = lat2 - lat1
    dlon = lon2 - lon1
    x = math.sin(dlat / 2) ** 2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon / 2) ** 2
    return 2 * r * math.asin(math.sqrt(x))

def weighted_point(rows):
    weights = [1.0 / max(float(r.get('acc') or 25.0), 1.0) for r in rows]
    total = sum(weights)
    return {
        'tst': sum(float(r['tst']) * w for r, w in zip(rows, weights)) / total,
        'received_at': sum(float(r.get('received_at') or r['tst']) * w for r, w in zip(rows, weights)) / total,
        'lat': sum(float(r['lat']) * w for r, w in zip(rows, weights)) / total,
        'lon': sum(float(r['lon']) * w for r, w in zip(rows, weights)) / total,
        'acc': min(float(r.get('acc') or 9999) for r in rows),
        'n': sum(int(r.get('n', 1)) for r in rows),
    }

def bin_readings(rows):
    bins = {}
    for row in rows:
        bucket = int((int(row['tst']) - start_of_day) // BIN_SECONDS)
        bins.setdefault(bucket, []).append(row)
    return [weighted_point(bins[k]) for k in sorted(bins)]

def speed_kph(a, b):
    dt = float(b['tst']) - float(a['tst'])
    if dt <= 0:
        return 0.0
    return haversine_m(a, b) / dt * 3.6

points = bin_readings(raw)
original_bins = len(points)
merge_count = 0
for _ in range(8):
    changed = False
    rebuilt = [points[0]]
    for point in points[1:]:
        prev = rebuilt[-1]
        if speed_kph(prev, point) > MAX_PLAUSIBLE_KPH:
            merged = weighted_point([prev, point])
            merged['n'] = int(prev.get('n', 1)) + int(point.get('n', 1))
            rebuilt[-1] = merged
            merge_count += 1
            changed = True
        else:
            rebuilt.append(point)
    points = rebuilt
    if not changed:
        break

if len(points) < 2:
    p = weighted_point(raw)
    points = [dict(p, tst=start_of_day), dict(p, tst=int(now.timestamp()))]

def east_north_m(a, b):
    mean_lat = math.radians((float(a['lat']) + float(b['lat'])) / 2.0)
    north = (float(b['lat']) - float(a['lat'])) * 111320.0
    east = (float(b['lon']) - float(a['lon'])) * 111320.0 * math.cos(mean_lat)
    return east, north

def direction_colour(east, north):
    mag = math.hypot(east, north)
    if mag <= 0:
        return (0.55, 0.55, 0.55, 1.0)
    red = max(east, 0.0) / mag
    blue = max(-east, 0.0) / mag
    green = abs(north) / mag
    peak = max(red, green, blue, 1e-9)
    # Keep colours visible while preserving the direction mix.
    return (0.18 + 0.82 * red / peak, 0.18 + 0.82 * green / peak, 0.18 + 0.82 * blue / peak, 1.0)

rows = []
segments = []
colours = []
max_speed = 0.0
for prev, cur in zip(points, points[1:]):
    dt = float(cur['tst']) - float(prev['tst'])
    if dt <= 0:
        continue
    east, north = east_north_m(prev, cur)
    dist = math.hypot(east, north)
    spd = min(dist / dt * 3.6, MAX_PLAUSIBLE_KPH)
    max_speed = max(max_speed, spd)
    mid_t = (float(prev['tst']) + float(cur['tst'])) / 2.0
    t = _dt.datetime.fromtimestamp(mid_t, _dt.UTC)
    hour = (mid_t - start_of_day) / 3600.0
    rows.append({'time': t, 'hour': hour, 'speed_kph': spd, 'east_m': east, 'north_m': north})

if not rows:
    rows = [{'time': start_of_day_dt, 'hour': 0.0, 'speed_kph': 0.0, 'east_m': 0.0, 'north_m': 0.0}]

df = pd.DataFrame(rows)
sns.set_theme(style='whitegrid')
fig, ax = plt.subplots(figsize=(11, 4.8))
# Seaborn provides the axes/style and light context; LineCollection gives per-segment colours.
sns.lineplot(data=df, x='hour', y='speed_kph', ax=ax, color='0.82', linewidth=1.0, estimator=None)
if len(df) >= 2:
    xy = df[['hour', 'speed_kph']].to_numpy()
    segments = [[xy[i], xy[i + 1]] for i in range(len(xy) - 1)]
    colours = [direction_colour(df.iloc[i + 1]['east_m'], df.iloc[i + 1]['north_m']) for i in range(len(df) - 1)]
    lc = LineCollection(segments, colors=colours, linewidths=3.0)
    ax.add_collection(lc)
else:
    ax.scatter(df['hour'], df['speed_kph'], color='0.55')
ax.scatter(df['hour'], df['speed_kph'], s=10, color='black', alpha=0.35, zorder=3)
ax.set_title(f'Speed through {target_day.isoformat()}, coloured by direction')
ax.set_xlabel('hours since midnight UTC')
ax.set_ylabel('speed (km/h)')
ax.set_xlim(0, max(24, float(df['hour'].max()) + 0.25))
threshold_kph = 125 * 1.609344
ax.axhline(threshold_kph, color='black', linestyle='--', linewidth=1.5, alpha=0.75)
ax.text(0.15, threshold_kph + 3, '125 mph / 201 km/h', color='black', fontsize=9)
ax.set_ylim(0, max(threshold_kph * 1.1, 1.0, float(df['speed_kph'].max()) * 1.2))
legend_items = [
    ('east', '#d62728'),
    ('west', '#1f77b4'),
    ('north/south', '#2ca02c'),
    ('diagonal = blended colour', '#666666'),
]
for i, (label, colour) in enumerate(legend_items):
    ax.plot([], [], color=colour, linewidth=3, label=label)
ax.legend(loc='upper right', frameon=True)
fig.tight_layout()

print(f'plot day: {target_day.isoformat()} UTC')
print(f'raw fixes: {len(raw)}')
print(f'5-minute bins after accuracy filtering: {original_bins}')
print(f'remediated points: {len(points)}')
print(f'insane-speed merges applied: {merge_count}')
print(f'max remaining segment speed: {max_speed:.1f} km/h')
show(fig)
plt.close(fig)