--- 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')