80 a = 6378137.0 |
80 a = 6378137.0 |
81 b = 6356752.314245179 |
81 b = 6356752.314245179 |
82 |
82 |
83 timeformat = '%Y-%m-%dT%H:%M:%SZ' |
83 timeformat = '%Y-%m-%dT%H:%M:%SZ' |
84 |
84 |
85 |
|
86 def distance(p1_, pm_, p2_, ele_weight=1.0): |
85 def distance(p1_, pm_, p2_, ele_weight=1.0): |
87 # returns distance of pm from line between p1 and p2 |
86 # returns distance of pm from line between p1 and p2 |
88 |
87 |
89 p1, pm, p2 = np.array(p1_), np.array(pm_), np.array(p2_) |
88 p1, pm, p2 = np.array(p1_), np.array(pm_), np.array(p2_) |
90 h1, hm, h2 = la.norm(p1), la.norm(pm), la.norm(p2) |
89 h1, hm, h2 = la.norm(p1), la.norm(pm), la.norm(p2) |
106 elif position > 1.0: |
105 elif position > 1.0: |
107 return la.norm(pm - p2) |
106 return la.norm(pm - p2) |
108 else: |
107 else: |
109 return la.norm(vm - line * position) |
108 return la.norm(vm - line * position) |
110 |
109 |
111 |
|
112 def rotate(x, y, phi): |
|
113 return x*cos(phi) - y*sin(phi), x*sin(phi) + y*cos(phi) |
|
114 |
|
115 |
|
116 def project_to_meters(lat, lon, latm, lonm): |
|
117 # azimuthal map projection centered at average track coordinate |
|
118 lon -= lonm |
|
119 xyz = latlonele_to_xyz(lat, lon, 0.0) |
|
120 zy = rotate(xyz[2], xyz[1], radians(90 - latm)) |
|
121 lat2 = atan2(zy[0], la.norm([zy[1], xyz[0]])) |
|
122 lon2 = atan2(xyz[0], -zy[1]) |
|
123 x_meters = rE * sin(lon2) * (pi / 2.0 - lat2) |
|
124 y_meters = -rE * cos(lon2) * (pi / 2.0 - lat2) |
|
125 return x_meters, y_meters |
|
126 |
|
127 |
|
128 def latlonele_to_xyz(lat, lon, ele): |
110 def latlonele_to_xyz(lat, lon, ele): |
129 s = sin(radians(lat)) |
111 s = sin(radians(lat)) |
130 c = cos(radians(lat)) |
112 c = cos(radians(lat)) |
131 r = ele + a * b / la.norm([s*a, c*b]) |
113 r = ele + a * b / la.norm([s*a, c*b]) |
132 lon = radians(lon) |
114 lon = radians(lon) |
133 return r * c * sin(lon), r * c * (-cos(lon)), r * s |
115 return r * c * sin(lon), r * c * (-cos(lon)), r * s |
134 |
|
135 |
|
136 def xyz_to_latlonele(x, y, z): |
|
137 r = la.norm([x, y, z]) |
|
138 if (r == 0): |
|
139 return 0.0, 0.0, 0.0 |
|
140 lat = degrees(atan2(z, la.norm([x, y]))) |
|
141 lon = degrees(atan2(x, -y)) |
|
142 ele = r * (1.0 - a * b / la.norm([a*z, b*x, b*y])) |
|
143 return lat, lon, ele |
|
144 |
116 |
145 |
117 |
146 def reduced_track_indices(coordinate_list, timesteps=None): |
118 def reduced_track_indices(coordinate_list, timesteps=None): |
147 # returns a list of indices of trackpoints that constitute the reduced track |
119 # returns a list of indices of trackpoints that constitute the reduced track |
148 # takes a list of kartesian coordinate tuples |
120 # takes a list of kartesian coordinate tuples |
301 |
273 |
302 |
274 |
303 ############################## main function ################################# |
275 ############################## main function ################################# |
304 for fname in args: |
276 for fname in args: |
305 # initialisations |
277 # initialisations |
306 tracksegs_old = [] |
278 ntot = 0 # total number of trackpoints (sum of segments) |
307 tracksegs_new = [] |
279 newtot = 0 # total number of trackpoints after removal |
308 sumx, sumy, sumz = 0.0, 0.0, 0.0 |
|
309 |
280 |
310 # import xml data from files |
281 # import xml data from files |
311 print 'opening file', fname |
282 print 'opening file', fname |
312 infile = open(fname) |
283 infile = open(fname) |
313 tree = etree.parse(infile) |
284 tree = etree.parse(infile) |
317 nsurl = gpx.tag.split ('}')[0][1:] # == 'http://www.topografix.com/GPX/1/1' |
288 nsurl = gpx.tag.split ('}')[0][1:] # == 'http://www.topografix.com/GPX/1/1' |
318 etree.register_namespace('', nsurl) # default namespace -> xmlns:.... in the output |
289 etree.register_namespace('', nsurl) # default namespace -> xmlns:.... in the output |
319 nsmap = '{' + nsurl + '}' # prefix for all tags in the tree |
290 nsmap = '{' + nsurl + '}' # prefix for all tags in the tree |
320 |
291 |
321 # extract data from xml |
292 # extract data from xml |
322 for trkseg in gpx.findall('.//' + nsmap + 'trkseg'): |
293 for si, trkseg in enumerate (gpx.findall('.//' + nsmap + 'trkseg')): |
323 trkpts = trkseg.findall(nsmap + 'trkpt') |
294 trkpts = trkseg.findall(nsmap + 'trkpt') |
324 n = len(trkpts) |
295 n = len(trkpts) |
325 |
296 |
326 # extract coordinate values |
297 # extract coordinate values |
327 lats = [float(trkpt.get('lat')) for trkpt in trkpts] |
298 lats = [float(trkpt.get('lat')) for trkpt in trkpts] |
329 eles = [float(trkpt.find(nsmap + 'ele').text) for trkpt in trkpts] |
300 eles = [float(trkpt.find(nsmap + 'ele').text) for trkpt in trkpts] |
330 try: |
301 try: |
331 times = [datetime.datetime.strptime(trkpt.find(nsmap + 'time' |
302 times = [datetime.datetime.strptime(trkpt.find(nsmap + 'time' |
332 ).text, timeformat) for trkpt in trkpts] |
303 ).text, timeformat) for trkpt in trkpts] |
333 except Exception as e: |
304 except Exception as e: |
334 print e |
305 print '-- trackpoint without time' |
335 times = None |
306 times = None |
336 |
307 |
337 # save original trackseg for plotting |
|
338 if options.plot: |
|
339 tracksegs_old.append([[lats[i], lons[i], eles[i]] for i in range(n)]) |
|
340 |
|
341 # calculate projected points to work on |
308 # calculate projected points to work on |
342 coords = [] |
309 coords = [] |
343 for i in range(n): |
310 for i in range(n): |
344 x, y, z = latlonele_to_xyz(lats[i], lons[i], eles[i]) |
311 x, y, z = latlonele_to_xyz (lats[i], lons[i], eles[i]) |
345 coords.append((x, y, z)) |
312 coords.append((x, y, z)) |
346 sumx += x |
|
347 sumy += y |
|
348 sumz += z |
|
349 |
313 |
350 # execute the reduction algorithm |
314 # execute the reduction algorithm |
351 final_pnums = reduced_track_indices(coords, times) |
315 final_pnums = reduced_track_indices(coords, times) |
352 |
316 |
353 n_new = len(final_pnums) |
317 n_new = len (final_pnums) |
354 print 'number of points:', n, '-', n - n_new, '=', n_new |
318 print 'segment %d, with %d - %d = %d points' % (si, n, n - n_new, n_new) |
355 |
319 ntot += n |
|
320 newtot += n_new |
|
321 |
356 # delete certain points from original data |
322 # delete certain points from original data |
357 delete_pnums = [i for i in range(n) if i not in final_pnums] |
323 delete_pnums = [i for i in range(n) if i not in final_pnums] |
358 for i in reversed(delete_pnums): |
324 for i in reversed(delete_pnums): |
359 trkseg.remove (trkpts[i]) # remove from the xml-tree |
325 trkseg.remove (trkpts[i]) # remove from the xml-tree |
360 del trkpts [i] # also remove from the list |
326 del trkpts [i] # also remove from the list |
366 for trkpnt in trkpts: |
332 for trkpnt in trkpts: |
367 for tag in taglist: |
333 for tag in taglist: |
368 e = trkpnt.find (nsmap + tag) |
334 e = trkpnt.find (nsmap + tag) |
369 if e != None: trkpnt.remove (e) |
335 if e != None: trkpnt.remove (e) |
370 |
336 |
371 # save reduced trackseg for plotting |
337 print 'total number of points: %d - %d = %d' % (ntot, ntot - newtot, newtot) |
372 if options.plot: |
338 |
373 tracksegs_new.append ([[float(trkpt.get('lat')), |
|
374 float(trkpt.get('lon')), float(trkpt.find(nsmap + 'ele').text)] |
|
375 for trkpt in trkseg.findall(nsmap + 'trkpt')]) |
|
376 |
|
377 # export data to file |
339 # export data to file |
378 if options.ofname != None: |
340 if options.ofname != None: |
379 ofname = options.ofname |
341 ofname = options.ofname |
380 elif fname.endswith('.gpx'): |
342 elif fname.endswith('.gpx'): |
381 ofname = fname[:-4] + '_reduced.gpx' |
343 ofname = fname[:-4] + '_reduced.gpx' |
383 ofname = fname + '_reduced.gpx' |
345 ofname = fname + '_reduced.gpx' |
384 outfile = open(ofname, 'w') |
346 outfile = open(ofname, 'w') |
385 tree.write (outfile, encoding="utf-8", xml_declaration=True, default_namespace=None, method="xml") |
347 tree.write (outfile, encoding="utf-8", xml_declaration=True, default_namespace=None, method="xml") |
386 outfile.close() |
348 outfile.close() |
387 print 'modified copy written to', ofname |
349 print 'modified copy written to', ofname |
388 |
|
389 # plot result to screen |
|
390 if options.plot: |
|
391 latm, lonm, elesum = xyz_to_latlonele(sumx, sumy, sumz) |
|
392 |
|
393 data_old = [] |
|
394 for trkseg in tracksegs_old: |
|
395 for trkpt in trkseg: |
|
396 xy = project_to_meters(trkpt[0], trkpt[1], latm, lonm) |
|
397 data_old.append (xy) |
|
398 |
|
399 data = [] |
|
400 for trkseg in tracksegs_new: |
|
401 for trkpt in trkseg: |
|
402 xy = project_to_meters(trkpt[0], trkpt[1], latm, lonm) |
|
403 data.append (xy) |
|
404 |
|
405 from subprocess import Popen, PIPE |
|
406 gnuPlotCmd = ['gnuplot'] |
|
407 plot = Popen (gnuPlotCmd, stdin=PIPE, stdout=PIPE, stderr=PIPE) |
|
408 plot.stdin.write (b"plot '-' with linespoints lc rgb 'red' pt 4, '-' with linespoints pt 6\n") |
|
409 plot.stdin.write ("\n".join ('%f %f' % d for d in data).encode ()) |
|
410 plot.stdin.write (b'\ne\n') |
|
411 plot.stdin.write ("\n".join ('%f %f' % d for d in data_old).encode ()) |
|
412 plot.stdin.write (b'\ne\n') |
|
413 plot.stdin.flush () |
|
414 raw_input ('druk') |
|