You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

1213 lines
40KB

  1. #!/usr/bin/env python
  2. #
  3. '''
  4. Collect and plot latency-profiling data from a running gpsd.
  5. Requires gnuplot, but gnuplot can be on another host.
  6. '''
  7. # This file is Copyright (c) 2010 by the GPSD project
  8. # SPDX-License-Identifier: BSD-2-clause
  9. #
  10. # Updated to conform with RCC-219-00, RCC/IRIG Standard 261-00
  11. # "STANDARD REPORT FORMAT FOR GLOBAL POSITIONING SYSTEM (GPS) RECEIVERS AND
  12. # SYSTEMS ACCURACY TESTS AND EVALUATIONS"
  13. #
  14. # TODO: put date from data on plot, not time of replot.
  15. # TODO: add lat/lon to polar plots
  16. #
  17. # This code runs compatibly under Python 2 and 3.x for x >= 2.
  18. # Preserve this property!
  19. from __future__ import absolute_import, print_function, division
  20. import copy
  21. import getopt
  22. import math
  23. import os
  24. import signal
  25. import socket
  26. import sys
  27. import time
  28. # pylint wants local modules last
  29. try:
  30. import gps
  31. except ImportError as e:
  32. sys.stderr.write(
  33. "gpsprof: can't load Python gps libraries -- check PYTHONPATH.\n")
  34. sys.stderr.write("%s\n" % e)
  35. sys.exit(1)
  36. gps_version = '3.20'
  37. if gps.__version__ != gps_version:
  38. sys.stderr.write("gpsprof: ERROR: need gps module version %s, got %s\n" %
  39. (gps_version, gps.__version__))
  40. sys.exit(1)
  41. debug = False
  42. def dist_2d(a, b):
  43. "calculate distance between a[x,y] and b[x,y]"
  44. # x and y are orthogonal, probably lat/lon in meters
  45. # ignore altitude change.
  46. return math.sqrt((a[0] - b[0]) ** 2 + (a[1] - b[1]) ** 2)
  47. def dist_3d(a, b):
  48. "calculate distance between a[x,y,z] and b[x,y,z]"
  49. # x, y, and z are othogonal, probably ECEF, probably in meters
  50. return math.sqrt((a[0] - b[0]) ** 2 +
  51. (a[1] - b[1]) ** 2 +
  52. (a[2] - b[2]) ** 2)
  53. def wgs84_to_ecef(wgs84):
  54. "Convert wgs84 coordinates to ECEF ones"
  55. # unpack args
  56. (lat, lon, alt) = wgs84
  57. # convert lat/lon/altitude in degrees and altitude in meters
  58. # to ecef x, y, z in meters
  59. # see
  60. # http://www.mathworks.de/help/toolbox/aeroblks/llatoecefposition.html
  61. lat = math.radians(lat)
  62. lon = math.radians(lon)
  63. rad = 6378137.0 # Radius of the Earth (in meters)
  64. f = 1.0 / 298.257223563 # Flattening factor WGS84 Model
  65. cosLat = math.cos(lat)
  66. sinLat = math.sin(lat)
  67. FF = (1.0 - f) ** 2
  68. C = 1 / math.sqrt((cosLat ** 2) + (FF * sinLat ** 2))
  69. S = C * FF
  70. x = (rad * C + alt) * cosLat * math.cos(lon)
  71. y = (rad * C + alt) * cosLat * math.sin(lon)
  72. z = (rad * S + alt) * sinLat
  73. return (x, y, z)
  74. class Baton(object):
  75. "Ship progress indication to stderr."
  76. def __init__(self, prompt, endmsg=None):
  77. self.stream = sys.stderr
  78. self.stream.write(prompt + "...")
  79. if os.isatty(self.stream.fileno()):
  80. self.stream.write(" \b")
  81. self.stream.flush()
  82. self.count = 0
  83. self.endmsg = endmsg
  84. self.time = time.time()
  85. def twirl(self, ch=None):
  86. "Twirl the baton"
  87. if self.stream is None:
  88. return
  89. if ch:
  90. self.stream.write(ch)
  91. elif os.isatty(self.stream.fileno()):
  92. self.stream.write("-/|\\"[self.count % 4])
  93. self.stream.write("\b")
  94. self.count = self.count + 1
  95. self.stream.flush()
  96. return
  97. def end(self, msg=None):
  98. "Write the end message"
  99. if msg is None:
  100. msg = self.endmsg
  101. if self.stream:
  102. self.stream.write("...(%2.2f sec) %s.\n"
  103. % (time.time() - self.time, msg))
  104. class stats(object):
  105. "Class for 1D stats: min, max, mean, sigma, skewness, kurtosis"
  106. def __init__(self):
  107. self.min = 0.0
  108. self.max = 0.0
  109. self.mean = 0.0
  110. self.median = 0.0
  111. self.sigma = 0.0
  112. self.skewness = 0.0
  113. self.kurtosis = 0.0
  114. def __str__(self):
  115. "return a nice string, for debug"
  116. return ("min %f, max %f, mean %f, median %f, sigma %f, skewedness %f, "
  117. "kurtosis %f" %
  118. (self.min, self.max, self.mean, self.median,
  119. self.sigma, self.skewness, self.kurtosis))
  120. def min_max_mean(self, fixes, index):
  121. "Find min, max, and mean of fixes[index]"
  122. if not fixes:
  123. return
  124. # might be fast to go through list once?
  125. if isinstance(fixes[0], tuple):
  126. self.mean = (sum([x[index] for x in fixes]) / len(fixes))
  127. self.min = min([x[index] for x in fixes])
  128. self.max = max([x[index] for x in fixes])
  129. else:
  130. # must be float
  131. self.mean = (sum([x for x in fixes]) / len(fixes))
  132. self.min = min([x for x in fixes])
  133. self.max = max([x for x in fixes])
  134. return
  135. def moments(self, fixes, index):
  136. "Find and set the (sigma, skewness, kurtosis) of fixes[index]"
  137. # The skewness of a random variable X is the third standardized
  138. # moment and is a dimension-less ratio. ntpviz uses the Pearson's
  139. # moment coefficient of skewness. Wikipedia describes it
  140. # best: "The qualitative interpretation of the skew is complicated
  141. # and unintuitive." A normal distribution has a skewness of zero.
  142. self.skewness = float('nan')
  143. # The kurtosis of a random variable X is the fourth standardized
  144. # moment and is a dimension-less ratio. Here we use the Pearson's
  145. # moment coefficient of kurtosis. A normal distribution has a
  146. # kurtosis of three. NIST describes a kurtosis over three as
  147. # "heavy tailed" and one under three as "light tailed".
  148. self.kurtosis = float('nan')
  149. if not fixes:
  150. return
  151. m3 = 0.0
  152. m4 = 0.0
  153. if isinstance(fixes[0], tuple):
  154. sum_squares = [(x[index] - self.mean) ** 2 for x in fixes]
  155. sigma = math.sqrt(sum(sum_squares) / (len(fixes) - 1))
  156. for fix in fixes:
  157. m3 += pow(fix[index] - sigma, 3)
  158. m4 += pow(fix[index] - sigma, 4)
  159. else:
  160. # must be float
  161. sum_squares = [(x - self.mean) ** 2 for x in fixes]
  162. sigma = math.sqrt(sum(sum_squares) / (len(fixes) - 1))
  163. for fix in fixes:
  164. m3 += pow(fix - sigma, 3)
  165. m4 += pow(fix - sigma, 4)
  166. self.sigma = sigma
  167. if sigma > 0.0001:
  168. self.skewness = m3 / (len(fixes) * pow(sigma, 3))
  169. self.kurtosis = m4 / (len(fixes) * pow(sigma, 4))
  170. return
  171. class plotter(object):
  172. "Generic class for gathering and plotting sensor statistics."
  173. def __init__(self):
  174. self.device = None
  175. self.fixes = []
  176. self.in_replot = False
  177. self.session = None
  178. self.start_time = int(time.time())
  179. self.watch = set(['TPV'])
  180. def whatami(self):
  181. "How do we identify this plotting run?"
  182. desc = "%s, %s, " % \
  183. (gps.misc.isotime(self.start_time),
  184. self.device.get('driver', "unknown"))
  185. if 'bps' in self.device:
  186. desc += "%d %dN%d, cycle %.3gs" % \
  187. (self.device['bps'], 9 - self.device['stopbits'],
  188. self.device['stopbits'], self.device['cycle'])
  189. else:
  190. desc += self.device['path']
  191. if 'subtype' in self.device:
  192. desc += "\\n%s" % self.device['subtype']
  193. return desc
  194. def collect(self, verb, log_fp=None):
  195. "Collect data from the GPS."
  196. try:
  197. self.session = gps.gps(host=host, port=port, verbose=verb)
  198. except socket.error:
  199. sys.stderr.write("gpsprof: gpsd unreachable.\n")
  200. sys.exit(1)
  201. # Initialize
  202. self.session.read()
  203. if self.session.version is None:
  204. sys.stderr.write("gpsprof: requires gpsd to speak new protocol.\n")
  205. sys.exit(1)
  206. # Set parameters
  207. flags = gps.WATCH_ENABLE | gps.WATCH_JSON
  208. if self.requires_time:
  209. flags |= gps.WATCH_TIMING
  210. if device:
  211. flags |= gps.WATCH_DEVICE
  212. try:
  213. signal.signal(signal.SIGUSR1,
  214. lambda empty, unused: sys.stderr.write(
  215. "%d of %d (%d%%)..."
  216. % (wait - countdown, wait,
  217. ((wait - countdown) * 100.0 / wait))))
  218. signal.siginterrupt(signal.SIGUSR1, False)
  219. self.session.stream(flags, device)
  220. baton = Baton("gpsprof: %d looking for fix" % os.getpid(), "done")
  221. countdown = wait
  222. basetime = time.time()
  223. while countdown > 0:
  224. if self.session.read() == -1:
  225. sys.stderr.write("gpsprof: gpsd has vanished.\n")
  226. sys.exit(1)
  227. baton.twirl()
  228. if self.session.data["class"] == "ERROR":
  229. sys.stderr.write(" ERROR: %s.\n"
  230. % self.session.data["message"])
  231. sys.exit(1)
  232. if self.session.data["class"] == "DEVICES":
  233. if len(self.session.data["devices"]) != 1 and not device:
  234. sys.stderr.write("ERROR: multiple devices connected, "
  235. "you must explicitly specify the "
  236. "device.\n")
  237. sys.exit(1)
  238. for i in range(len(self.session.data["devices"])):
  239. self.device = copy.copy(
  240. self.session.data["devices"][i])
  241. if self.device['path'] == device:
  242. break
  243. if self.session.data["class"] == "WATCH":
  244. if ((self.requires_time and
  245. not self.session.data.get("timing"))):
  246. sys.stderr.write("timing is not enabled.\n")
  247. sys.exit(1)
  248. # Log before filtering - might be good for post-analysis.
  249. if log_fp:
  250. log_fp.write(self.session.response)
  251. # Ignore everything but what we're told to
  252. if self.session.data["class"] not in self.watch:
  253. continue
  254. # We can get some funky artifacts at start of self.session
  255. # apparently due to RS232 buffering effects. Ignore
  256. # them.
  257. if ((threshold and
  258. time.time() - basetime < self.session.cycle * threshold)):
  259. continue
  260. if self.session.fix.mode <= gps.MODE_NO_FIX:
  261. continue
  262. if self.sample():
  263. if countdown == wait:
  264. sys.stderr.write("first fix in %.2fsec, gathering %d "
  265. "samples..."
  266. % (time.time() - basetime, wait))
  267. countdown -= 1
  268. baton.end()
  269. finally:
  270. self.session.stream(gps.WATCH_DISABLE | gps.WATCH_TIMING)
  271. signal.signal(signal.SIGUSR1, signal.SIG_DFL)
  272. def replot(self, infp):
  273. "Replot from a JSON log file."
  274. self.in_replot = True
  275. baton = Baton("gpsprof: replotting", "done")
  276. self.session = gps.gps(host=None)
  277. for line in infp:
  278. baton.twirl()
  279. self.session.unpack(line)
  280. if self.session.data["class"] == "DEVICES":
  281. self.device = copy.copy(self.session.data["devices"][0])
  282. elif self.session.data["class"] not in self.watch:
  283. continue
  284. self.sample()
  285. baton.end()
  286. def dump(self):
  287. "Dump the raw data for post-analysis."
  288. return self.header() + self.data()
  289. class spaceplot(plotter):
  290. "Spatial scattergram of fixes."
  291. name = "space"
  292. requires_time = False
  293. def __init__(self):
  294. "Initialize class spaceplot"
  295. plotter.__init__(self)
  296. self.centroid = None
  297. self.centroid_ecef = None
  298. self.recentered = []
  299. def sample(self):
  300. "Grab samples"
  301. # Watch out for the NaN value from gps.py.
  302. if (((self.in_replot or self.session.valid) and
  303. self.session.data["class"] == "TPV")):
  304. # get sat used count
  305. sats_used = 0
  306. for sat in self.session.satellites:
  307. if sat.used:
  308. sats_used += 1
  309. if 'altHAE' not in self.session.data:
  310. self.session.data['altHAE'] = gps.NaN
  311. self.fixes.append((self.session.data['lat'],
  312. self.session.data['lon'],
  313. self.session.data['altHAE'], sats_used))
  314. return True
  315. def header(self):
  316. "Return header"
  317. return "\n# Position uncertainty, %s\n" % self.whatami()
  318. def postprocess(self):
  319. "Postprocess the sample data"
  320. return
  321. def data(self):
  322. "Format data for dump"
  323. res = ""
  324. for i in range(len(self.recentered)):
  325. (lat, lon) = self.recentered[i][:2]
  326. (raw1, raw2, alt) = self.fixes[i]
  327. res += "%.9f\t%.9f\t%.9f\t%.9f\t%.9f\n" \
  328. % (lat, lon, raw1, raw2, alt)
  329. return res
  330. def plot(self):
  331. "Plot the data"
  332. stat_lat = stats()
  333. stat_lon = stats()
  334. stat_alt = stats()
  335. stat_used = stats()
  336. # recentered stats
  337. stat_lat_r = stats()
  338. stat_lon_r = stats()
  339. stat_alt_r = stats()
  340. sats_used = []
  341. for x in self.fixes:
  342. # skip missing sats, if any, often missing at start
  343. if x[3] != 0:
  344. sats_used.append(x[3])
  345. # calc sats used data: mean, min, max, sigma
  346. stat_used.min_max_mean(sats_used, 0)
  347. stat_used.moments(sats_used, 0)
  348. # find min, max and mean of lat/lon
  349. stat_lat.min_max_mean(self.fixes, 0)
  350. stat_lon.min_max_mean(self.fixes, 1)
  351. # centroid is just arithmetic avg of lat,lon
  352. self.centroid = (stat_lat.mean, stat_lon.mean)
  353. # Sort fixes by distance from centroid
  354. # sorted to make getting CEP() easy
  355. self.fixes.sort(key=lambda p: dist_2d(self.centroid, p[:2]))
  356. # compute min/max as meters, ignoring altitude
  357. # EarthDistance always returns a positve value
  358. lat_min_o = -gps.EarthDistance((stat_lat.min, self.centroid[1]),
  359. self.centroid[:2])
  360. lat_max_o = gps.EarthDistance((stat_lat.max, self.centroid[1]),
  361. self.centroid[:2])
  362. lon_min_o = -gps.EarthDistance((self.centroid[0], stat_lon.min),
  363. self.centroid[:2])
  364. lon_max_o = gps.EarthDistance((self.centroid[0], stat_lon.max),
  365. self.centroid[:2])
  366. # Convert fixes to offsets from centroid in meters
  367. self.recentered = [
  368. gps.MeterOffset(fix[:2], self.centroid) for fix in self.fixes]
  369. stat_lat_r.min_max_mean(self.recentered, 0)
  370. stat_lon_r.min_max_mean(self.recentered, 1)
  371. # compute sigma, skewness and kurtosis of lat/lon
  372. stat_lat_r.moments(self.recentered, 0)
  373. stat_lon_r.moments(self.recentered, 1)
  374. # CEP(50) calculated per RCC 261-00, Section 3.1.1
  375. calc_cep = 0.5887 * (stat_lat_r.sigma + stat_lon_r.sigma)
  376. # 2DRMS calculated per RCC 261-00, Section 3.1.4
  377. calc_2drms = 2 * math.sqrt(stat_lat_r.sigma ** 2 +
  378. stat_lon_r.sigma ** 2)
  379. # Compute measured CEP(50%)
  380. # same as median distance from centroid, 50% closer, 50% further
  381. cep_meters = gps.misc.EarthDistance(
  382. self.centroid[:2], self.fixes[int(len(self.fixes) * 0.50)][:2])
  383. # Compute measured CEP(95%)
  384. # distance from centroid, 95% closer, 5% further
  385. cep95_meters = gps.misc.EarthDistance(
  386. self.centroid[:2], self.fixes[int(len(self.fixes) * 0.95)][:2])
  387. # Compute measured CEP(99%)
  388. # distance from centroid, 99% closer, 1% further
  389. cep99_meters = gps.misc.EarthDistance(
  390. self.centroid[:2], self.fixes[int(len(self.fixes) * 0.99)][:2])
  391. # Compute CEP(100%)
  392. # max distance from centroid, 100% closer, 0% further
  393. cep100_meters = gps.misc.EarthDistance(
  394. self.centroid[:2], self.fixes[len(self.fixes) - 1][:2])
  395. # init altitude data
  396. alt_ep = gps.NaN
  397. alt_ep95 = gps.NaN
  398. alt_ep99 = gps.NaN
  399. dist_3d_max = 0.0
  400. alt_fixes = []
  401. alt_fixes_r = []
  402. latlon_data = ""
  403. alt_data = ""
  404. # init calcualted hep, sep and mrse
  405. calc_hep = gps.NaN
  406. calc_sep = gps.NaN
  407. calc_mrse = gps.NaN
  408. # grab and format the fixes as gnuplot will use them
  409. for i in range(len(self.recentered)):
  410. # grab valid lat/lon data, recentered and raw
  411. (lat, lon) = self.recentered[i][:2]
  412. alt = self.fixes[i][2]
  413. latlon_data += "%.9f\t%.9f\n" % (lat, lon)
  414. if not math.isnan(alt):
  415. # only keep good fixes
  416. alt_fixes.append(alt)
  417. # micro meters should be good enough
  418. alt_data += "%.6f\n" % (alt)
  419. if alt_fixes:
  420. # got altitude data
  421. # Convert fixes to offsets from avg in meters
  422. alt_data_centered = ""
  423. # find min, max and mean of altitude
  424. stat_alt.min_max_mean(alt_fixes, 0)
  425. for alt in alt_fixes:
  426. alt_fixes_r.append(alt - stat_alt.mean)
  427. alt_data_centered += "%.6f\n" % (alt - stat_alt.mean)
  428. stat_alt_r.min_max_mean(alt_fixes_r, 0)
  429. stat_alt_r.moments(alt_fixes_r, 0)
  430. # centroid in ECEF
  431. self.centroid_ecef = wgs84_to_ecef([stat_lat.mean,
  432. stat_lon.mean,
  433. stat_alt.mean])
  434. # once more through the data, looking for 3D max
  435. for fix_lla in self.fixes:
  436. if not math.isnan(fix_lla[2]):
  437. fix_ecef = wgs84_to_ecef(fix_lla[:3])
  438. dist3d = dist_3d(self.centroid_ecef, fix_ecef)
  439. if dist_3d_max < dist3d:
  440. dist_3d_max = dist3d
  441. # Sort fixes by distance from average altitude
  442. alt_fixes_r.sort(key=lambda a: abs(a))
  443. # so we can rank fixes for EPs
  444. alt_ep = abs(alt_fixes_r[int(len(alt_fixes_r) * 0.50)])
  445. alt_ep95 = abs(alt_fixes_r[int(len(alt_fixes_r) * 0.95)])
  446. alt_ep99 = abs(alt_fixes_r[int(len(alt_fixes_r) * 0.99)])
  447. # HEP(50) calculated per RCC 261-00, Section 3.1.2
  448. calc_hep = 0.6745 * stat_alt_r.sigma
  449. # SEP(50) calculated per RCC 261-00, Section 3.1.3 (3)
  450. calc_sep = 0.51 * (stat_lat_r.sigma +
  451. stat_lon_r.sigma +
  452. stat_alt_r.sigma)
  453. # MRSE calculated per RCC 261-00, Section 3.1.5
  454. calc_mrse = math.sqrt(stat_lat_r.sigma ** 2 +
  455. stat_lon_r.sigma ** 2 +
  456. stat_alt_r.sigma ** 2)
  457. fmt_lab11a = ('hep = %.3f meters\\n'
  458. 'sep = %.3f meters\\n'
  459. 'mrse = %.3f meters\\n'
  460. ) % (calc_hep, calc_sep, calc_mrse)
  461. if self.centroid[0] < 0.0:
  462. latstring = "%.9fS" % -self.centroid[0]
  463. elif stat_lat.mean > 0.0:
  464. latstring = "%.9fN" % self.centroid[0]
  465. else:
  466. latstring = "0.0"
  467. if self.centroid[1] < 0.0:
  468. lonstring = "%.9fW" % -self.centroid[1]
  469. elif stat_lon.mean > 0.0:
  470. lonstring = "%.9fE" % self.centroid[1]
  471. else:
  472. lonstring = "0.0"
  473. # oh, this is fun, mixing gnuplot and python string formatting
  474. # Grrr, python implements %s max width or precision incorrectly...
  475. # and the old and new styles also disagree...
  476. fmt = ('set xlabel "Meters east from %s"\n'
  477. 'set ylabel "Meters north from %s"\n'
  478. 'cep=%.9f\n'
  479. 'cep95=%.9f\n'
  480. 'cep99=%.9f\n'
  481. ) % (lonstring, latstring,
  482. cep_meters, cep95_meters, cep99_meters)
  483. fmt += ('set autoscale\n'
  484. 'set multiplot\n'
  485. # plot to use 95% of width
  486. # set x and y scales to same distance
  487. 'set size ratio -1 0.95,0.7\n'
  488. # leave room at bottom for computed variables
  489. 'set origin 0.025,0.30\n'
  490. 'set format x "%.3f"\n'
  491. 'set format y "%.3f"\n'
  492. 'set key left at screen 0.6,0.30 vertical\n'
  493. 'set key noautotitle\n'
  494. 'set style line 2 pt 1\n'
  495. 'set style line 3 pt 2\n'
  496. 'set style line 5 pt 7 ps 1\n'
  497. 'set xtic rotate by -45\n'
  498. 'set border 15\n'
  499. # now the CEP stuff
  500. 'set parametric\n'
  501. 'set trange [0:2*pi]\n'
  502. 'cx(t, r) = sin(t)*r\n'
  503. 'cy(t, r) = cos(t)*r\n'
  504. 'chlen = cep/20\n'
  505. # what do the next two lines do??
  506. 'set arrow from -chlen,0 to chlen,0 nohead\n'
  507. 'set arrow from 0,-chlen to 0,chlen nohead\n')
  508. fmt += ('set label 11 at screen 0.01, screen 0.30 '
  509. '"RCC 261-00\\n'
  510. 'cep = %.3f meters\\n'
  511. '2drms = %.3f meters\\n%s'
  512. '2d max = %.3f meters\\n'
  513. '3d max = %.3f meters"\n'
  514. ) % (calc_cep, calc_2drms, fmt_lab11a, cep100_meters,
  515. dist_3d_max)
  516. # row labels
  517. fmt += ('set label 12 at screen 0.01, screen 0.12 '
  518. '"RCC 261-00\\n'
  519. '\\n'
  520. 'Lat\\n'
  521. 'Lon\\n'
  522. 'AltHAE\\n'
  523. 'Used"\n')
  524. # mean
  525. fmt += ('set label 13 at screen 0.06, screen 0.12 '
  526. '"\\n'
  527. ' mean\\n'
  528. '%s\\n'
  529. '%s\\n'
  530. '%s\\n'
  531. '%s"\n'
  532. ) % ('{:>15}'.format(latstring),
  533. '{:>15}'.format(lonstring),
  534. '{:>15.3f}'.format(stat_alt.mean),
  535. '{:>15.1f}'.format(stat_used.mean))
  536. fmt += ('set label 14 at screen 0.23, screen 0.12 '
  537. '"\\n'
  538. ' min max sigma '
  539. 'skewness kurtosis\\n'
  540. '%s %s %s meters %s %s\\n'
  541. '%s %s %s meters %s %s\\n'
  542. '%s %s %s meters %s %s\\n'
  543. '%12d %12d %s sats"\n'
  544. ) % ('{:>10.3f}'.format(lat_min_o),
  545. '{:>10.3f}'.format(lat_max_o),
  546. '{:>10.3f}'.format(stat_lat_r.sigma),
  547. '{:>10.1f}'.format(stat_lat_r.skewness),
  548. '{:>10.1f}'.format(stat_lat_r.kurtosis),
  549. '{:>10.3f}'.format(lon_min_o),
  550. '{:>10.3f}'.format(lon_max_o),
  551. '{:>10.3f}'.format(stat_lon_r.sigma),
  552. '{:>10.1f}'.format(stat_lon_r.skewness),
  553. '{:>10.1f}'.format(stat_lon_r.kurtosis),
  554. '{:>10.3f}'.format(stat_alt_r.min),
  555. '{:>10.3f}'.format(stat_alt_r.max),
  556. '{:>10.3f}'.format(stat_alt_r.sigma),
  557. '{:>10.1f}'.format(stat_alt_r.skewness),
  558. '{:>10.1f}'.format(stat_alt_r.kurtosis),
  559. stat_used.min,
  560. stat_used.max,
  561. '{:>10.1f}'.format(stat_used.sigma))
  562. if debug:
  563. fmt += ('set label 15 at screen 0.6, screen 0.12 '
  564. '"\\n'
  565. ' min\\n'
  566. '%s\\n'
  567. '%s\\n'
  568. '%s"\n'
  569. ) % ('{:>15.9f}'.format(stat_lat_r.min),
  570. '{:>15.9f}'.format(stat_lon_r.min),
  571. '{:>15.3f}'.format(stat_alt.min))
  572. fmt += ('set label 16 at screen 0.75, screen 0.12 '
  573. '"\\n'
  574. ' max\\n'
  575. '%s\\n'
  576. '%s\\n'
  577. '%s"\n'
  578. ) % ('{:>15.9f}'.format(stat_lat_r.max),
  579. '{:>15.9f}'.format(stat_lon_r.max),
  580. '{:>15.3f}'.format(stat_alt.max))
  581. if len(self.fixes) > 1000:
  582. plot_style = 'dots'
  583. else:
  584. plot_style = 'points'
  585. # got altitude data?
  586. if not math.isnan(stat_alt.mean):
  587. fmt += ('set ytics nomirror\n'
  588. 'set y2tics\n'
  589. 'set format y2 "%.3f"\n')
  590. fmt += (('set y2label "AltHAE from %.3f meters"\n') %
  591. (stat_alt.mean))
  592. # add ep(50)s
  593. altitude_x = cep100_meters * 1.2
  594. fmt += ('$EPData << EOD\n'
  595. '%.3f %.3f\n'
  596. '%.3f %.3f\n'
  597. 'EOD\n'
  598. ) % (altitude_x, alt_ep,
  599. altitude_x, -alt_ep)
  600. fmt += ('$EP95Data << EOD\n'
  601. '%.3f %.3f\n'
  602. '%.3f %.3f\n'
  603. 'EOD\n'
  604. ) % (altitude_x, alt_ep95,
  605. altitude_x, -alt_ep95)
  606. fmt += ('$EP99Data << EOD\n'
  607. '%.3f %.3f\n'
  608. '%.3f %.3f\n'
  609. 'EOD\n'
  610. ) % (altitude_x, alt_ep99,
  611. altitude_x, -alt_ep99)
  612. # setup now done, plot it!
  613. fmt += ('plot "-" using 1:2 with %s ls 3 title "%d GPS fixes" '
  614. ', cx(t,cep),cy(t,cep) ls 1 title "CEP (50%%) = %.3f meters"'
  615. ', cx(t,cep95),cy(t,cep95) title "CEP (95%%) = %.3f meters"'
  616. ', cx(t,cep99),cy(t,cep99) title "CEP (99%%) = %.3f meters"'
  617. ) % (plot_style, len(self.fixes),
  618. cep_meters, cep95_meters, cep99_meters)
  619. if not math.isnan(stat_alt.mean):
  620. # add plot of altitude
  621. fmt += (', "-" using ( %.3f ):( $1 - %.3f ) '
  622. 'axes x1y2 with points ls 2 lc "green"'
  623. ' title " %d AltHAE fixes"'
  624. ) % (cep100_meters * 1.1, stat_alt.mean, len(alt_fixes))
  625. # altitude EPs
  626. fmt += (', $EPData using 1:2 '
  627. 'axes x1y2 with points ls 5 lc "dark-green"'
  628. ' title " EP(50%%) = %.3f meters"'
  629. ) % (alt_ep)
  630. fmt += (', $EP95Data using 1:2 '
  631. 'axes x1y2 with points ls 5 lc "blue"'
  632. ' title " EP(95%%) = %.3f meters"'
  633. ) % (alt_ep95)
  634. fmt += (', $EP99Data using 1:2 '
  635. 'axes x1y2 with points ls 5 lc "red"'
  636. ' title " EP(99%%) = %.3f meters"'
  637. ) % (alt_ep99)
  638. fmt += self.header() + latlon_data
  639. if not math.isnan(stat_alt.mean):
  640. # add altitude samples
  641. fmt += 'e\n' + alt_data
  642. return fmt
  643. class polarplot(plotter):
  644. "Polar plot of signal strength"
  645. name = "polar"
  646. requires_time = False
  647. seen_used = [] # count of seen and used in each SKY
  648. def __init__(self):
  649. plotter.__init__(self)
  650. self.watch = set(['SKY'])
  651. def sample(self):
  652. "Grab samples"
  653. if self.session.data["class"] == "SKY":
  654. sats = self.session.data['satellites']
  655. seen = 0
  656. used = 0
  657. for sat in sats:
  658. seen += 1
  659. # u'ss': 42, u'el': 15, u'PRN': 18, u'az': 80, u'used': True
  660. if sat['used'] is True:
  661. used += 1
  662. if 'polarunused' == self.name:
  663. continue
  664. if (('polarused' == self.name) and (sat['used'] is False)):
  665. continue
  666. self.fixes.append((sat['PRN'], sat['ss'], sat['az'],
  667. sat['el'], sat['used']))
  668. self.seen_used.append((seen, used))
  669. return True
  670. def header(self):
  671. "Return header"
  672. return "# Polar plot of signal strengths, %s\n" % self.whatami()
  673. def postprocess(self):
  674. "Postprocess the sample data"
  675. return
  676. def data(self):
  677. "Format data for dump"
  678. res = ""
  679. for (prn, ss, az, el, used) in self.fixes:
  680. res += "%d\t%d\t%d\t%d\t%s\n" % (prn, ss, az, el, used)
  681. return res
  682. def plot(self):
  683. "Format data for dump"
  684. # calc SNR: mean, min, max, sigma
  685. stat_ss = stats()
  686. stat_ss.min_max_mean(self.fixes, 1)
  687. stat_ss.moments(self.fixes, 1)
  688. # calc sats seen data: mean, min, max, sigma
  689. stat_seen = stats()
  690. stat_seen.min_max_mean(self.seen_used, 0)
  691. stat_seen.moments(self.seen_used, 0)
  692. # calc sats used data: mean, min, max, sigma
  693. stat_used = stats()
  694. stat_used.min_max_mean(self.seen_used, 1)
  695. stat_used.moments(self.seen_used, 1)
  696. fmt = '''\
  697. unset border
  698. set polar
  699. set angles degrees # set gnuplot on degrees instead of radians
  700. set style line 10 lt 1 lc 0 lw 0.3 #redefine a new line style for the grid
  701. set grid polar 45 #set the grid to be displayed every 45 degrees
  702. set grid ls 10
  703. # x is angle, go from 0 to 360 degrees
  704. # y is radius, go from 90 at middle to 0 at edge
  705. set xrange [0:360]
  706. set rrange [90:0] # 90 at center
  707. set yrange [-90:90]
  708. # set xtics axis #display the xtics on the axis instead of on the border
  709. # set ytics axis
  710. set xtics axis nomirror; set ytics axis nomirror
  711. # "remove" the tics so that only the y tics are displayed
  712. set xtics scale 0
  713. # set the xtics only go from 0 to 90 with increment of 30
  714. # but do not display anything. This has to be done otherwise the grid
  715. # will not be displayed correctly.
  716. set xtics ("" 90, "" 60, "" 30,)
  717. # make the ytics go from the center (0) to 360 with incrment of 90
  718. # set ytics 0, 45, 360
  719. set ytics scale 0
  720. # set the ytics only go from 0 to 90 with increment of 30
  721. # but do not display anything. This has to be done otherwise the grid
  722. # will not be displayed correctly.
  723. set ytics ("" 90, "" 60, "" 30,)
  724. set size square
  725. set key lmargin
  726. # this places a compass label on the outside
  727. set_label(x, text) = sprintf("set label '%s' at (93*cos(%f)), (93*sin(%f)) center", text, x, x)
  728. # here all labels are created
  729. # we compute North (0) at top, East (90) at right
  730. # bug gnuplot puts 0 at right, 90 at top
  731. eval set_label(0, "E")
  732. eval set_label(90, "N")
  733. eval set_label(180, "W")
  734. eval set_label(270, "S")
  735. set style line 11 pt 2 ps 2 #set the line style for the plot
  736. set style fill transparent solid 0.8 noborder
  737. # set rmargin then put colorbox in the margin.
  738. set lmargin at screen .08
  739. set rmargin at screen .85
  740. set cbrange [10:60]
  741. set palette defined (100 "blue", 200 "green", 300 "red")
  742. set colorbox user origin .92, .15 size .03, .6
  743. set label 10 at screen 0.89, screen 0.13 "SNR dBHz"
  744. '''
  745. count = len(self.fixes)
  746. fmt += '''\
  747. set label 11 at screen 0.01, screen 0.15 "%s plot, samples %d"
  748. set label 12 at screen 0.01, screen 0.10 "\\nSS\\nSeen\\nUsed"
  749. ''' % (self.name, count)
  750. fmt += '''\
  751. set label 13 at screen 0.11, screen 0.10 "min\\n%d\\n%d\\n%d" right
  752. ''' % (stat_ss.min, stat_seen.min, stat_used.min)
  753. fmt += '''\
  754. set label 14 at screen 0.21, screen 0.10 "max\\n%d\\n%d\\n%d" right
  755. ''' % (stat_ss.max, stat_seen.max, stat_used.max)
  756. fmt += '''\
  757. set label 15 at screen 0.31, screen 0.10 "mean\\n%.1f\\n%.1f\\n%.1f" right
  758. ''' % (stat_ss.mean, stat_seen.mean, stat_used.mean)
  759. fmt += '''\
  760. set label 16 at screen 0.41, screen 0.10 "sigma\\n%.1f\\n%.1f\\n%.1f" right
  761. ''' % (stat_ss.sigma, stat_seen.sigma, stat_used.sigma)
  762. fmt += '''\
  763. # and finally the plot
  764. # flip azimuth to plot north up, east right
  765. # plot "-" u (90 - $3):4 t "Sat" with points ls 11
  766. plot "-" u (90 - $3):4:(1):($2) t "Sat" w circles lc palette
  767. '''
  768. # return fmt + self.header() + self.data()
  769. return self.header() + fmt + self.data()
  770. class polarplotunused(polarplot):
  771. "Polar plot of unused sats signal strength"
  772. name = "polarunused"
  773. class polarplotused(polarplot):
  774. "Polar plot of used sats signal strength"
  775. name = "polarused"
  776. class timeplot(plotter):
  777. "Time drift against PPS."
  778. name = "time"
  779. requires_time = True
  780. def __init__(self):
  781. plotter.__init__(self)
  782. self.watch = set(['PPS'])
  783. def sample(self):
  784. "Grab samples"
  785. if self.session.data["class"] == "PPS":
  786. self.fixes.append((self.session.data['real_sec'],
  787. self.session.data['real_nsec'],
  788. self.session.data['clock_sec'],
  789. self.session.data['clock_nsec']))
  790. return True
  791. def header(self):
  792. "Return header"
  793. return "# Time drift against PPS, %s\n" % self.whatami()
  794. def postprocess(self):
  795. "Postprocess the sample data"
  796. return
  797. def data(self):
  798. "Format data for dump"
  799. res = ""
  800. for (real_sec, real_nsec, clock_sec, clock_nsec) in self.fixes:
  801. res += "%d\t%d\t%d\t%d\n" % (real_sec, real_nsec, clock_sec,
  802. clock_nsec)
  803. return res
  804. def plot(self):
  805. "Format data for dump"
  806. fmt = '''\
  807. set autoscale
  808. set key below
  809. set ylabel "System clock delta from GPS time (nsec)"
  810. plot "-" using 0:((column(1)-column(3))*1e9 + (column(2)-column(4))) \
  811. title "Delta" with impulses
  812. '''
  813. return fmt + self.header() + self.data()
  814. class uninstrumented(plotter):
  815. "Total times without instrumentation."
  816. name = "uninstrumented"
  817. requires_time = False
  818. def __init__(self):
  819. plotter.__init__(self)
  820. def sample(self):
  821. "Grab samples"
  822. if self.session.fix.time:
  823. seconds = time.time() - gps.misc.isotime(self.session.data.time)
  824. self.fixes.append(seconds)
  825. return True
  826. return False
  827. def header(self):
  828. "Return header"
  829. return "# Uninstrumented total latency, " + self.whatami() + "\n"
  830. def postprocess(self):
  831. "Postprocess the sample data"
  832. return
  833. def data(self):
  834. "Format data for dump"
  835. res = ""
  836. for seconds in self.fixes:
  837. res += "%2.6lf\n" % seconds
  838. return res
  839. def plot(self):
  840. "Plot the data"
  841. fmt = '''\
  842. set autoscale
  843. set key below
  844. set key title "Uninstrumented total latency"
  845. plot "-" using 0:1 title "Total time" with impulses
  846. '''
  847. return fmt + self.header() + self.data()
  848. class instrumented(plotter):
  849. "Latency as analyzed by instrumentation."
  850. name = "instrumented"
  851. requires_time = True
  852. def __init__(self):
  853. "Initialize class instrumented()"
  854. plotter.__init__(self)
  855. def sample(self):
  856. "Grab the samples"
  857. if 'rtime' in self.session.data:
  858. self.fixes.append((gps.misc.isotime(self.session.data['time']),
  859. self.session.data["chars"],
  860. self.session.data['sats'],
  861. self.session.data['sor'],
  862. self.session.data['rtime'],
  863. time.time()))
  864. return True
  865. return False
  866. def header(self):
  867. "Return the header"
  868. res = "# Analyzed latency, " + self.whatami() + "\n"
  869. res += "#-- Fix time -- - Chars - -- Latency - RS232- " \
  870. "Analysis - Recv -\n"
  871. return res
  872. def postprocess(self):
  873. "Postprocess the sample data"
  874. return
  875. def data(self):
  876. "Format data for dump"
  877. res = ""
  878. for (fix_time, chars, sats, start, xmit, recv) in self.fixes:
  879. rs232_time = (chars * 10.0) / self.device['bps']
  880. res += "%.3f %9u %2u %.6f %.6f %.6f %.6f\n" \
  881. % (fix_time, chars, sats, start - fix_time,
  882. (start - fix_time) + rs232_time, xmit - fix_time,
  883. recv - fix_time)
  884. return res
  885. def plot(self):
  886. "Do the plot"
  887. legends = (
  888. "Reception delta",
  889. "Analysis time",
  890. "RS232 time",
  891. "Fix latency",
  892. )
  893. fmt = '''\
  894. set autoscale
  895. set key title "Analyzed latency"
  896. set key below
  897. plot \\\n'''
  898. for (i, legend) in enumerate(legends):
  899. j = len(legends) - i + 4
  900. fmt += ' "-" using 0:%d title "%s" with impulses, \\\n' \
  901. % (j, legend)
  902. fmt = fmt[:-4] + "\n"
  903. return fmt + self.header() + (self.data() + "e\n") * len(legends)
  904. formatters = (polarplot, polarplotunused, polarplotused, spaceplot, timeplot,
  905. uninstrumented, instrumented)
  906. def usage():
  907. "Print help, then exit"
  908. sys.stderr.write('''\
  909. usage: gpsprof [OPTION]... [server[:port[:device]]]
  910. [-D debuglevel]
  911. [-d dumpfile]
  912. [-f {%s}]
  913. [-h]
  914. [-l logfile]
  915. [-m threshold]
  916. [-n samplecount]
  917. [-r]
  918. [-s speed]
  919. [-S subtitle]
  920. [-T terminal]
  921. [-t title]
  922. [-V]
  923. ''' % ("|".join([x.name for x in formatters])))
  924. sys.exit(0)
  925. if __name__ == '__main__':
  926. try:
  927. (options, arguments) = getopt.getopt(sys.argv[1:],
  928. "d:f:hl:m:n:rs:t:D:S:T:V")
  929. plotmode = "space"
  930. raw = False
  931. title = None
  932. subtitle = None
  933. threshold = 0
  934. wait = 100
  935. verbose = 0
  936. terminal = None
  937. dumpfile = None
  938. logfp = None
  939. redo = False
  940. for (switch, val) in options:
  941. if switch == '-d':
  942. dumpfile = val
  943. elif switch == '-D':
  944. # set debug level, 0 off, 1 medium, 2 high
  945. verbose = int(val)
  946. elif switch == '-f':
  947. plotmode = val
  948. elif switch == '-h':
  949. # print help, then exit
  950. usage()
  951. elif switch == '-l':
  952. logfp = open(val, "w")
  953. elif switch == '-m':
  954. threshold = int(val)
  955. elif switch == '-n':
  956. if val[-1] == 'h':
  957. wait = int(val[:-1]) * 360
  958. else:
  959. wait = int(val)
  960. elif switch == '-r':
  961. redo = True
  962. elif switch == '-t':
  963. # replace title
  964. title = val
  965. elif switch == '-S':
  966. # add sub title
  967. subtitle = val
  968. elif switch == '-T':
  969. terminal = val
  970. elif switch == '-V':
  971. sys.stderr.write("gpsprof: Version %s\n" % gps_version)
  972. sys.exit(0)
  973. (host, port, device) = ("localhost", "2947", None)
  974. if arguments:
  975. args = arguments[0].split(":")
  976. if len(args) >= 1:
  977. host = args[0]
  978. if len(args) >= 2:
  979. port = args[1]
  980. if len(args) >= 3:
  981. device = args[2]
  982. # Select the plotting mode
  983. if plotmode:
  984. for formatter in formatters:
  985. if formatter.name == plotmode:
  986. plot = formatter()
  987. break
  988. else:
  989. sys.stderr.write("gpsprof: no such formatter.\n")
  990. sys.exit(1)
  991. # Get fix data from the GPS
  992. if redo:
  993. plot.replot(sys.stdin)
  994. else:
  995. plot.collect(verbose, logfp)
  996. plot.postprocess()
  997. # Save the timing data (only) for post-analysis if required.
  998. if dumpfile:
  999. with open(dumpfile, "w") as fp:
  1000. fp.write(plot.dump())
  1001. if logfp:
  1002. logfp.close()
  1003. # Ship the plot to standard output
  1004. if not title:
  1005. title = plot.whatami()
  1006. # escape " for gnuplot
  1007. title = title.replace('"', '\\"')
  1008. if subtitle:
  1009. title += '\\n' + subtitle
  1010. if terminal:
  1011. sys.stdout.write("set terminal %s truecolor enhanced size"
  1012. " 800,950\n"
  1013. % terminal)
  1014. # double quotes on title so \n is parsed by gnuplot
  1015. sys.stdout.write('set title noenhanced "%s\\n\\n"\n' % title)
  1016. sys.stdout.write(plot.plot())
  1017. except getopt.GetoptError as error:
  1018. print(error)
  1019. print("")
  1020. usage()
  1021. exit(1)
  1022. except KeyboardInterrupt:
  1023. pass
  1024. # The following sets edit modes for GNU EMACS
  1025. # Local Variables:
  1026. # mode:python
  1027. # End: