split-gpx

#!/usr/bin/python

import xml.etree.ElementTree as et
import argparse
from datetime import datetime, timedelta
from pathlib import Path
import re
import math

# lenghts/distances are in km


name_timeformat = '%Y-%m-%d %H:%M'

time_units = dict(
    ms='milliseconds',
    s='seconds',
    m='minutes',
    h='hours',
    d='days',
    w='weeks',
)
default_time_unit = 'm'


duration_re = re.compile(r'^(\d+)\s*(' + '|'.join(list(time_units.keys()) + list(time_units.values())) + r')?$')


def duration(s):
    s = s.strip()

    m = duration_re.match(s)
    if not m:
        raise ValueError(f'malformed duration: {s}. Use {duration_re}.')

    num = m.group(1)
    unit = m.group(2)

    num = float(num)
    if unit in time_units:
        unit = time_units[unit]
    kwargs = {unit: num}
    return timedelta(**kwargs)



def vel(d, dt):
    s = dt.seconds
    if s == 0:
        assert d < 1
        return 0
    return d / s * 60 * 60


def haversine(lon1, lat1, lon2, lat2):
    """
    Calculate the great circle distance between two points
    on the earth (specified in decimal degrees).
    Source: http://gis.stackexchange.com/a/56589/15183
    """
    # convert decimal degrees to radians
    lon1, lat1, lon2, lat2 = map(math.radians, [lon1, lat1, lon2, lat2])
    # haversine formula
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    a = (
        math.sin(dlat / 2) ** 2
        + math.cos(lat1) * math.cos(lat2) * math.sin(dlon / 2) ** 2
    )
    c = 2 * math.asin(math.sqrt(a))
    km = 6367 * c
    return km


if __name__ == '__main__':
    parser = argparse.ArgumentParser(
        description="Detect pauses in a monolithic gpx track and split it into segments."
    )
    parser.add_argument('ingpx', help="Input gpx")
    parser.add_argument('outgpx', help="Output gpx")
    parser.add_argument('--pauses', help="Include pause tracks", action='store_true')
    parser.add_argument(
        '--excluded',
        help="Include tracks that are too short, too near, too slow",
        action='store_true',
    )
    parser.add_argument(
        '--split', help="Each track into its own file", action='store_true'
    )

    parser.add_argument('--pause-duration-min', help='Pauses need to be at least this long.', default=timedelta(hours=1.5), type=duration)
    parser.add_argument('--pause-velocity-max', help='While GPS signal is lost, recognize pause only when speed between points is below this value. (km/h)', default=2, type=float)
    points_moved_dist_default = 0.15
    parser.add_argument('--pause-distance-max', help='Max distance between consecutive points of a pause. Tune this according to your GPS precision.', default=points_moved_dist_default, type=float)
    # distance between points to check
    # (tune this according to your gps precision)
    # TODO better docu
    parser.add_argument('--route-distance-max', help='Maximum allowed distance for any point to check.', default=points_moved_dist_default + 0.5, type=float)
    parser.add_argument('--route-max-gap-duration', help='Time gap of this length forces start of a new track', default=timedelta(hours=4), type=duration)
    parser.add_argument('--route-length-min', help='Exclude tracks shorter than this value. (km)', default=0.4, type=float)
    parser.add_argument('--route-duration-min', help='Exclude tracks shorter than this value. (duration)', default=timedelta(minutes=5), type=duration)
    parser.add_argument('--route-velocity-min', help='Exclude tracks slower than this value. (km/h)', default=0.4, type=float)


    args = parser.parse_args()

    points_max_dist = args.pause_distance_max
    points_moved_dist = args.route_distance_max
    pause_min_timedelta = args.pause_duration_min
    pause_max_velocity = args.pause_velocity_max
    route_max_timedelta = args.route_max_gap_duration
    route_min_dist = args.route_length_min
    route_min_duration = args.route_duration_min
    route_min_velocity = args.route_velocity_min

    ns = {'gpx': 'http://www.topografix.com/GPX/1/1'}

    early = datetime(year=1970, day=1, month=1)

    tree = et.parse(args.ingpx)
    et.register_namespace('', ns['gpx'])
    gpx = tree.getroot()
    tracks = []
    start_end = []
    distances = []
    for trk in gpx.findall('gpx:trk', ns):
        for seg in trk.findall('gpx:trkseg', ns):
            start = 'beginning'
            route = []
            route_start = None
            last = None
            pts = seg.findall('gpx:trkpt', ns)
            for ipt, pt in enumerate(pts):
                t = pt.find('gpx:time', ns).text
                #'2016-08-14T12:36:24.312Z'
                # fixme: parse float seconds
                assert t[-1] == 'Z'
                assert t[-5] == '.'
                t = datetime.strptime(t[:-5], '%Y-%m-%dT%H:%M:%S')

                pos = [float(x) for x in [pt.attrib['lon'], pt.attrib['lat']]]

                route_d = haversine(*route[-1][0], *pos) if route else 999999

                # we have enough distance from last point
                moved = route_d > points_moved_dist
                pause_start = None
                if route and moved:
                    dt = t - route[-1][1]
                    assert dt.seconds >= 0
                    velocity = vel(route_d, dt)
                    # last point was long ago, and we didnt travel during this time
                    if (
                        dt > pause_min_timedelta and velocity < pause_max_velocity
                    ) or dt > route_max_timedelta:
                        pause_start = route[-1]

                if pause_start:
                    # print('pause of %s [from %s to %s] (%s km, %s km/h) %s points' % (
                    #     dt
                    #     , pause_start[1], t
                    #     , round(pause_start[3], 2)
                    #     , round(velocity, 1)
                    #     , ipt - pause_start[2]
                    #     ))

                    def add_track(comment, start, end):
                        # include last and pause point
                        ps = pts[start[2] : end[2] + 1]
                        tracks.append(ps)
                        start_end.append((comment, start[1], end[1]))
                        dist = end[3] - start[3]
                        distances.append(dist)

                        dt = end[1] - start[1]
                        assert dt.seconds >= 0
                        velocity = vel(dist, dt)

                        reason = '%s\t(%s km,\t%s,\t%s km/h)\twith %s points' % (
                            comment,
                            round(dist, 2),
                            dt,
                            round(velocity, 2),
                            len(ps),
                        )

                        print(reason)

                    if True:  # maxdist > remove_diff_dist:
                        comment = ''
                        d = pause_start[3] - route_start[3]
                        dt = pause_start[1] - route_start[1]
                        if d < route_min_dist:
                            comment = 'near'
                        elif dt < route_min_duration:
                            comment = 'short'
                        elif vel(d, dt) < route_min_velocity:
                            comment = 'slow'
                        if not comment or args.excluded:
                            add_track(comment, route_start, pause_start)

                        if last[2] > pause_start[2] and args.pauses:
                            add_track('pause', pause_start, last)

                    route_start = last
                    route = [last]
                else:
                    # remove moved route points
                    for io, o in enumerate(route):
                        d = haversine(*o[0], *pos)
                        dt = t - o[1]
                        if d < points_max_dist or dt < pause_min_timedelta:
                            # remove moved away points unless they're recent
                            route = route[io:]
                            break

                dist = 0
                if last:
                    d = haversine(*last[0], *pos)
                    dist = last[3] + d

                last = pos, t, ipt, dist
                if moved:
                    route.append(last)
                if not route_start:
                    route_start = last

            if route_start[2] < len(pts):
                add_track('', route_start, last)
        gpx.remove(trk)
    print(len(tracks))
    assert len(tracks) == len(start_end)
    for i, (pts, time, dist) in enumerate(zip(tracks, start_end, distances)):
        trk = et.SubElement(gpx, 'trk')
        name = et.SubElement(trk, 'name')
        name.text = (
            time[0]
            + ' - '.join([t.strftime(name_timeformat) for t in time[1:]])
            + ' '
            + str(round(dist, 2))
            + ' km'
        )
        seg = et.SubElement(trk, 'trkseg')
        for pt in pts:
            seg.append(pt)
        if args.split:
            outf = Path(args.outgpx)
            outf = Path(outf.stem + str(i) + str(time[1])).with_suffix(outf.suffix)
            print(outf)
            tree.write(str(outf), encoding='utf-8', xml_declaration=True)
            gpx.remove(trk)

    if not args.split:
        tree.write(args.outgpx, encoding='utf-8', xml_declaration=True)