|
1 #!/usr/bin/env python |
|
2 # -*- coding: utf8 -*- |
|
3 |
|
4 ''' |
|
5 This program is free software: you can redistribute it and/or modify |
|
6 it under the terms of the GNU General Public License as published by |
|
7 the Free Software Foundation, either version 3 of the License, or |
|
8 (at your option) any later version. |
|
9 |
|
10 This program is distributed in the hope that it will be useful, |
|
11 but WITHOUT ANY WARRANTY; without even the implied warranty of |
|
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|
13 GNU General Public License for more details. |
|
14 |
|
15 You should have received a copy of the GNU General Public License |
|
16 along with this program. If not, see <http://www.gnu.org/licenses/>. |
|
17 ''' |
|
18 |
|
19 import datetime |
|
20 import sys |
|
21 import time |
|
22 import numpy as np |
|
23 import numpy.linalg as la |
|
24 from math import * |
|
25 import xml.etree.ElementTree as etree |
|
26 from optparse import OptionParser |
|
27 |
|
28 |
|
29 parser = OptionParser('usage: %prog [options] input-file.gpx') |
|
30 (options, args) = parser.parse_args() |
|
31 |
|
32 |
|
33 if len(args) < 1: |
|
34 parser.print_usage() |
|
35 exit(2) |
|
36 |
|
37 |
|
38 # use the WGS-84 ellipsoid |
|
39 rE = 6356752.314245 # earth's radius |
|
40 a = 6378137.0 |
|
41 b = 6356752.314245179 |
|
42 |
|
43 timeformat = '%Y-%m-%dT%H:%M:%SZ' |
|
44 |
|
45 def rotate(x, y, phi): |
|
46 return x*cos(phi) - y*sin(phi), x*sin(phi) + y*cos(phi) |
|
47 |
|
48 |
|
49 def project_to_meters(lat, lon, latm, lonm): |
|
50 # azimuthal map projection centered at average track coordinate |
|
51 lon -= lonm |
|
52 xyz = latlonele_to_xyz(lat, lon, 0.0) |
|
53 zy = rotate(xyz[2], xyz[1], radians(90 - latm)) |
|
54 lat2 = atan2(zy[0], la.norm([zy[1], xyz[0]])) |
|
55 lon2 = atan2(xyz[0], -zy[1]) |
|
56 x_meters = rE * sin(lon2) * (pi / 2.0 - lat2) |
|
57 y_meters = -rE * cos(lon2) * (pi / 2.0 - lat2) |
|
58 return x_meters, y_meters |
|
59 |
|
60 |
|
61 def latlonele_to_xyz(lat, lon, ele): |
|
62 s = sin(radians(lat)) |
|
63 c = cos(radians(lat)) |
|
64 r = ele + a * b / la.norm([s*a, c*b]) |
|
65 lon = radians(lon) |
|
66 return r * c * sin(lon), r * c * (-cos(lon)), r * s |
|
67 |
|
68 |
|
69 def xyz_to_latlonele(x, y, z): |
|
70 r = la.norm([x, y, z]) |
|
71 if (r == 0): |
|
72 return 0.0, 0.0, 0.0 |
|
73 lat = degrees(atan2(z, la.norm([x, y]))) |
|
74 lon = degrees(atan2(x, -y)) |
|
75 ele = r * (1.0 - a * b / la.norm([a*z, b*x, b*y])) |
|
76 return lat, lon, ele |
|
77 |
|
78 |
|
79 ############################## main function ################################# |
|
80 tracks = [] |
|
81 npoints = [] |
|
82 for fname in args: |
|
83 # initialisations |
|
84 tracksegs = [] |
|
85 sumx, sumy, sumz = 0.0, 0.0, 0.0 |
|
86 ntot = 0 # total number of trackpoints (sum of segments) |
|
87 |
|
88 # import xml data from files |
|
89 print 'opening file', fname |
|
90 infile = open(fname) |
|
91 tree = etree.parse(infile) |
|
92 infile.close() |
|
93 |
|
94 gpx = tree.getroot() |
|
95 nsurl = gpx.tag.split ('}')[0][1:] # == 'http://www.topografix.com/GPX/1/1' |
|
96 etree.register_namespace('', nsurl) # default namespace -> xmlns:.... in the output |
|
97 nsmap = '{' + nsurl + '}' # prefix for all tags in the tree |
|
98 |
|
99 # extract data from xml |
|
100 for si, trkseg in enumerate (gpx.findall('.//' + nsmap + 'trkseg')): |
|
101 trkpts = trkseg.findall(nsmap + 'trkpt') |
|
102 n = len(trkpts) |
|
103 |
|
104 # extract coordinate values |
|
105 lats = [float(trkpt.get('lat')) for trkpt in trkpts] |
|
106 lons = [float(trkpt.get('lon')) for trkpt in trkpts] |
|
107 eles = [float(trkpt.find(nsmap + 'ele').text) for trkpt in trkpts] |
|
108 try: |
|
109 times = [datetime.datetime.strptime(trkpt.find(nsmap + 'time' |
|
110 ).text, timeformat) for trkpt in trkpts] |
|
111 except Exception as e: |
|
112 print '-- trackpoint without time' |
|
113 times = None |
|
114 |
|
115 # save original trackseg for plotting |
|
116 tracksegs.append ([ (lats[i], lons[i], eles[i]) for i in range(n) ]) |
|
117 |
|
118 # calculate projected points to work on |
|
119 for i in range(n): |
|
120 x, y, z = latlonele_to_xyz (lats[i], lons[i], eles[i]) |
|
121 sumx += x |
|
122 sumy += y |
|
123 sumz += z |
|
124 |
|
125 print 'segment %d, with %d points:' % (si, n) |
|
126 ntot += n |
|
127 |
|
128 tracks.append (tracksegs) |
|
129 npoints.append ((fname.replace ('_','\_'), ntot)) # underscore -> subscripts in gnuplot |
|
130 print 'number of points in track %s = %d:' % (fname, ntot) |
|
131 |
|
132 latm, lonm, elesum = xyz_to_latlonele (sumx, sumy, sumz) |
|
133 |
|
134 data = [] |
|
135 xmin, xmax = ymin, ymax = float ('inf'), float ('-inf') |
|
136 for tracksegs in tracks: |
|
137 for trkseg in tracksegs: |
|
138 for lat, lon, ele in trkseg: |
|
139 x, y = project_to_meters (lat, lon, latm, lonm) |
|
140 data.append ('%f %f' % (x, y)) |
|
141 xmin, ymin = min (xmin, x), min (ymin, y) # determine the x range |
|
142 xmax, ymax = max (xmax, x), max (ymax, y) # and the y range |
|
143 data.append ('e') |
|
144 |
|
145 dx, dy = xmax - xmin, ymax - ymin # make x and y ranges equal to the largest |
|
146 if dx > dy: # and keep ranges centered |
|
147 dr = (dx - dy) / 2 |
|
148 ymax += dr |
|
149 ymin -= dr |
|
150 else: |
|
151 dr = (dy - dx) / 2 |
|
152 xmax += dr |
|
153 xmin -= dr |
|
154 |
|
155 from subprocess import Popen, PIPE |
|
156 gnuPlotCmd = ['gnuplot'] |
|
157 plot = Popen (gnuPlotCmd, stdin=PIPE, stdout=PIPE, stderr=PIPE) |
|
158 |
|
159 range = 'set xrange [%f:%f]\nset yrange [%f:%f]\n' % (xmin, xmax, ymin, ymax) |
|
160 plot.stdin.write (range) |
|
161 |
|
162 curves = ','.join ("'-' with linespoints ti '%s: %d punten'" % t for t in npoints) |
|
163 plot.stdin.write ('plot ' + curves + '\n') |
|
164 |
|
165 plot.stdin.write ("\n".join (data)) |
|
166 plot.stdin.write ('\n') |
|
167 plot.stdin.flush () |
|
168 raw_input ('druk') |