#!/usr/bin/env python2
"""
Convert GPX track logs into KML with some processing.

By default, each track in the GPX file generates a new track in the KML. In
addition, tracks are split at waypoints (e.g. if the button on the logger has
been pressed).

By defining a "split file", the tracks can be customized. If the split file
does not exist when the program is run, a new file with the split times derived
with the standard rules is generated. On subsequent runs, only the settings in
the split file are relevant and the standard splitting rules are not applied.

The split file is a text file where each line contains a timestamp in ISO
format, UTC time. Lines starting with '#' are comments. By adding timestamps,
tracks can be split; by removing timestamps, tracks can be joined. The
timestamps must be ordered chronologically. Track points before the first
timestamp will be discarded. If the timestamp is followed by whitespace and a
hyphen ('-'), the track segment will be discarded as well. If the timestamp is
followed by a caret ('^'), the track segment will not be turned into a track of
its own, but appended to the previous track that has not been discarded.

Following the timestamp, a line may contain a name of the track segment that
starts at that time. To specify the character set used for the titles, write a
line like '\encoding iso8859-1' or '\encoding utf-8' at the beginning of the
file.

Tracks can be sorted into KML subfolders by specifying the target folder name
on a '\folder <name>' line. All following splits will be put into a folder of
that name. A '\folder' line without a name disables folder sorting.

The split file can also be used to customize the graphical appearance of the
generated tracks in the KML file. To this extent, line styles can be defined
with the special command '\defstyle .name <xml>', where 'name' is the name of
the style and the remainder of the line contains a KML line style specification.
An example is included automatically when a new split file is generated.
To make use of the styles, the name (including the leading dot) has to be
written on a line, directly after the timestamp, before the title.

A 'full-featured' example of a line in a split file might be:
    2012-12-21T13:37:42Z .myStyle This is an example track

It is also possible to include other split files from a central split file by
usging the '\include <filename>' command. Some limits apply: Files can only be
included once, and while the defined styles are shared globally, each file has
its own independent encoding and folder settings.

In addition to this advanced track splitting, the tool supports automatic
simplification of tracks using a defined threshold, and generating placemarks
for JPEG files. If the latter feature is used, please not that the timestamps
of the photos will *not* be read from Exif or other metadata; instead, the
program relies on the file modification time to be set properly. The images
will be referenced directly in the KML file; no embedding or transcoding of
any kind will be done.
"""
__version__ = "1.1.0"
__author__ = "Martin J. Fiedler <martin.fiedler@gmx.net>"
__copyright__ = """
This software is published under the terms of KeyJ's Research License,
version 0.2. Usage of this software is subject to the following conditions:
0. There's no warranty whatsoever. The author(s) of this software can not
   be held liable for any damages that occur when using this software.
1. This software may be used freely for both non-commercial and commercial
   purposes.
2. This software may be redistributed freely as long as no fees are charged
   for the distribution and this license information is included.
3. This software may be modified freely except for this license information,
   which must not be changed in any way.
4. If anything other than configuration, indentation or comments have been
   altered in the code, the original author(s) must receive a copy of the
   modified code.
"""
import sys, os, stat, re, time, calendar, optparse, xml.parsers.expat
from math import *

DEFAULT_STYLE_NAME = 'defaultTrackStyle'
PHOTO_STYLE_NAME = 'photoStyle'
DEFAULT_STYLES = {
    DEFAULT_STYLE_NAME: '<LineStyle><color>ffff0000</color><width>3</width></LineStyle><LabelStyle><color>00ffffff</color><scale>0</scale></LabelStyle><IconStyle><color>00ffffff</color><scale>0</scale><Icon></Icon></IconStyle>',
    PHOTO_STYLE_NAME: '<LabelStyle><scale>0</scale></LabelStyle><IconStyle><scale>0.7</scale><Icon><href>http://maps.google.com/mapfiles/kml/shapes/camera.png</href></Icon></IconStyle>',
}
DEFAULT_PHOTO_DIAG = 1000

###############################################################################

def IsValidTime(t):
    return (t > 86400)

def Time2Human(t):
    return time.strftime("%Y-%m-%d %H:%M:%S", time.localtime(t))

def Time2ISO(t):
    return time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime(t))

re_isotime = re.compile(r'(?P<dy>\d{4})-(?P<dm>\d\d?)-(?P<dd>\d\d?)[tT](?P<th>\d\d?):(?P<tm>\d\d?):(?P<ts>\d\d?)(\.\d+)?([zZ]|(?P<tzh>[+-]\d\d)(?P<tzm>\d\d)?)?$')
def ISO2Time(t):
    m = re_isotime.match(t)
    if not m:
        raise ValueError("invalid ISO 8601 time stamp")
    t = calendar.timegm(map(int, [m.group('dy'), m.group('dm'), m.group('dd'), m.group('th'), m.group('tm'), m.group('ts')]) + [0, 0, -1])
    try:
        t -= int(m.group('tzh')) * 3600 + int(m.group('tzm')) * 60
    except (TypeError, ValueError):
        pass
    return t

def X(s):
    return s.encode('utf-8', 'replace').replace('&', '&amp;').replace('<', '&lt;').replace('>', '&gt;')

def sqr(x): return x*x

def calculate_dop(a, b):
    try:
        x = b*b
        if b < 0.0:
            x = -x
        x = sqrt(a*a + x)
    except ValueError:
        x = 1.0
    return max(0.1, min(99.99, x))

def hermite(a, b, c, d, t):
    return (((3.0 * (b - c) - a + b) * t + 2.0 * a - 5.0 * b + 4.0 * c - d) * t - a + c) * t + 2.0 * b

###############################################################################

DEFAULT_VINCENTY_PRECISION = 1.0E-12

def GeoDistance(lat1, lon1, lat2, lon2, precision=DEFAULT_VINCENTY_PRECISION):
    """implementation of the Vincenty formula (C) 2002-2012 Chris Veness
    (http://www.movable-type.co.uk/scripts/latlong-vincenty.html)
    ported from JavaScript to Python"""
    a = 6378137
    b = 6356752.314245
    f = 1/298.257223563
    L = radians(lon2-lon1)
    U1 = atan((1-f) * tan(radians(lat1)))
    U2 = atan((1-f) * tan(radians(lat2)))
    sinU1 = sin(U1)
    cosU1 = cos(U1)
    sinU2 = sin(U2)
    cosU2 = cos(U2)
    la = L
    iterLimit = 100
    while True:
        sinla = sin(la)
        cosla = cos(la)
        sinSigma = sqrt((cosU2*sinla) * (cosU2*sinla) + (cosU1*sinU2-sinU1*cosU2*cosla) * (cosU1*sinU2-sinU1*cosU2*cosla))
        if not sinSigma:
            return 0.0
        cosSigma = sinU1*sinU2 + cosU1*cosU2*cosla
        sigma = atan2(sinSigma, cosSigma)
        sinAlpha = cosU1 * cosU2 * sinla / sinSigma
        cosSqAlpha = 1 - sinAlpha*sinAlpha
        if cosSqAlpha:
            cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha
        else:
            cos2SigmaM = 0.0
        C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha))
        laP = la
        la = L + (1-C) * f * sinAlpha * (sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)))
        iterLimit -= 1
        if (abs(la-laP) < precision) or not(iterLimit):
            break
    if not iterLimit:
        return 0.0
    uSq = cosSqAlpha * (a*a - b*b) / (b*b);
    A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)))
    B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)))
    deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)))
    return b*A*(sigma-deltaSigma)

def ApproxGeoDistance(lat1, lon1, lat2, lon2, precision=DEFAULT_VINCENTY_PRECISION, tolerance=0.01):
    """approximation of the distance between two points;
    only valid for very small differences, doesn't work near
    the poles and near the date line; falls back to Vincenty
    automatically if the values are problematic"""
    dlat = abs(lat1 - lat2)
    dlon = abs(lon1 - lon2)
    if (dlat > tolerance) or (dlon > tolerance):
        return GeoDistance(lat1, lon1, lat2, lon2, precision)
    cos_lat = cos(radians(0.5 * (lat1 + lat2)))
    dlat *= 110946.25761306843 + cos_lat * 373.23318020513398
    dlon *= 111319.49079327358 * cos_lat
    return sqrt(dlat * dlat + dlon * dlon)

###############################################################################

class TrackPoint:
    def __init__(self, lat, lon, alt=0, t=0, track_start=False):
        self.lat = float(lat)
        self.lon = float(lon)
        self.alt = float(alt)
        self.t = float(t)
        self.hdop = None
        self.vdop = None
        self.pdop = None
        self.track_start = track_start

    def has_time(self):
        return IsValidTime(self.t)

    def finish(self):
        if not self.vdop:
            if self.hdop and self.pdop:
                self.vdop = calculate_dop(self.pdop, -self.hdop)
            else:
                self.vdop = 1.0
        if not self.hdop:
            if self.vdop and self.pdop:
                self.hdop = calculate_dop(self.pdop, -self.vdop)
            else:
                self.hdop = 1.0
        if not self.pdop:
            self.pdop = calculate_dop(self.hdop, self.vdop)

    def __str__(self):
        return "[N %s E %s, %s m, %s%s, PDOP %.2f]" % (self.lat, self.lon, self.alt, Time2Human(self.t), (" *" if self.track_start else ""), self.pdop)


class Track:
    def __init__(self):
        self.data = []
        self.start_time = None
        self.style = DEFAULT_STYLE_NAME
        self.name = None
        self.folder = None

    def append(self, point):
        self.data.append(point)
        if not(self.start_time) and point.has_time():
            self.start_time = point.t

    def __getitem__(self, i):
        return self.data[i]

    def __len__(self):
        return len(self.data)

    def count_subtracks(self):
        return 1 + len([None for p in self.data[1:] if p.track_start])

    def get_name(self):
        if self.name:
            return unicode(self.name)
        if self.start_time:
            return unicode(Time2Human(self.start_time))
        return u'???'

    def sort(self):
        self.data.sort(key=lambda p: p.t)

    def simplify_once(self, tolerance=5.0, hdop_influence=0.5, max_hdop_factor=5.0):
        points_removed = 0
        i = 0
        while (i + 2) < len(self.data):
            i += 1
            l = self.data[i-1]
            p = self.data[i]
            r = self.data[i+1]
            if not(l.has_time()) or not(p.has_time()) or not(r.has_time()) \
            or p.track_start or r.track_start:
                continue
            a = float(p.t - l.t) / float(r.t - l.t)
            clat = l.lat + a * (r.lat - l.lat)
            clon = l.lon + a * (r.lon - l.lon)
            d = ApproxGeoDistance(clat, clon, p.lat, p.lon)
            # compute "fresnel zone" radius
            maxd = tolerance * sqrt(4.0 * (a - a*a))
            if hdop_influence and p.hdop:
                maxd *= min(pow(p.hdop, hdop_influence), max_hdop_factor)
            if d < maxd:
                points_removed += 1
                del self.data[i]
                i -= 1
        return points_removed

    def simplify(self, tolerance=5.0, hdop_influence=0.5, max_hdop_factor=5.0):
        total_points_removed = 0
        res = 1
        while res:
            res = self.simplify_once(tolerance, hdop_influence, max_hdop_factor)
            total_points_removed += res
        return total_points_removed

    def point_time(self, index):
        """determine the timestamp of a point with a specific index;
        if the point doesn't have a timestamp, it will be approximated;
        returns None if no timestamp is available at all"""
        if (index < 0) or (index >= len(self.data)):
            return None
        t = self.data[index].t
        if IsValidTime(t):
            return t

        # no direct hit -> search left neighbor
        il = index
        while True:
            il -= 1
            if il < 0:
                return None  # no left neighbor
            tl = self.data[il].t
            if IsValidTime(tl):
                break

        # search right neighbor
        ir = index
        while True:
            ir += 1
            if ir >= len(self.data):
                return None  # no right neighbor
            tl = self.data[il].t
            if IsValidTime(tl):
                break

        # interpolate
        return (tl * (ir - index) + tr * (index - il)) / float(ir - il)

    def point_index_at(self, t):
        """determine the index of the point just before a specific timestamp;
        returns None if no point can be found;
        point list must be sorted for this to work"""
        t0 = self.point_time(0)
        t1 = self.point_time(len(self.data) - 1)
        if (t0 and (t < t0)) or (t1 and (t > t1)):
            return None
        a = 0
        b = len(self.data)
        while (a + 1) < b:
            c = (a + b) >> 1
            tc = self.point_time(c)
            if (abs(tc - t) <= 0.001):
                return c
            if t < tc:
                b = c
            else:
                a = c
        return a

    def point_at(self, t, interpolate=1):
        """determine the point's position at a specific timestamp as tuple(lat, lon);
        returns None if no point can be found;
        interpolate=0 => no interpolation, return nearest trackpoint;
        interpolate=1 => linear interpolation;
        interpolate=2 => cubic interpolation"""
        idx = self.point_index_at(t)
        if idx is None:
            return None
        if idx == (len(self.data) - 1):
            return (self.data[-1].lat, self.data[-1].lon)

        # compute relative position
        t0 = self.point_time(idx)
        t1 = self.point_time(idx+1)
        if t0 >= t1:
            rel = 0.0
        else:
            rel = float(t - t0) / float(t1 - t0)

        # perform interpolation
        if (interpolate > 1) and ((idx < 1) or (idx >= (len(self.data) - 2))):
            interpolate = 1
        if not interpolate:
            if rel > 0.5:
                rel = 1.0
            else:
                rel = 0.0
        if rel < 0.01:
            return (self.data[idx].lat, self.data[idx].lon)
        if rel > 0.99:
            return (self.data[idx+1].lat, self.data[idx+1].lon)

        # cubic interpolation
        if interpolate == 2:
            return (hermite(self.data[idx-1].lat, self.data[idx].lat, self.data[idx+1].lat, self.data[idx+2].lat, rel), \
                    hermite(self.data[idx-1].lon, self.data[idx].lon, self.data[idx+1].lon, self.data[idx+2].lon, rel))

        # else: linear interpolation
        return (self.data[idx].lat + rel * (self.data[idx+1].lat - self.data[idx].lat), \
                self.data[idx].lon + rel * (self.data[idx+1].lon - self.data[idx].lon))

################################################################################

class GPXParser:
    point_types = set(['trkpt', 'wpt'])

    def __init__(self, filename, waypoints=None, track=None):
        if waypoints is None:
            waypoints = Track()
        if track is None:
            track = Track()
        self.waypoints = waypoints
        self.track = track
        self.current_point = None
        self.current_point_type = None
        self.current_tag = None
        self.current_data = None
        self.first_point = True
        self.filename = filename
        self.n = 0
        p = xml.parsers.expat.ParserCreate()
        p.StartElementHandler = self.start_element
        p.EndElementHandler = self.end_element
        p.CharacterDataHandler = self.cdata
        f = open(filename, "rb")
        p.ParseFile(f)
        f.close()
        self.track.sort()

    def start_element(self, name, attrs):
        if name in self.point_types:
            self.current_point_type = name
            try:
                self.current_point = TrackPoint(attrs['lat'], attrs['lon'], track_start=self.first_point)
            except KeyError:
                self.current_point = None
            except ValueError:
                print >>sys.stderr, "Invalid latitude/longitude in GPX file: %s / %s" % (attrs.get('lat', '?'), attrs.get('lon', '?'))
        else:
            self.current_tag = name
            self.current_data = ""

    def end_element(self, name):
        if (name == self.current_point_type) and self.current_point:
            self.current_point.finish()
            if (name == 'wpt'):
                self.waypoints.append(self.current_point)
            else:
                self.track.append(self.current_point)
            self.n += 1
            if not(self.n % 1000):
                sys.stderr.write("%s: %d points ...\r" % (self.filename, self.n))
                sys.stderr.flush()
            self.current_point = None
            self.first_point = False
        elif (name == self.current_tag) and self.current_point and self.current_data:
            self.handle_data()
        elif name == 'trk':
            self.first_point = True

    def cdata(self, data):
        self.current_data += data

    def get_data(self, name="value"):
        try:
            return float(self.current_data.strip())
        except ValueError:
            print >>sys.stderr, "Invalid", name, "in GPX file:", self.current_data.strip()

    def handle_data(self):
        if self.current_tag == 'time':
            try:
                self.current_point.t = ISO2Time(self.current_data)
            except ValueError:
                print >>sys.stderr, "Invalid time in GPX file:", self.current_data.strip()
        elif self.current_tag == 'ele':
            self.current_point.alt = self.get_data("altitude")
        elif self.current_tag == 'hdop':
            self.current_point.hdop = self.get_data("HDOP")
        elif self.current_tag == 'vdop':
            self.current_point.vdop = self.get_data("VDOP")
        elif self.current_tag == 'pdop':
            self.current_point.pdop = self.get_data("PDOP")

def ParseGPX(filename):
    gpx = GPXParser(filename)
    return gpx.waypoints, gpx.track

###############################################################################

def Waypoints2SplitInfo(waypoints):
    return [(p.t, None, None, None) for p in waypoints.data if p.has_time()]

def SplitTracks(source_track, info, use_subtrack_splits=True):
    if not source_track: return []
    tracks = []
    style = None
    current_track = None
    folder = None
    info.sort()
    source_track[0].track_start = use_subtrack_splits
    def append_current_track():
        if not(current_track) or (current_track.name == '-'):
            return
        if tracks and (current_track.name == '^'):
            tracks[-1].data.extend(current_track.data)
        else:
            tracks.append(current_track)
    for p in source_track:
        if info and (p.t >= info[0][0]):
            t, style, name, folder = info.pop(0)
            track_start = True
        else:
            t = None
            name = None
            track_start = p.track_start and use_subtrack_splits
        if track_start:
            append_current_track()
            current_track = Track()
            if t:      current_track.start_time = t
            if style:  current_track.style = style
            if name:   current_track.name = name
            if folder: current_track.folder = folder
        current_track.append(p)
    append_current_track()
    return tracks

###############################################################################

def SaveSplitFile(filename, tracks):
    f = open(filename, "w")
    for name, content in styles.iteritems():
        print >>f, '\\defstyle .%s %s' % (name, content)
    print >>f, '# format: timestamp [.style] [title]; if title is "-", track is skipped'
    for track in tracks:
        if track.start_time:
            print >>f, Time2ISO(track.start_time)
    f.close()

re_include = re.compile(r'\\include\s+(.*)')
re_folder = re.compile(r'\\folder\s*(.*)')
re_defstyle = re.compile(r'\\defstyle\s*\.([a-zA-Z0-9_-]+)\s(.*)')
re_split = re.compile(r'(\d\d\d\d-\d\d-\d\dT\d\d:\d\d:\d\dZ)\s*(\.[a-zA-Z0-9_-]+)?\s*(.*)')
re_encoding = re.compile(r'\\encoding\s*([a-zA-Z0-9_-]+)$')

def LoadSplitFile(filename, styles, include_visited=None, encoding=None, folder=None):
    if not include_visited:
        include_visited = set([])
    fullpath = os.path.abspath(filename)
    if fullpath in include_visited:
        return []
    include_visited.add(fullpath)
    info = []
    try:
        f = open(filename, "r")
    except IOError:
        print >>sys.stderr, "%s: could not open split list file" % filename
        return []
    if not encoding:
        encoding = sys.getfilesystemencoding()
    n = 0
    for line in f:
        n += 1
        line = line.split('#', 1)[0].strip()
        if not line:
            continue
        m = re_split.match(line)
        if m:
            style = m.group(2)
            if style:
                style = style[1:]
            info.append((ISO2Time(m.group(1)), style, unicode(m.group(3), encoding, 'replace'), folder))
            continue
        m = re_defstyle.match(line)
        if m:
            styles[m.group(1)] = m.group(2).strip()
            continue
        m = re_encoding.match(line)
        if m:
            encoding = m.group(1)
            continue
        m = re_include.match(line)
        if m:
            fn = os.path.normpath(os.path.join(os.path.dirname(filename), m.group(1)))
            info.extend(LoadSplitFile(fn, styles, include_visited, folder))
            continue
        m = re_folder.match(line)
        if m:
            folder = m.group(1).strip()
            if folder in ('', '-'):
                folder = None
            else:
                folder = unicode(folder, encoding, 'replace')
            continue
        print >>sys.stderr, "%s:%d: invalid line '%s'" % (filename, n, line)
    f.close()
    return info

###############################################################################

# photo list structure: tuple(path, width, height, lat, lon, time)

def GetJPEGSize(filename):
    f = open(filename, "rb")
    if f.read(2) != "\xFF\xD8":
        raise IOError("not a JPEG file")
    while True:
        marker = map(ord, f.read(4))
        if (len(marker) != 4) or (marker[0] != 0xFF):
            raise IOError("error in JPEG file")
        size = (marker[2] << 8) + marker[3] - 2
        if size < 0:
            raise IOError("invalid marker size in JPEG file")
        data = f.read(size)
        if len(data) != size:
            raise IOError("unexpected EOF in JPEG file")
        if marker[1] in (0xC5, 0xC6, 0xC7, 0xC8, 0xCD, 0xCE, 0xCF):
            raise IOError("unsupported JPEG format")
        if marker[1] in (0xC0, 0xC1, 0xC2, 0xC3, 0xC9, 0xCA, 0xCB):
            if size < 6:
                raise IOError("corrupted JPEG frame header")
            break
    f.close()
    return ((ord(data[3]) << 8) + ord(data[4]), (ord(data[1]) << 8) + ord(data[2]))

def ImportPhotos(photos, track, photodir, time_offset=0.0):
    dirpath = os.path.normpath(photodir)
    try:
        items = os.listdir(dirpath)
    except OSError, e:
        print >>sys.stderr, "Error: could not access directory", dirpath, "-", e
        return 0
    total_photos = 0
    for item in items:
        if item.startswith('.') or not(os.path.splitext(item)[-1].lower() in (".jpg", ".jpeg")):
            continue
        fullpath = os.path.normpath(os.path.join(photodir, item))
        try:
            s = os.stat(fullpath)
        except OSError, e:
            print >>sys.stderr, "Error: could not access", fullpath, "-", e
            continue
        if not stat.S_ISREG(s.st_mode):
            continue
        try:
            w, h = GetJPEGSize(fullpath)
        except IOError, e:
            print >>sys.stderr, "Error: could not read", fullpath, "-", e
            continue
        t = s.st_mtime + time_offset
        pos = track.point_at(t)
        if pos:
            photos.append((fullpath.replace(os.path.sep, '/'), w, h, pos[0], pos[1], t))
        total_photos += 1
    return total_photos

###############################################################################

def GenerateKML(filename, tracks, waypoints=None, photos=None, styles=DEFAULT_STYLES, photo_diag=DEFAULT_PHOTO_DIAG):
    f = open(filename, "w")
    print >>f, '<?xml version="1.0" encoding="UTF-8"?>'
    print >>f, '<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">'
    print >>f, '<Document>'
    print >>f, '\t<name>%s</name>' % os.path.splitext(os.path.basename(filename))[0]

    for name, content in styles.iteritems():
        print >>f, '\t<Style id="%s">' % name
        print >>f, '\t\t' + content
        print >>f, '\t</Style>'

    if waypoints:
        print >>f, '\t<Folder>'
        print >>f, '\t\t<name>Waypoints</name>'
        for p in waypoints:
            print >>f, '\t\t<Placemark>'
            if p.has_time():
                print >>f, '\t\t\t<name>%s</name>' % Time2Human(p.t)
            print >>f, '\t\t\t<Point><coordinates>%s,%s,%s</coordinates></Point>' % (p.lon, p.lat, p.alt)
            print >>f, '\t\t</Placemark>'
        print >>f, '\t</Folder>'

    print >>f, '\t<Folder>'
    print >>f, '\t\t<name>Track Polygons</name>'
    folder = None
    for track in tracks:
        if not track: continue
        if track.folder != folder:
            if folder:
                print >>f, '\t\t</Folder>'
            folder = track.folder
            if folder:
                print >>f, '\t\t<Folder><name>%s</name>' % X(folder)
        print >>f, '\t\t<Placemark>'
        print >>f, '\t\t\t<name>%s</name>' % X(track.get_name())
        print >>f, '\t\t\t<styleUrl>#%s</styleUrl>' % track.style
        print >>f, '\t\t\t<LineString>'
        print >>f, '\t\t\t\t<tessellate>1</tessellate>'
        print >>f, '\t\t\t\t<coordinates>'
        for p in track:
            print >>f, '\t\t\t\t\t%s,%s,%s' % (p.lon, p.lat, p.alt)
        print >>f, '\t\t\t\t</coordinates>'
        print >>f, '\t\t\t</LineString>'
        print >>f, '\t\t</Placemark>'
    if folder:
        print >>f, '\t\t</Folder>'
    print >>f, '\t</Folder>'

    print >>f, '\t<Folder>'
    print >>f, '\t\t<name>Tracks</name>'
    folder = None
    for track in tracks:
        if not track: continue
        if track.folder != folder:
            if folder:
                print >>f, '\t\t</Folder>'
            folder = track.folder
            if folder:
                print >>f, '\t\t<Folder><name>%s</name>' % X(folder)
        print >>f, '\t\t<Placemark>'
        print >>f, '\t\t\t<name>%s</name>' % X(track.get_name())
        print >>f, '\t\t\t<styleUrl>#%s</styleUrl>' % track.style
        print >>f, '\t\t\t<gx:Track>'
        for p in track:
            if p.has_time():
                print >>f, '\t\t\t\t<when>%s</when><gx:coord>%s %s %s</gx:coord>' % (Time2ISO(p.t), p.lon, p.lat, p.alt)
        print >>f, '\t\t\t</gx:Track>'
        print >>f, '\t\t</Placemark>'
    if folder:
        print >>f, '\t\t</Folder>'
    print >>f, '\t</Folder>'

    if photos:
        print >>f, '\t<Folder>'
        print >>f, '\t\t<name>Photos</name>'
        for p in photos:
            scale = float(photo_diag) / hypot(p[1], p[2])
            w = int(p[1] * scale + 0.5)
            h = int(p[2] * scale + 0.5)
            print >>f, '\t\t<Placemark>'
            print >>f, '\t\t\t<name>%s</name>' % (p[0].rsplit('/', 1)[-1].rsplit('.', 1)[0])
            print >>f, '\t\t\t<styleUrl>#' + PHOTO_STYLE_NAME + '</styleUrl>'
            print >>f, '\t\t\t<Point><coordinates>%s,%s,%s</coordinates></Point>' % (p[4], p[3], 0)
            print >>f, '\t\t\t<description><![CDATA['
            print >>f, '\t\t\t\t<img src="%s" width="%d" height="%d" />' % (p[0], w, h)
            print >>f, '\t\t\t]]></description>'
            print >>f, '\t\t</Placemark>'
        print >>f, '\t</Folder>'

    print >>f, '</Document>'
    print >>f, '</kml>'
    f.close()

###############################################################################

def _help(option, opt, value, parser):
    print __doc__.lstrip()
    parser.print_help()
    sys.exit(0)

if __name__ == "__main__":
    parser = optparse.OptionParser("%prog [OPTIONS] <input.gpx> <output.kml> [<splits.txt>]", add_help_option=False)
    parser.add_option("-h", "--help", help="show this help message and exit",
                      action="callback", callback=_help)
    parser.add_option("-s", "--simplify", dest="simplify_tolerance",
                      metavar="TOLERANCE", type="float", default=0.0,
                      help="simplify track with a tolerance of <TOLERANCE> meters")
    parser.add_option("-H", "--hdop", dest="hdop_influence",
                      metavar="PERCENT", type="int", default=50,
                      help="amount of influence of the HDOP value")
    parser.add_option("-p", "--photos", dest="photodirs",
                      metavar="DIR", action="append", type="string", default=[],
                      help="specify a directory with JPEG photos; may be specified multiple times; photos must have exact modification time stamps")
    parser.add_option("-t", "--offset", dest="offset",
                      metavar="SECONDS", type="float", default=0.0,
                      help="offset to add to photo times")
    parser.add_option("-d", "--diag", dest="diag",
                      metavar="PIXELS", type="int", default=DEFAULT_PHOTO_DIAG,
                      help="diagonal of photos in pop-ups [default: %default]")
    opts, args = parser.parse_args()
    
    gpx = []
    info_in = []
    info_out = None
    kml = None
    for arg in args:
        ext = os.path.splitext(arg)[-1].lower()
        if ext == ".gpx":
            gpx.append(arg)
        elif ext == ".kml":
            if kml:
                parser.error("multiple output KML files specified")
            kml = arg
        elif not os.path.isfile(arg):
            if info_in or info_out:
                parser.error("multiple output split files specified")
            info_out = arg
        else:
            if info_out:
                parser.error("multiple output split files specified")
            info_in.append(arg)
    if not gpx:
        parser.error("no input GPX file(s) specified")
    if not kml:
        parser.error("no output KML file specfied")

    waypoints = Track()
    full_track = Track()
    try:
        for fn in gpx:
            GPXParser(fn, waypoints, full_track)
    except IOError, e:
        print >>sys.stderr, "Input error:", e
        sys.exit(1)
    except KeyboardInterrupt:
        print >>sys.stderr, "Import cancelled. "
        sys.exit(0)
    print "imported %d waypoints and %d trackpoints in %d tracks" % (len(waypoints), len(full_track), full_track.count_subtracks())

    styles = DEFAULT_STYLES
    if info_in:
        print "splitting tracks by control file(s) ..."
        track_info = []
        for fn in info_in:  
            track_info.extend(LoadSplitFile(fn, styles))
        tracks = SplitTracks(full_track, track_info, use_subtrack_splits=False)
    else:
        print "splitting tracks by waypoints ..."
        tracks = SplitTracks(full_track, Waypoints2SplitInfo(waypoints))
    print "final result: %d tracks" % len(tracks)
    if info_out:
        print "saving split file '%s' ..." % info_out
        SaveSplitFile(info_out, tracks)

    if opts.simplify_tolerance:
        print "simplifying tracks ..."
        old_total = sum(map(len, tracks))
        for track in tracks:
            track.simplify(opts.simplify_tolerance, 0.01 * opts.hdop_influence)
        new_total = sum(map(len, tracks))
        if old_total == new_total:
            print "no points removed"
        else:
            delta = old_total - new_total
            print "reduced from %d to %d points (%d points = %.1f%% removed)" % (old_total, new_total, delta, 100.0 * delta / old_total)

    photos = []
    total_photos = 0
    for photodir in opts.photodirs:
        print "importing photos from", photodir, "..."
        total_photos += ImportPhotos(photos, full_track, photodir, opts.offset)
    if opts.photodirs:
        photos.sort(key=lambda x: x[5])
        print "imported", len(photos), "out of", total_photos, "photos"

    print "generating KML file '%s' ..." % kml
    try:
        GenerateKML(kml, tracks, waypoints, photos, styles, photo_diag=opts.diag)
    except IOError, e:
        print >>sys.stderr, "Output error:", e
        sys.exit(1)
    print "done."