4 ''' |
4 ''' |
5 gpx_reduce v1.8: removes points from gpx-files to reduce filesize and |
5 gpx_reduce v1.8: removes points from gpx-files to reduce filesize and |
6 tries to keep introduced distortions to the track at a minimum. |
6 tries to keep introduced distortions to the track at a minimum. |
7 Copyright (C) 2011,2012,2013,2015,2016,2017 travelling_salesman on OpenStreetMap |
7 Copyright (C) 2011,2012,2013,2015,2016,2017 travelling_salesman on OpenStreetMap |
8 |
8 |
9 changelog: v1.2: clarity refractoring + speedup for identical points |
9 changelog: v1.3: removal of all dependencies (lxml, scipy, numpy, pylab) |
|
10 v1.2: clarity refractoring + speedup for identical points |
10 v1.3: new track weighting functions, progress display |
11 v1.3: new track weighting functions, progress display |
11 v1.4: new track weighting function, restructuring for memory saving |
12 v1.4: new track weighting function, restructuring for memory saving |
12 v1.5: algorithm speedup by roughly a factor of 2 by eliminating some cases. |
13 v1.5: algorithm speedup by roughly a factor of 2 by eliminating some cases. |
13 v1.6: presets for train etc. |
14 v1.6: presets for train etc. |
14 v1.7: introduced weighting function for elevation errors |
15 v1.7: introduced weighting function for elevation errors |
26 |
27 |
27 You should have received a copy of the GNU General Public License |
28 You should have received a copy of the GNU General Public License |
28 along with this program. If not, see <http://www.gnu.org/licenses/>. |
29 along with this program. If not, see <http://www.gnu.org/licenses/>. |
29 ''' |
30 ''' |
30 |
31 |
31 import datetime |
32 import datetime, sys, time, operator |
32 import sys |
|
33 import time |
|
34 import numpy as np |
|
35 import numpy.linalg as la |
|
36 from math import * |
33 from math import * |
37 import xml.etree.ElementTree as etree |
34 import xml.etree.ElementTree as etree |
38 from optparse import OptionParser |
35 from optparse import OptionParser |
39 |
36 |
40 |
37 |
41 parser = OptionParser('usage: %prog [options] input-file.gpx') |
38 parser = OptionParser('usage: %prog [options] input-file.gpx') |
42 parser.add_option('-v', '--verbose', action='store', type='int', |
39 parser.add_option('-v', '--verbose', action='store', type='int', |
43 dest='verbose', default=1, help='verbose=[0,1]') |
40 dest='verbose', default=1, help='verbose=[0,1]') |
44 parser.add_option('-p', '--plot', action='store_true', dest='plot', |
|
45 default=False, help='Show a plot of the result at the end.') |
|
46 parser.add_option('-d', '--dist', action='store', type='float', dest='max_dist', |
41 parser.add_option('-d', '--dist', action='store', type='float', dest='max_dist', |
47 default=0.5, help='Maximum distance of line from original points in meters') |
42 default=0.5, help='Maximum distance of line from original points in meters') |
48 parser.add_option('-o', '--out', action='store', type='string', |
43 parser.add_option('-o', '--out', action='store', type='string', |
49 dest='ofname', default=None, help='Output file name') |
44 dest='ofname', default=None, help='Output file name') |
50 parser.add_option('-m', '--maxsep', action='store', type='float', dest='max_sep', |
45 parser.add_option('-m', '--maxsep', action='store', type='float', dest='max_sep', |
70 parser.add_option ('-r', dest='remove', default='extensions,hdop', help='remove tags T1,T2,..,Tn from every trackpoint', metavar='T1 T2 Tn') |
65 parser.add_option ('-r', dest='remove', default='extensions,hdop', help='remove tags T1,T2,..,Tn from every trackpoint', metavar='T1 T2 Tn') |
71 |
66 |
72 options, args = parser.parse_args() |
67 options, args = parser.parse_args() |
73 |
68 |
74 if len(args) < 1: |
69 if len(args) < 1: |
75 parser.print_usage() |
70 parser.print_help () |
76 exit(2) |
71 sys.exit (2) |
77 |
72 |
78 # use the WGS-84 ellipsoid |
73 # use the WGS-84 ellipsoid |
79 rE = 6356752.314245 # earth's radius |
74 rE = 6356752.314245 # earth's radius |
80 a = 6378137.0 |
75 a = 6378137.0 |
81 b = 6356752.314245179 |
76 b = 6356752.314245179 |
82 |
77 |
83 timeformat = '%Y-%m-%dT%H:%M:%SZ' |
78 timeformat = '%Y-%m-%dT%H:%M:%SZ' |
84 |
79 |
85 def distance(p1_, pm_, p2_, ele_weight=1.0): |
80 # the linear algebra with lists |
|
81 dot = lambda p1, p2: reduce (operator.add, [a * b for a, b in zip (p1, p2)]) |
|
82 vmin = lambda p1, p2: [a - b for a, b in zip (p1, p2)] |
|
83 smul = lambda p, c: [a * c for a in p] |
|
84 sdiv = lambda p, c: [a / c for a in p] |
|
85 norm = lambda p: sqrt (sum (a * a for a in p)) |
|
86 |
|
87 def distance (p1, pm, p2, ele_weight=1.0): |
86 # returns distance of pm from line between p1 and p2 |
88 # returns distance of pm from line between p1 and p2 |
87 |
89 |
88 p1, pm, p2 = np.array(p1_), np.array(pm_), np.array(p2_) |
90 h1, hm, h2 = norm (p1), norm (pm), norm (p2) |
89 h1, hm, h2 = la.norm(p1), la.norm(pm), la.norm(p2) |
|
90 if ele_weight != 1.0 and min(h1, hm, h2) > 0.0: |
91 if ele_weight != 1.0 and min(h1, hm, h2) > 0.0: |
91 hmean = (h1 + hm + h2) / 3.0 |
92 hmean = (h1 + hm + h2) / 3.0 |
92 p1 *= (ele_weight + (1.0 - ele_weight) * hmean / h1) |
93 p1 = smul (p1, ele_weight + (1.0 - ele_weight) * hmean / h1) |
93 pm *= (ele_weight + (1.0 - ele_weight) * hmean / hm) |
94 pm = smul (pm, ele_weight + (1.0 - ele_weight) * hmean / hm) |
94 p2 *= (ele_weight + (1.0 - ele_weight) * hmean / h2) |
95 p2 = smul (p2, ele_weight + (1.0 - ele_weight) * hmean / h2) |
95 line = p2 - p1 |
96 line = vmin (p2, p1) |
96 linel = la.norm(line) |
97 linel = norm (line) |
97 vm = pm - p1 |
98 vm = vmin (pm, p1) |
98 if linel == 0.0: |
99 if linel == 0.0: |
99 return la.norm(vm) |
100 return norm (vm) |
100 linem = line / linel |
101 linem = sdiv (line, linel) |
101 |
102 |
102 position = np.dot(vm, linem) / linel |
103 position = dot (vm, linem) / linel |
103 if position < 0.0: |
104 if position < 0.0: |
104 return la.norm(vm) |
105 return norm (vm) |
105 elif position > 1.0: |
106 elif position > 1.0: |
106 return la.norm(pm - p2) |
107 return norm (vmin (pm, p2)) |
107 else: |
108 else: |
108 return la.norm(vm - line * position) |
109 return norm (vmin (vm, smul (line, position))) |
109 |
110 |
110 def latlonele_to_xyz(lat, lon, ele): |
111 def latlonele_to_xyz (lat, lon, ele): |
111 s = sin(radians(lat)) |
112 s = sin (radians (lat)) |
112 c = cos(radians(lat)) |
113 c = cos (radians (lat)) |
113 r = ele + a * b / la.norm([s*a, c*b]) |
114 r = ele + a * b / norm ([s*a, c*b]) |
114 lon = radians(lon) |
115 lon = radians (lon) |
115 return r * c * sin(lon), r * c * (-cos(lon)), r * s |
116 return r * c * sin (lon), r * c * (-cos (lon)), r * s |
116 |
117 |
117 |
118 |
118 def reduced_track_indices(coordinate_list, timesteps=None): |
119 def reduced_track_indices(coordinate_list, timesteps=None): |
119 # returns a list of indices of trackpoints that constitute the reduced track |
120 # returns a list of indices of trackpoints that constitute the reduced track |
120 # takes a list of kartesian coordinate tuples |
121 # takes a list of kartesian coordinate tuples |
151 for i2 in range(1, n): |
152 for i2 in range(1, n): |
152 penalties = {} |
153 penalties = {} |
153 imin = None |
154 imin = None |
154 costmin = float('inf') |
155 costmin = float('inf') |
155 for i1 in reversed(range(i2)): |
156 for i1 in reversed(range(i2)): |
156 p1 = np.array(points[i1]['p']) |
157 p1 = points[i1]['p'] |
157 p2 = np.array(points[i2]['p']) |
158 p2 = points[i2]['p'] |
158 seglength = la.norm(p2 - p1) |
159 seglength = norm (vmin (p2, p1)) |
159 |
160 |
160 # estimate speed between p1 and p2 |
161 # estimate speed between p1 and p2 |
161 if timesteps != None: |
162 if timesteps != None: |
162 dt = (points[i2]['t'] - points[i1]['t']).total_seconds() |
163 dt = (points[i2]['t'] - points[i1]['t']).total_seconds() |
163 v = seglength / max(0.1, dt) |
164 v = seglength / max(0.1, dt) |
186 distance_squaremax = 0.0 |
187 distance_squaremax = 0.0 |
187 distance_squaresum = 0.0 |
188 distance_squaresum = 0.0 |
188 distances_squared = [] |
189 distances_squared = [] |
189 # iterate all medium points between i1 and i2 |
190 # iterate all medium points between i1 and i2 |
190 for im in range(i1+1, i2): |
191 for im in range(i1+1, i2): |
191 pm = np.array(points[im]['p']) |
192 pm = points[im]['p'] |
192 d = distance(p1, pm, p2, options.ele_weight) |
193 d = distance (p1, pm, p2, options.ele_weight) |
193 if (d <= options.max_dist): |
194 if (d <= options.max_dist): |
194 d_sq = (d / options.max_dist) ** 2 |
195 d_sq = (d / options.max_dist) ** 2 |
195 distance_squaremax = max(distance_squaremax, d_sq) |
196 distance_squaremax = max(distance_squaremax, d_sq) |
196 distance_squaresum += points[im]['weight'] * d_sq |
197 distance_squaresum += points[im]['weight'] * d_sq |
197 distances_squared.append(d_sq) |
198 distances_squared.append(d_sq) |
198 else: |
199 else: |
199 i1_i2_segment_valid = False |
200 i1_i2_segment_valid = False |
200 |
201 |
201 # check if connection to any further point i1 is impossible |
202 # check if connection to any further point i1 is impossible |
202 d1 = np.dot(p1 - p2, p1 - p2) |
203 d1 = dot (vmin (p1, p2), vmin (p1, p2)) |
203 d2 = np.dot(pm - p2, pm - p2) |
204 d2 = dot (vmin (pm, p2), vmin (pm, p2)) |
204 dd = options.max_dist ** 2 |
205 dd = options.max_dist ** 2 |
205 d1d2 = np.dot(p1 - p2, pm - p2) |
206 d1d2 = dot (vmin (p1, p2), vmin (pm, p2)) |
206 # formula from cosines of point separation angle and cone-opening angles around points |
207 # formula from cosines of point separation angle and cone-opening angles around points |
207 if (d1 > dd and d2 > dd and (d1d2 + dd)**2 < (d2 - dd) * (d1 - dd)): |
208 if (d1 > dd and d2 > dd and (d1d2 + dd)**2 < (d2 - dd) * (d1 - dd)): |
208 lower_i1_possible = False |
209 lower_i1_possible = False |
209 break |
210 break |
210 |
211 |
227 penalties[i1] = 0.0 |
228 penalties[i1] = 0.0 |
228 |
229 |
229 # add a penalty for kinks |
230 # add a penalty for kinks |
230 if options.bend > 0.: |
231 if options.bend > 0.: |
231 if points[i1]['prev'] != -1: |
232 if points[i1]['prev'] != -1: |
232 p0 = np.array(points[points[i1]['prev']]['p']) |
233 p0 = points [points [i1]['prev']] ['p'] |
233 v0 = p1 - p0 |
234 v0 = vmin (p1, p0) |
234 v1 = p2 - p1 |
235 v1 = vmin (p2, p1) |
235 if la.norm(v0) > 0. and la.norm(v1) > 0.: |
236 if norm (v0) > 0. and norm (v1) > 0.: |
236 v0 /= la.norm(v0) |
237 v0 = sdiv (v0, norm (v0)) |
237 v1 /= la.norm(v1) |
238 v1 = sdiv (v1, norm (v1)) |
238 kink = (1.0 - np.dot(v0, v1)) / 2.0 |
239 kink = (1.0 - dot (v0, v1)) / 2.0 |
239 penalties[i1] += options.bend * kink |
240 penalties[i1] += options.bend * kink |
240 |
241 |
241 # find best predecessor |
242 # find best predecessor |
242 imin = None |
243 imin = None |
243 costmin = float('inf') |
244 costmin = float('inf') |