Journey Speed Plot

import psycopg2, datetime

# Change days_ago to 0 for today, 1 for yesterday, etc.
days_ago = 0

day_start = (datetime.datetime.now(datetime.timezone.utc).replace(hour=0, minute=0, second=0, microsecond=0)
             - datetime.timedelta(days=days_ago))
day_end = day_start + datetime.timedelta(days=1)
start = int(day_start.timestamp())
end = int(day_end.timestamp())

# 'who=john' maps to this device UUID via /etc/owntracks/config.json on gravlax
DEVICE_UUID = '53AECFB9-BA35-44B8-BB70-9A35D2B32500'

# Reads via the [postgres] stunnel tunnel (local :5432 -> gravlax:15432 -> owntracks Unix socket)
conn = psycopg2.connect(host='127.0.0.1', port=5432, dbname='owntracks', user='owntracks_ro')
cur = conn.cursor()
cur.execute('''
    SELECT extract(epoch from received_at)::bigint, tst, lat, lon, acc, batt, topic
    FROM locations
    WHERE topic LIKE %s AND tst >= %s AND tst <= %s
    ORDER BY tst ASC
''', (f'%{DEVICE_UUID}%', start, end))
rows = cur.fetchall()
cur.close()
conn.close()

locations = [
    {'received_at': r[0], 'tst': r[1], 'lat': r[2], 'lon': r[3], 'acc': r[4], 'batt': r[5], 'topic': r[6]}
    for r in rows
]
print(f'Fetched {len(locations)} fixes for {day_start.strftime("%Y-%m-%d")} (days_ago={days_ago}) via postgres')

import math, pandas as pd, numpy as np

BST = datetime.timezone(datetime.timedelta(hours=1))
R = 6371000

def east_north_m(lat1, lon1, lat2, lon2):
    p = math.pi / 180
    mean_lat = (lat1 + lat2) / 2 * p
    north = (lat2 - lat1) * p * R
    east  = (lon2 - lon1) * p * R * math.cos(mean_lat)
    return east, north

rows = []
prev = None
for loc in locations:
    t = datetime.datetime.fromtimestamp(loc['tst'], tz=BST)
    speed_mph = east_m = north_m = None
    if prev:
        dt = loc['tst'] - prev['tst']
        if 0 < dt < 300 and loc['acc'] <= 200:
            east, north = east_north_m(prev['lat'], prev['lon'], loc['lat'], loc['lon'])
            dist = math.hypot(east, north)
            speed_mph = (dist / dt) * 2.23694
            if speed_mph > 200:
                speed_mph = east_m = north_m = None
            else:
                east_m, north_m = east, north
    rows.append({'time': t, 'speed_mph': speed_mph, 'east_m': east_m, 'north_m': north_m})
    prev = loc

df = pd.DataFrame(rows).dropna(subset=['speed_mph'])
print(f'{len(df)} speed points computed')

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import seaborn as sns
import numpy as np, math, datetime
import psycopg2

RED  = np.array([0.85, 0.08, 0.08])
BLUE = np.array([0.08, 0.08, 0.85])
GREY = np.array([0.55, 0.55, 0.55])

def direction_colour(east, north):
    try:
        e, n = float(east), float(north)
    except (TypeError, ValueError):
        return tuple(float(x) for x in GREY)
    if not (math.isfinite(e) and math.isfinite(n)):
        return tuple(float(x) for x in GREY)
    mag = math.sqrt(e*e + n*n)
    if mag < 1e-3:
        return tuple(float(x) for x in GREY)
    t = e / mag
    rgb = GREY + t * (RED - GREY) if t >= 0 else GREY + (-t) * (BLUE - GREY)
    return tuple(float(x) for x in rgb)

def hav_m(lat1, lon1, lat2, lon2):
    R = 6371000; p = math.pi/180
    a = math.sin((lat2-lat1)*p/2)**2 + math.cos(lat1*p)*math.cos(lat2*p)*math.sin((lon2-lon1)*p/2)**2
    return 2*R*math.asin(math.sqrt(a))

# --- Match fixes against known stations (gazetteer table, via the
# [postgres] tunnel — see location-db/gazetteer). Query only the bounding
# box of today's route (with a margin) so we don't pull all 2606 stations
# every run, then do exact nearest-approach matching in Python.
lats = [loc['lat'] for loc in locations]
lons = [loc['lon'] for loc in locations]
margin = 0.05  # degrees, ~5 km

conn = psycopg2.connect(host='127.0.0.1', port=5432, dbname='owntracks', user='owntracks_ro')
cur = conn.cursor()
cur.execute('''
    SELECT name, lat, lon FROM gazetteer
    WHERE category = 'station'
      AND geom && ST_Expand(ST_MakeEnvelope(%s, %s, %s, %s, 4326), %s)
''', (min(lons), min(lats), max(lons), max(lats), margin))
stations = cur.fetchall()
cur.close()
conn.close()
print(f'{len(stations)} candidate stations in the route bounding box')

STATION_THRESHOLD_M = 500
BST = datetime.timezone(datetime.timedelta(hours=1))
visits = {}  # name -> (time_bst, dist_m)
for loc in locations:
    if loc.get('acc', 999) > 300:
        continue
    for name, slat, slon in stations:
        d = hav_m(loc['lat'], loc['lon'], slat, slon)
        if d < STATION_THRESHOLD_M and (name not in visits or d < visits[name][1]):
            t = datetime.datetime.fromtimestamp(loc['tst'], tz=BST)
            visits[name] = (t, d)

label_map = sorted(((t, name) for name, (t, d) in visits.items()), key=lambda x: x[0])
print(f'Found {len(label_map)} station(s) near the route:')
for t, name in label_map:
    print(f'  {t.strftime("%H:%M")}  {name}')

# --- Plot ---
colours = [direction_colour(r['east_m'], r['north_m']) for _, r in df.iterrows()]
times_num = mdates.date2num(df['time'].tolist())
widths = np.empty(len(times_num))
widths[:-1] = np.diff(times_num) * 0.85
widths[-1] = widths[-2] if len(widths) > 1 else 1/288

sns.set_theme(style='whitegrid')
fig, ax = plt.subplots(figsize=(14, 7.5))

ax.bar(times_num, df['speed_mph'].values, width=widths, color=colours, align='center', edgecolor='none')
ax.xaxis_date()
ax.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M'))
ax.xaxis.set_major_locator(mdates.AutoDateLocator())
ax.set_ylim(bottom=0)
ax.set_ylabel('Speed (mph)')
ax.set_xlabel('Time (BST)')
ax.set_title(f'Speed over time — {day_start.strftime("%Y-%m-%d")}')

from matplotlib.patches import Patch
ax.legend(handles=[
    Patch(color=tuple(RED),  label='East'),
    Patch(color=tuple(BLUE), label='West'),
    Patch(color=tuple(GREY), label='N / S'),
    Patch(color=direction_colour(1, 1),  label='NE/SE'),
    Patch(color=direction_colour(-1, 1), label='NW/SW'),
], loc='upper left', frameon=True, fontsize=8)

# Station labels above the plot: small font, rotated 45°, staggered across
# a few heights to cut down on overlap between neighbouring labels. x is a
# real date (data coords), y is axes-fraction (via get_xaxis_transform), so
# this works regardless of the speed scale.
trans = ax.get_xaxis_transform()
LEVELS = [1.05, 1.22, 1.39, 1.56]
for i, (t, name) in enumerate(label_map):
    tn = mdates.date2num(t)
    y = LEVELS[i % len(LEVELS)]
    ax.axvline(tn, color='#999', linewidth=0.6, ymin=0.0, ymax=y, clip_on=False, alpha=0.5)
    ax.annotate(
        name, xy=(tn, y), xycoords=trans,
        rotation=45, ha='left', va='bottom',
        fontsize=6, color='#333', clip_on=False, annotation_clip=False,
    )

# Headroom above the axes for the staggered labels.
fig.subplots_adjust(top=0.55, bottom=0.1)

show(fig)
plt.close(fig)