diff -r f94099a0277a -r 1b96bb9de1f3 gpx_plot.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/gpx_plot.py	Wed Aug 23 21:25:58 2017 +0200
@@ -0,0 +1,168 @@
+#!/usr/bin/env python
+# -*- coding: utf8 -*-
+
+'''
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+'''
+
+import datetime
+import sys
+import time
+import numpy as np
+import numpy.linalg as la
+from math import *
+import xml.etree.ElementTree as etree
+from optparse import OptionParser
+
+
+parser = OptionParser('usage: %prog [options] input-file.gpx')
+(options, args) = parser.parse_args()
+
+
+if len(args) < 1:
+    parser.print_usage()
+    exit(2)
+
+
+# use the WGS-84 ellipsoid
+rE = 6356752.314245 # earth's radius
+a = 6378137.0
+b = 6356752.314245179
+
+timeformat = '%Y-%m-%dT%H:%M:%SZ'
+
+def rotate(x, y, phi):
+    return x*cos(phi) - y*sin(phi), x*sin(phi) + y*cos(phi)
+
+
+def project_to_meters(lat, lon, latm, lonm):
+    # azimuthal map projection centered at average track coordinate
+    lon -= lonm
+    xyz = latlonele_to_xyz(lat, lon, 0.0)
+    zy = rotate(xyz[2], xyz[1], radians(90 - latm))
+    lat2 = atan2(zy[0], la.norm([zy[1], xyz[0]]))
+    lon2 = atan2(xyz[0], -zy[1])
+    x_meters = rE * sin(lon2) * (pi / 2.0 - lat2)
+    y_meters = -rE * cos(lon2) * (pi / 2.0 - lat2)
+    return x_meters, y_meters
+
+
+def latlonele_to_xyz(lat, lon, ele):
+    s = sin(radians(lat))
+    c = cos(radians(lat))
+    r = ele + a * b / la.norm([s*a, c*b])
+    lon = radians(lon)
+    return r * c * sin(lon), r * c * (-cos(lon)), r * s
+
+
+def xyz_to_latlonele(x, y, z):
+    r = la.norm([x, y, z])
+    if (r == 0):
+        return 0.0, 0.0, 0.0
+    lat = degrees(atan2(z, la.norm([x, y])))
+    lon = degrees(atan2(x, -y))
+    ele = r * (1.0 - a * b / la.norm([a*z, b*x, b*y]))
+    return lat, lon, ele
+
+
+############################## main function #################################
+tracks = []
+npoints = []
+for fname in args:
+    # initialisations
+    tracksegs = []
+    sumx, sumy, sumz = 0.0, 0.0, 0.0
+    ntot = 0    # total number of trackpoints (sum of segments)
+    
+    # import xml data from files
+    print 'opening file', fname
+    infile = open(fname)
+    tree = etree.parse(infile)
+    infile.close()
+
+    gpx = tree.getroot()
+    nsurl = gpx.tag.split ('}')[0][1:]  # == 'http://www.topografix.com/GPX/1/1'
+    etree.register_namespace('', nsurl) # default namespace -> xmlns:.... in the output
+    nsmap = '{' + nsurl + '}'           # prefix for all tags in the tree
+
+    # extract data from xml
+    for si, trkseg in enumerate (gpx.findall('.//' + nsmap + 'trkseg')):
+        trkpts = trkseg.findall(nsmap + 'trkpt')
+        n = len(trkpts)
+        
+        # extract coordinate values
+        lats = [float(trkpt.get('lat')) for trkpt in trkpts]
+        lons = [float(trkpt.get('lon')) for trkpt in trkpts]
+        eles = [float(trkpt.find(nsmap + 'ele').text) for trkpt in trkpts]
+        try:
+            times = [datetime.datetime.strptime(trkpt.find(nsmap + 'time'
+                                       ).text, timeformat) for trkpt in trkpts]
+        except Exception as e:
+            print '-- trackpoint without time'
+            times = None
+
+        # save original trackseg for plotting
+        tracksegs.append ([ (lats[i], lons[i], eles[i]) for i in range(n) ])
+
+        # calculate projected points to work on
+        for i in range(n):
+            x, y, z = latlonele_to_xyz (lats[i], lons[i], eles[i])
+            sumx += x
+            sumy += y
+            sumz += z
+
+        print 'segment %d, with %d points:' % (si, n)
+        ntot += n
+
+    tracks.append (tracksegs)
+    npoints.append ((fname.replace ('_','\_'), ntot))   # underscore -> subscripts in gnuplot
+    print 'number of points in track %s = %d:' % (fname, ntot)
+
+latm, lonm, elesum = xyz_to_latlonele (sumx, sumy, sumz)
+
+data = []
+xmin, xmax = ymin, ymax = float ('inf'), float ('-inf')
+for tracksegs in tracks:
+    for trkseg in tracksegs:
+        for lat, lon, ele in trkseg:
+            x, y = project_to_meters (lat, lon, latm, lonm)
+            data.append ('%f %f' % (x, y))
+            xmin, ymin = min (xmin, x), min (ymin, y)   # determine the x range
+            xmax, ymax = max (xmax, x), max (ymax, y)   # and the y range
+    data.append ('e')
+
+dx, dy = xmax - xmin, ymax - ymin   # make x and y ranges equal to the largest
+if dx > dy:                         # and keep ranges centered
+    dr = (dx - dy) / 2
+    ymax += dr
+    ymin -= dr
+else:
+    dr = (dy - dx) / 2
+    xmax += dr
+    xmin -= dr
+
+from subprocess import Popen, PIPE
+gnuPlotCmd = ['gnuplot']
+plot = Popen (gnuPlotCmd, stdin=PIPE, stdout=PIPE, stderr=PIPE)
+
+range = 'set xrange [%f:%f]\nset yrange [%f:%f]\n' % (xmin, xmax, ymin, ymax)
+plot.stdin.write (range)
+
+curves = ','.join ("'-' with linespoints ti '%s: %d punten'" % t for t in npoints)
+plot.stdin.write ('plot ' + curves + '\n')
+
+plot.stdin.write ("\n".join (data))
+plot.stdin.write ('\n')
+plot.stdin.flush ()
+raw_input ('druk')