--- a/gpx_reduce.py	Wed Aug 23 21:25:58 2017 +0200
+++ b/gpx_reduce.py	Wed Aug 23 22:02:19 2017 +0200
@@ -82,7 +82,6 @@
 
 timeformat = '%Y-%m-%dT%H:%M:%SZ'
 
-
 def distance(p1_, pm_, p2_, ele_weight=1.0):
     # returns distance of pm from line between p1 and p2
 
@@ -108,23 +107,6 @@
     else:
         return la.norm(vm - line * position)
 
-
-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))
@@ -133,16 +115,6 @@
     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
-
-
 def reduced_track_indices(coordinate_list, timesteps=None):
     # returns a list of indices of trackpoints that constitute the reduced track
     # takes a list of kartesian coordinate tuples
@@ -303,9 +275,8 @@
 ############################## main function #################################
 for fname in args:
     # initialisations
-    tracksegs_old = []
-    tracksegs_new = []
-    sumx, sumy, sumz = 0.0, 0.0, 0.0
+    ntot = 0    # total number of trackpoints (sum of segments)
+    newtot = 0 # total number of trackpoints after removal
     
     # import xml data from files
     print 'opening file', fname
@@ -319,7 +290,7 @@
     nsmap = '{' + nsurl + '}'           # prefix for all tags in the tree
 
     # extract data from xml
-    for trkseg in gpx.findall('.//' + nsmap + 'trkseg'):
+    for si, trkseg in enumerate (gpx.findall('.//' + nsmap + 'trkseg')):
         trkpts = trkseg.findall(nsmap + 'trkpt')
         n = len(trkpts)
         
@@ -331,28 +302,23 @@
             times = [datetime.datetime.strptime(trkpt.find(nsmap + 'time'
                                        ).text, timeformat) for trkpt in trkpts]
         except Exception as e:
-            print e
+            print '-- trackpoint without time'
             times = None
 
-        # save original trackseg for plotting
-        if options.plot:
-            tracksegs_old.append([[lats[i], lons[i], eles[i]] for i in range(n)])
-        
         # calculate projected points to work on
         coords = []
         for i in range(n):
-            x, y, z = latlonele_to_xyz(lats[i], lons[i], eles[i])
+            x, y, z = latlonele_to_xyz (lats[i], lons[i], eles[i])
             coords.append((x, y, z))
-            sumx += x
-            sumy += y
-            sumz += z
         
         # execute the reduction algorithm
         final_pnums = reduced_track_indices(coords, times)
-        
-        n_new = len(final_pnums)
-        print 'number of points:', n, '-', n - n_new, '=', n_new
-        
+
+        n_new = len (final_pnums)
+        print 'segment %d, with %d - %d = %d points' % (si, n, n - n_new, n_new)
+        ntot += n
+        newtot += n_new
+
         # delete certain points from original data
         delete_pnums = [i for i in range(n) if i not in final_pnums]
         for i in reversed(delete_pnums):
@@ -368,12 +334,8 @@
                 e = trkpnt.find (nsmap + tag)
                 if e != None: trkpnt.remove (e)
 
-        # save reduced trackseg for plotting
-        if options.plot:
-            tracksegs_new.append ([[float(trkpt.get('lat')),
-                float(trkpt.get('lon')), float(trkpt.find(nsmap + 'ele').text)]
-                for trkpt in trkseg.findall(nsmap + 'trkpt')])
-    
+    print 'total number of points: %d - %d = %d' % (ntot, ntot - newtot, newtot)
+
     # export data to file
     if options.ofname != None:
         ofname = options.ofname
@@ -385,30 +347,3 @@
     tree.write (outfile, encoding="utf-8", xml_declaration=True, default_namespace=None, method="xml")
     outfile.close()
     print 'modified copy written to', ofname
-
-    # plot result to screen
-    if options.plot:
-        latm, lonm, elesum = xyz_to_latlonele(sumx, sumy, sumz)
-
-        data_old = []
-        for trkseg in tracksegs_old:
-            for trkpt in trkseg:
-                xy = project_to_meters(trkpt[0], trkpt[1], latm, lonm)
-                data_old.append (xy)
-        
-        data = []
-        for trkseg in tracksegs_new:
-            for trkpt in trkseg:
-                xy = project_to_meters(trkpt[0], trkpt[1], latm, lonm)
-                data.append (xy)
-
-        from subprocess import Popen, PIPE
-        gnuPlotCmd = ['gnuplot']
-        plot = Popen (gnuPlotCmd, stdin=PIPE, stdout=PIPE, stderr=PIPE)
-        plot.stdin.write (b"plot '-' with linespoints lc rgb 'red' pt 4, '-' with linespoints pt 6\n")
-        plot.stdin.write ("\n".join ('%f %f' % d for d in data).encode ())
-        plot.stdin.write (b'\ne\n')
-        plot.stdin.write ("\n".join ('%f %f' % d for d in data_old).encode ())
-        plot.stdin.write (b'\ne\n')
-        plot.stdin.flush ()
-        raw_input ('druk')