--- a/gpx_reduce.py Wed Aug 23 22:02:19 2017 +0200
+++ b/gpx_reduce.py Mon Aug 28 16:24:53 2017 +0200
@@ -6,7 +6,8 @@
tries to keep introduced distortions to the track at a minimum.
Copyright (C) 2011,2012,2013,2015,2016,2017 travelling_salesman on OpenStreetMap
-changelog: v1.2: clarity refractoring + speedup for identical points
+changelog: v1.3: removal of all dependencies (lxml, scipy, numpy, pylab)
+ v1.2: clarity refractoring + speedup for identical points
v1.3: new track weighting functions, progress display
v1.4: new track weighting function, restructuring for memory saving
v1.5: algorithm speedup by roughly a factor of 2 by eliminating some cases.
@@ -28,11 +29,7 @@
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
+import datetime, sys, time, operator
from math import *
import xml.etree.ElementTree as etree
from optparse import OptionParser
@@ -41,8 +38,6 @@
parser = OptionParser('usage: %prog [options] input-file.gpx')
parser.add_option('-v', '--verbose', action='store', type='int',
dest='verbose', default=1, help='verbose=[0,1]')
-parser.add_option('-p', '--plot', action='store_true', dest='plot',
- default=False, help='Show a plot of the result at the end.')
parser.add_option('-d', '--dist', action='store', type='float', dest='max_dist',
default=0.5, help='Maximum distance of line from original points in meters')
parser.add_option('-o', '--out', action='store', type='string',
@@ -72,8 +67,8 @@
options, args = parser.parse_args()
if len(args) < 1:
- parser.print_usage()
- exit(2)
+ parser.print_help ()
+ sys.exit (2)
# use the WGS-84 ellipsoid
rE = 6356752.314245 # earth's radius
@@ -82,37 +77,43 @@
timeformat = '%Y-%m-%dT%H:%M:%SZ'
-def distance(p1_, pm_, p2_, ele_weight=1.0):
+# the linear algebra with lists
+dot = lambda p1, p2: reduce (operator.add, [a * b for a, b in zip (p1, p2)])
+vmin = lambda p1, p2: [a - b for a, b in zip (p1, p2)]
+smul = lambda p, c: [a * c for a in p]
+sdiv = lambda p, c: [a / c for a in p]
+norm = lambda p: sqrt (sum (a * a for a in p))
+
+def distance (p1, pm, p2, ele_weight=1.0):
# returns distance of pm from line between p1 and p2
- p1, pm, p2 = np.array(p1_), np.array(pm_), np.array(p2_)
- h1, hm, h2 = la.norm(p1), la.norm(pm), la.norm(p2)
+ h1, hm, h2 = norm (p1), norm (pm), norm (p2)
if ele_weight != 1.0 and min(h1, hm, h2) > 0.0:
hmean = (h1 + hm + h2) / 3.0
- p1 *= (ele_weight + (1.0 - ele_weight) * hmean / h1)
- pm *= (ele_weight + (1.0 - ele_weight) * hmean / hm)
- p2 *= (ele_weight + (1.0 - ele_weight) * hmean / h2)
- line = p2 - p1
- linel = la.norm(line)
- vm = pm - p1
+ p1 = smul (p1, ele_weight + (1.0 - ele_weight) * hmean / h1)
+ pm = smul (pm, ele_weight + (1.0 - ele_weight) * hmean / hm)
+ p2 = smul (p2, ele_weight + (1.0 - ele_weight) * hmean / h2)
+ line = vmin (p2, p1)
+ linel = norm (line)
+ vm = vmin (pm, p1)
if linel == 0.0:
- return la.norm(vm)
- linem = line / linel
+ return norm (vm)
+ linem = sdiv (line, linel)
- position = np.dot(vm, linem) / linel
+ position = dot (vm, linem) / linel
if position < 0.0:
- return la.norm(vm)
+ return norm (vm)
elif position > 1.0:
- return la.norm(pm - p2)
+ return norm (vmin (pm, p2))
else:
- return la.norm(vm - line * position)
+ return norm (vmin (vm, smul (line, position)))
-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 latlonele_to_xyz (lat, lon, ele):
+ s = sin (radians (lat))
+ c = cos (radians (lat))
+ r = ele + a * b / norm ([s*a, c*b])
+ lon = radians (lon)
+ return r * c * sin (lon), r * c * (-cos (lon)), r * s
def reduced_track_indices(coordinate_list, timesteps=None):
@@ -153,9 +154,9 @@
imin = None
costmin = float('inf')
for i1 in reversed(range(i2)):
- p1 = np.array(points[i1]['p'])
- p2 = np.array(points[i2]['p'])
- seglength = la.norm(p2 - p1)
+ p1 = points[i1]['p']
+ p2 = points[i2]['p']
+ seglength = norm (vmin (p2, p1))
# estimate speed between p1 and p2
if timesteps != None:
@@ -188,8 +189,8 @@
distances_squared = []
# iterate all medium points between i1 and i2
for im in range(i1+1, i2):
- pm = np.array(points[im]['p'])
- d = distance(p1, pm, p2, options.ele_weight)
+ pm = points[im]['p']
+ d = distance (p1, pm, p2, options.ele_weight)
if (d <= options.max_dist):
d_sq = (d / options.max_dist) ** 2
distance_squaremax = max(distance_squaremax, d_sq)
@@ -199,10 +200,10 @@
i1_i2_segment_valid = False
# check if connection to any further point i1 is impossible
- d1 = np.dot(p1 - p2, p1 - p2)
- d2 = np.dot(pm - p2, pm - p2)
+ d1 = dot (vmin (p1, p2), vmin (p1, p2))
+ d2 = dot (vmin (pm, p2), vmin (pm, p2))
dd = options.max_dist ** 2
- d1d2 = np.dot(p1 - p2, pm - p2)
+ d1d2 = dot (vmin (p1, p2), vmin (pm, p2))
# formula from cosines of point separation angle and cone-opening angles around points
if (d1 > dd and d2 > dd and (d1d2 + dd)**2 < (d2 - dd) * (d1 - dd)):
lower_i1_possible = False
@@ -229,13 +230,13 @@
# add a penalty for kinks
if options.bend > 0.:
if points[i1]['prev'] != -1:
- p0 = np.array(points[points[i1]['prev']]['p'])
- v0 = p1 - p0
- v1 = p2 - p1
- if la.norm(v0) > 0. and la.norm(v1) > 0.:
- v0 /= la.norm(v0)
- v1 /= la.norm(v1)
- kink = (1.0 - np.dot(v0, v1)) / 2.0
+ p0 = points [points [i1]['prev']] ['p']
+ v0 = vmin (p1, p0)
+ v1 = vmin (p2, p1)
+ if norm (v0) > 0. and norm (v1) > 0.:
+ v0 = sdiv (v0, norm (v0))
+ v1 = sdiv (v1, norm (v1))
+ kink = (1.0 - dot (v0, v1)) / 2.0
penalties[i1] += options.bend * kink
# find best predecessor