count.py

#!/usr/bin/env python3
"""
count.py - read GPSLogger .gpx files, work out which country each day was spent in,
drop physically impossible "teleport" points, and write data/days.json.

Reads:  $DACO_GPS_DIR  (default /home/ubuntu/gpslogger)
Writes: <project>/data/days.json
"""

import os, glob, json, math
from collections import defaultdict, Counter
from datetime import datetime, timezone
import xml.etree.ElementTree as ET
import reverse_geocoder as rg

GPS_DIR = os.environ.get("DACO_GPS_DIR", "/home/ubuntu/gpslogger")
PROJECT = os.path.dirname(os.path.abspath(__file__))
DATA    = os.path.join(PROJECT, "data")

MAX_KMH = 1100          # faster than any airliner -> a "teleport", almost certainly a bad fix / spoof

SCHENGEN = {"AT","BE","CZ","DK","EE","FI","FR","DE","GR","HU","IS","IT","LV","LI","LT",
            "LU","MT","NL","NO","PL","PT","SK","SI","ES","SE","CH","HR","RO","BG"}
NAMES = {"CY":"Cyprus","IL":"Israel","UA":"Ukraine","ES":"Spain","GR":"Greece","IT":"Italy",
         "FR":"France","DE":"Germany","GB":"United Kingdom","NL":"Netherlands","PT":"Portugal",
         "CH":"Switzerland","AT":"Austria","TR":"Turkey","AE":"UAE","US":"United States",
         "PL":"Poland","RO":"Romania","BG":"Bulgaria","HR":"Croatia"}


def localname(tag): return tag.rsplit("}", 1)[-1]

def parse_time(s):
    if not s: return None
    try: return datetime.fromisoformat(s.replace("Z", "+00:00"))
    except ValueError: return None

def read_points():
    pts = []
    for f in sorted(glob.glob(os.path.join(GPS_DIR, "*.gpx"))):
        try:
            root = ET.fromstring(open(f, encoding="utf-8").read())
        except ET.ParseError:
            continue
        for el in root.iter():
            if localname(el.tag) != "trkpt": continue
            lat, lon = el.get("lat"), el.get("lon")
            if lat is None or lon is None: continue
            t = None
            for c in el:
                if localname(c.tag) == "time": t = parse_time((c.text or "").strip())
            if t is None: continue
            pts.append({"lat": float(lat), "lon": float(lon), "t": t, "file": os.path.basename(f)})
    pts.sort(key=lambda p: p["t"])
    return pts

def haversine(a, b):
    R = 6371.0
    p1, p2 = math.radians(a[0]), math.radians(b[0])
    dphi = math.radians(b[0]-a[0]); dl = math.radians(b[1]-a[1])
    h = math.sin(dphi/2)**2 + math.cos(p1)*math.cos(p2)*math.sin(dl/2)**2
    return 2*R*math.asin(min(1, math.sqrt(h)))

def flag_teleports(pts):
    """mark a point as anomaly if both the jump in and the jump out exceed MAX_KMH (isolated spike)."""
    n = len(pts)
    for p in pts: p["anomaly"] = False
    for i in range(n):
        def speed(j, k):
            dt_h = abs((pts[k]["t"] - pts[j]["t"]).total_seconds()) / 3600.0
            if dt_h == 0: return float("inf")
            return haversine((pts[j]["lat"], pts[j]["lon"]), (pts[k]["lat"], pts[k]["lon"])) / dt_h
        sp_in  = speed(i-1, i) if i > 0 else 0
        sp_out = speed(i, i+1) if i < n-1 else 0
        if sp_in > MAX_KMH and sp_out > MAX_KMH:
            pts[i]["anomaly"] = True
    return pts

def main():
    os.makedirs(DATA, exist_ok=True)
    pts = read_points()
    if not pts:
        json.dump({"days": [], "generated": datetime.now(timezone.utc).isoformat(), "gps_dir": GPS_DIR,
                   "point_count": 0}, open(os.path.join(DATA, "days.json"), "w"))
        print("no points found in", GPS_DIR); return

    geo = rg.search([(p["lat"], p["lon"]) for p in pts], mode=1)
    for p, r in zip(pts, geo):
        p["cc"] = r["cc"]
        p["admin1"] = r.get("admin1", "") or ""   # region / US state / province
    flag_teleports(pts)

    by_day = defaultdict(list)
    for p in pts:
        by_day[p["t"].date().isoformat()].append(p)

    days = []
    for d in sorted(by_day):
        dpts = by_day[d]
        clean = [p for p in dpts if not p["anomaly"]] or dpts
        counts = Counter(p["cc"] for p in clean)
        dom = counts.most_common(1)[0][0]
        adm = Counter(p["admin1"] for p in clean if p["cc"] == dom and p["admin1"])
        admin1 = adm.most_common(1)[0][0] if adm else ""
        seen = sorted(set(p["cc"] for p in dpts))
        flagged = any(p["anomaly"] for p in dpts) or len(set(p["cc"] for p in clean)) > 1
        days.append({"date": d, "cc": dom, "country": NAMES.get(dom, dom), "admin1": admin1,
                     "schengen": dom in SCHENGEN, "flagged": flagged,
                     "countries_seen": seen, "points": len(dpts)})

    out = {"generated": datetime.now(timezone.utc).isoformat(), "gps_dir": GPS_DIR,
           "point_count": len(pts), "days": days}
    json.dump(out, open(os.path.join(DATA, "days.json"), "w"), ensure_ascii=False, indent=1)
    print(f"points: {len(pts)} | days: {len(days)} | flagged days: {sum(d['flagged'] for d in days)}")
    print("wrote", os.path.join(DATA, "days.json"))

if __name__ == "__main__":
    main()