Cut, reverse and combine GPS tracks

When I plan a cycling trip I search for tracks in the target area and try to maximize the “fun factor”. Thus I often end up with a bunch of tracks from which I want to use parts. The difficult task is to cut, sometime reverse and combine the source tracks into a planned track. Most of the times this seems a too much waste of time since it involves a number of operations in qlandkarte (cutting and combining) and gpsbabel (reversing).

For this purpose here is a python script which takes as input a kml file containing all the tracks and some extra waypoints which define the cutting points then creates a combined track. Here is an example of how some original tracks look in google earth while planning a 2 day cycling trip in Parang mountains:

Original tracks

Original tracks

My goal is to specify that the final track will contain portions from track A from waypoint start to p1, from track B from waypoint p1 to p2 by reversing the original track and so on. The waypoints are manually specified so they don’t always sit on the track, the script will find the closest matching point on the tracks.

The format of the command is:

./parsePlan.py inputKml outputKml "track name/start waypoint/stop waypoint/[rev]" "track name/start waypoint/stop waypoint/[rev]" ...

Here are 2 example runs:

./parsePlan.py ~/parang.kml ~/plan-1.kml "Transalpina de la L. Oasa la Novaci/start/p1/rev" "Lacul Cilcescu - Pasul Urdele - Curmatura Oltetului - Lacul Petrimanu/p1/p2/rev" "parang-pejos/p2/p3/rev" "Transalpina de la L. Oasa la Novaci/p3/start/"

./parsePlan.py ~/parang.kml ~/plan-2.kml "Circuit Parang/start/p4/" "Path/p4/p5/" "Circuit Parang/p5/p6/" "Lacul Cilcescu - Pasul Urdele - Curmatura Oltetului - Lacul Petrimanu/p6/p1/rev" "Transalpina de la L. Oasa la Novaci/p1/start/"

and the resulted tracks

Cut, reverse, unify tracks

Cut, reverse, unify tracks

And here is the python script. It handles kml files and supports both the Placemark/LineString/coordinates and the Placemark/gx:Track/gx:coord formats. It’s quick and dirty but it worked for me and saved me more time than I required to write it

#!/usr/bin/env python
#$Revision: 1.2 $
# Copyright (C) 2011 Marilen Corciovei
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
 
import libxml2, sys, math, re, getopt, os
 
#http://www.johndcook.com/python_longitude_latitude.html
def distance_on_unit_sphere(lat1, long1, lat2, long2):
 
    # Convert latitude and longitude to 
    # spherical coordinates in radians.
    degrees_to_radians = math.pi/180.0
 
    # phi = 90 - latitude
    phi1 = (90.0 - lat1)*degrees_to_radians
    phi2 = (90.0 - lat2)*degrees_to_radians
 
    # theta = longitude
    theta1 = long1*degrees_to_radians
    theta2 = long2*degrees_to_radians
 
    # Compute spherical distance from spherical coordinates.
 
    # For two locations in spherical coordinates 
    # (1, theta, phi) and (1, theta, phi)
    # cosine( arc length ) = 
    #    sin phi sin phi' cos(theta-theta') + cos phi cos phi'
    # distance = rho * arc length
 
    cos = (math.sin(phi1)*math.sin(phi2)*math.cos(theta1 - theta2) + 
           math.cos(phi1)*math.cos(phi2))
    arc = math.acos( cos )
 
    # Remember to multiply arc by the radius of the earth 
    # in your favorite set of units to get length.
    return arc*6373
 
def distance(c1, c2):
    c1x = float(c1[0])
    c1y = float(c1[1])
    c2x = float(c2[0])
    c2y = float(c2[1])
    return distance_on_unit_sphere(c1x, c1y, c2x, c2y)
 
# removes default namespace
def clean_data(file):
    data = open(file, 'r').read();
    data = re.sub('xmlns="[^"]*"\s*', '', data)
    return data
 
def parse_cp(cp_file, day):
    global cps
    doc = libxml2.parseDoc(clean_data(cp_file))
    ctxt = doc.xpathNewContext()
    ctxt.xpathRegisterNs('kml', "http://www.opengis.net/kml/2.2")
    res = ctxt.xpathEval("//Placemark")
    for e in res:
        n = e.xpathEval('name')[0]
        name = n.content.strip()
        if name.startswith('cp') and name.endswith(day):
            coord = e.xpathEval('Point/coordinates')[0].content.split(',')
            cps.append([name, coord])
 
    doc.freeDoc()
    ctxt.xpathFreeContext()
 
    cps = sorted(cps, key = lambda cp : cp[0])
    print cps
 
def parse_track_file(track_file, precision, cp_stats):
    doc = libxml2.parseDoc(clean_data(track_file))
    ctxt = doc.xpathNewContext()
    res = ctxt.xpathEval("//Placemark/Point/coordinates")
    for e in res:
        pos = e.content.strip().split(',')
        for c in cp_stats:
            if not c[1]:
                dst = distance(c[0][1], pos)
                if dst < precision:
                    print 'Near cp', c[0][0], 'by', dst
                    c[1] = True
    doc.freeDoc()
    ctxt.xpathFreeContext()
 
def findWaypoint(ctxt, name):
    if name and len(name) > 0:
        res = ctxt.xpathEval('//Placemark[name="%s"]' % name)
        if res[0]:
            point = res[0].xpathEval('Point/coordinates')[0].content.split(',') 
            print point
            return point
    return None
 
def findClosestPoint(coordinates, point):
    if point:
        track_point = None
        i = 0
        for c in coordinates:
            coord = c.split(',')
            dist = distance(coord, point)
            if (not track_point) or (track_point and dist < track_point[0]):
                track_point = [dist, coord, i]
            i = i+1
        return track_point
    return [0, coordinates[0], 0]
 
def parse_file(ikml, opts):
    doc = libxml2.parseDoc(clean_data(ikml))
    ctxt = doc.xpathNewContext()
    ctxt.xpathRegisterNs("gx", "http://www.google.com/kml/ext/2.2")
 
    okml = """<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/
ext/2.2" xmlns:kml="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.or
g/2005/Atom">
<Document>
        <name>result.kml</name>
        <Style id="lineStyle0">
                <LineStyle>
                        <color>ff00aaff</color>
                        <width>6</width>
                </LineStyle>
        </Style>
<Placemark>
    <name>output</name>
    <styleUrl>#lineStyle0</styleUrl>
    <LineString>
    <tessellate>1</tessellate>
    <coordinates>"""
 
    for o in opts:
        (track,start,stop,direction) = o.split('/')
        #TODO verify length
        startPoint = findWaypoint(ctxt, start)
        stopPoint = findWaypoint(ctxt, stop)
 
        res = ctxt.xpathEval('//Placemark[name="%s"]' % track)
        if res[0]:
            coordinates = []
            #find LineString/coordinates
            coordElem = res[0].xpathEval('LineString/coordinates')
            if coordElem:
                coordinates = coordElem[0].content.strip().split(' ')
            else:
                #find gx:Track
                ctxt.setContextNode(res[0])
                gxTrack = ctxt.xpathEval('gx:Track/gx:coord')
                for g in gxTrack:
                    coordinates.append(g.content.replace(' ',','))
            if 'rev' == direction:
                coordinates.reverse()
            #print coordinates
            startTrack = findClosestPoint(coordinates, startPoint)
            stopTrack = findClosestPoint(coordinates, stopPoint)
            print startTrack, stopTrack
            okml = okml + ' '.join(coordinates[startTrack[2]:stopTrack[2]]) + ' '
    okml = okml + """</coordinates>
                        </LineString>
                </Placemark>
</Document>
</kml>"""
    return okml
    doc.freeDoc()
    ctxt.xpathFreeContext()
 
def main():
    ikml = sys.argv[1]
    okml = sys.argv[2]
    opts = sys.argv[3:]
    okmlData = parse_file(ikml, opts)
    open(okml,'w').write(okmlData)
 
if __name__=='__main__':
    main()

 

Leave a Reply

*