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
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
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
<pre lang="python">
#!/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 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
<kml xmlns="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.or
g/2005/Atom" xmlns:gx="http://www.google.com/kml/
ext/2.2" xmlns:kml="http://www.opengis.net/kml/2.2">
<document>
<name>result.kml</name>
<style id="lineStyle0">
<LineStyle>
<color>ff00aaff
<width>6
</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()