From 29cdbc7199b53dee0c36873fc288acf49dc4db58 Mon Sep 17 00:00:00 2001 From: "Gary E. Miller" Date: Mon, 25 Jun 2018 20:55:47 -0700 Subject: gpsprof: Add statistics to scatterplot. Add satellite heat maps. Plus the man page. This work supported by the kind generatosity of Virgin Orbit. --- gpsprof | 784 +++++++++++++++++++++++++++++++++++++++++++++++++++--------- gpsprof.xml | 201 +++++++++++----- 2 files changed, 804 insertions(+), 181 deletions(-) diff --git a/gpsprof b/gpsprof index fecbd47a..fb01fecb 100755 --- a/gpsprof +++ b/gpsprof @@ -6,6 +6,10 @@ # Collect and plot latency-profiling data from a running gpsd. # Requires gnuplot. # +# Updated to conform with RCC-219-00, RCC/IRIG Standard 261-00 +# "STANDARD REPORT FORMAT FOR GLOBAL POSITIONING SYSTEM (GPS) RECEIVERS AND +# SYSTEMS ACCURACY TESTS AND EVALUATIONS" +# # This code runs compatibly under Python 2 and 3.x for x >= 2. # Preserve this property! from __future__ import absolute_import, print_function, division @@ -21,6 +25,50 @@ import time import gps +debug = False + + +def dist_2d(a, b): + # calculate distance between a[x,y] and b[x,y] + # x and y are orthogonal, probably lat/lon in meters + # ignore altitude change. + return math.sqrt((a[0] - b[0]) ** 2 + (a[1] - b[1]) ** 2) + + +def dist_3d(a, b): + # calculate distance between a[x,y,z] and b[x,y,z] + # x, y, and z are othogonal, probably ECEF, probably in meters + return math.sqrt((a[0] - b[0]) ** 2 + + (a[1] - b[1]) ** 2 + + (a[2] - b[2]) ** 2) + + +def wgs84_to_ecef(wgs84): + "Convert wgs84 coordinates to ECEF ones" + + # unpack args + (lat, lon, alt) = wgs84 + # convert lat/lon/altitude in degrees and altitude in meters + # to ecef x, y, z in meters + # see + # http://www.mathworks.de/help/toolbox/aeroblks/llatoecefposition.html + lat = math.radians(lat) + lon = math.radians(lon) + + rad = 6378137.0 # Radius of the Earth (in meters) + f = 1.0 / 298.257223563 # Flattening factor WGS84 Model + cosLat = math.cos(lat) + sinLat = math.sin(lat) + FF = (1.0 - f) ** 2 + C = 1 / math.sqrt((cosLat ** 2) + (FF * sinLat ** 2)) + S = C * FF + + x = (rad * C + alt) * cosLat * math.cos(lon) + y = (rad * C + alt) * cosLat * math.sin(lon) + z = (rad * S + alt) * sinLat + + return (x, y, z) + class Baton(object): "Ship progress indication to stderr." @@ -57,11 +105,87 @@ class Baton(object): return +class stats(object): + "Class for 1D stats: min, max, mean, sigma, skewness, kurtosis" + + def __init__(self): + self.min = 0.0 + self.max = 0.0 + self.mean = 0.0 + self.median = 0.0 + self.sigma = 0.0 + self.skewness = 0.0 + self.kurtosis = 0.0 + + def min_max_mean(self, fixes, index): + "Find min, max, and mean of fixes[index]" + + if 0 == len(fixes): + return + + total = 0.0 + # might be fast to go through list once? + if type(fixes[0]) == tuple: + self.mean = (sum([x[index] for x in fixes]) / len(fixes)) + self.min = min([x[index] for x in fixes]) + self.max = max([x[index] for x in fixes]) + else: + # must be float + self.mean = (sum([x for x in fixes]) / len(fixes)) + self.min = min([x for x in fixes]) + self.max = max([x for x in fixes]) + + return + + def moments(self, fixes, index): + "Find and set the (sigma, skewness, kurtosis) of fixes[index]" + + # The skewness of a random variable X is the third standardized + # moment and is a dimension-less ratio. ntpviz uses the Pearson's + # moment coefficient of skewness. Wikipedia describes it + # best: "The qualitative interpretation of the skew is complicated + # and unintuitive." A normal distribution has a skewness of zero. + skewness = float('nan') + + # The kurtosis of a random variable X is the fourth standardized + # moment and is a dimension-less ratio. Here we use the Pearson's + # moment coefficient of kurtosis. A normal distribution has a + # kurtosis of three. NIST describes a kurtosis over three as + # "heavy tailed" and one under three as "light tailed". + kurtosis = float('nan') + + if 0 == len(fixes): + return + + m3 = 0.0 + m4 = 0.0 + if type(fixes[0]) == tuple: + sum_squares = [x[index] ** 2 for x in fixes] + sigma = math.sqrt(sum(sum_squares) / len(fixes)) + for fix in fixes: + m3 += pow(fix[index] - sigma, 3) + m4 += pow(fix[index] - sigma, 4) + else: + # must be float + sum_squares = [x ** 2 for x in fixes] + sigma = math.sqrt(sum(sum_squares) / len(fixes)) + for fix in fixes: + m3 += pow(fix - sigma, 3) + m4 += pow(fix - sigma, 4) + + self.sigma = sigma + self.skewness = m3 / (len(fixes) * pow(sigma, 3)) + self.kurtosis = m4 / (len(fixes) * pow(sigma, 4)) + + return + + class plotter(object): "Generic class for gathering and plotting sensor statistics." def __init__(self): self.fixes = [] + self.in_replot = False self.start_time = int(time.time()) self.watch = set(['TPV']) @@ -71,11 +195,14 @@ class plotter(object): (gps.misc.isotime(self.start_time), self.device.get('driver', "unknown")) if 'bps' in self.device: - desc += "%d %dN%d, cycle %ds" % \ + desc += "%d %dN%d, cycle %.3gs" % \ (self.device['bps'], 9 - self.device['stopbits'], self.device['stopbits'], self.device['cycle']) else: desc += self.device['path'] + if 'subtype' in self.device: + desc += "\\n%s" % self.device['subtype'] + return desc def collect(self, verbose, logfp=None): @@ -128,8 +255,8 @@ class plotter(object): if self.device['path'] == device: break if self.session.data["class"] == "WATCH": - if ((self.requires_time - and not self.session.data.get("timing"))): + if ((self.requires_time and + not self.session.data.get("timing"))): sys.stderr.write("timing is not enabled.\n") sys.exit(1) # Log before filtering - might be good for post-analysis. @@ -159,6 +286,7 @@ class plotter(object): def replot(self, infp): "Replot from a JSON log file." + self.in_replot = True baton = Baton("gpsprof: replotting", "done") self.session = gps.gps(host=None) for line in infp: @@ -185,33 +313,29 @@ class spaceplot(plotter): plotter.__init__(self) self.recentered = [] - def d(self, a, b): - return math.sqrt((a[0] - b[0]) ** 2 + (a[1] - b[1]) ** 2) - def sample(self): # Watch out for the NaN value from gps.py. - if (self.session.valid & - (gps.ONLINE_SET | gps.PACKET_SET | gps.LATLON_SET)): - self.fixes.append((self.session.fix.latitude, - self.session.fix.longitude, - self.session.fix.altitude)) + if (((self.in_replot or self.session.valid) and + self.session.data["class"] == "TPV")): + # get sat used count + sats_used = 0 + for sat in self.session.satellites: + if sat.used: + sats_used += 1 + + self.fixes.append((self.session.data['lat'], + self.session.data['lon'], + self.session.data['alt'], sats_used)) return True def header(self): - return "# Position uncertainty, %s\n" % self.whatami() + return "\n# Position uncertainty, %s\n" % self.whatami() def postprocess(self): - if not self.recentered: - # centroid is just arithmetic avg of lat,lon - self.centroid = (sum([x[0] for x in self.fixes]) / len(self.fixes), - sum([x[1] for x in self.fixes]) / len(self.fixes)) - # Sort fixes by distance from centroid - self.fixes.sort(key=lambda p: self.d(self.centroid, p)) - # Convert fixes to offsets from centroid in meters - self.recentered = [ - gps.MeterOffset(self.centroid, fix[:2]) for fix in self.fixes] + pass def data(self): + # dump of data res = "" for i in range(len(self.recentered)): (lat, lon) = self.recentered[i][:2] @@ -221,102 +345,492 @@ class spaceplot(plotter): return res def plot(self): - # Compute CEP(50%) + stat_lat = stats() + stat_lon = stats() + stat_alt = stats() + stat_used = stats() + + sats_used = [] + for x in self.fixes: + # skip missing sats, if any, often missing at start + if x[3] != 0: + sats_used.append(x[3]) + + # calc sats used data: mean, min, max, sigma + stat_used.min_max_mean(sats_used, 0) + stat_lat.moments(sats_used, 0) + + # find min, max and mean of lat/lon + stat_lat.min_max_mean(self.fixes, 0) + stat_lon.min_max_mean(self.fixes, 1) + + # centroid is just arithmetic avg of lat,lon + self.centroid = (stat_lat.mean, stat_lon.mean) + + # Sort fixes by distance from centroid + # sorted to make getting CEP() easy + self.fixes.sort(key=lambda p: dist_2d(self.centroid, p[:2])) + + # compute min/max as meters, ignoring altitude + # EarthDistance always returns a positve value + lat_min_o = -gps.EarthDistance((stat_lat.min, self.centroid[1]), + self.centroid[:2]) + lat_max_o = gps.EarthDistance((stat_lat.max, self.centroid[1]), + self.centroid[:2]) + + lon_min_o = -gps.EarthDistance((self.centroid[0], stat_lon.min), + self.centroid[:2]) + lon_max_o = gps.EarthDistance((self.centroid[0], stat_lon.max), + self.centroid[:2]) + + # Convert fixes to offsets from centroid in meters + self.recentered = [ + gps.MeterOffset(fix[:2], self.centroid) for fix in self.fixes] + + # compute sigma, skewness and kurtosis of lat/lon + stat_lat.moments(self.recentered, 0) + stat_lon.moments(self.recentered, 1) + + # CEP(50) calculated per RCC 261-00, Section 3.1.1 + calc_cep = 0.5887 * (stat_lat.sigma + stat_lon.sigma) + + # 2DRMS calculated per RCC 261-00, Section 3.1.4 + calc_2drms = 2 * math.sqrt(stat_lat.sigma ** 2 + stat_lon.sigma ** 2) + + # Compute measured CEP(50%) + # same as median distance from centroid, 50% closer, 50% further cep_meters = gps.misc.EarthDistance( self.centroid[:2], self.fixes[int(len(self.fixes) * 0.50)][:2]) + + # Compute measured CEP(95%) + # distance from centroid, 95% closer, 5% further cep95_meters = gps.misc.EarthDistance( self.centroid[:2], self.fixes[int(len(self.fixes) * 0.95)][:2]) + + # Compute measured CEP(99%) + # distance from centroid, 99% closer, 1% further cep99_meters = gps.misc.EarthDistance( self.centroid[:2], self.fixes[int(len(self.fixes) * 0.99)][:2]) - alt_sum = 0.0 - alt_num = 0 + + # Compute CEP(100%) + # max distance from centroid, 100% closer, 0% further + cep100_meters = gps.misc.EarthDistance( + self.centroid[:2], self.fixes[len(self.fixes) - 1][:2]) + + # init altitude data + alt_ep = gps.NaN + alt_ep95 = gps.NaN + alt_ep99 = gps.NaN + fmt_lan11 = '' + dist_3d_max = 0.0 alt_fixes = [] - lon_max = -9999 + latlon_data = "" + alt_data = "" + + # grab and format the fixes as gnuplot will use them for i in range(len(self.recentered)): - (_lat, lon) = self.recentered[i][:2] - (_raw1, _raw2, alt) = self.fixes[i] - if not gps.isnan(alt): - alt_sum += alt + # grab valid lat/lon data, recentered and raw + (lat, lon) = self.recentered[i][:2] + alt = self.fixes[i][2] + latlon_data += "%.9f\t%.9f\n" % (lat, lon) + + if not math.isnan(alt): + # only keep good fixes alt_fixes.append(alt) - alt_num += 1 - if lon > lon_max: - lon_max = lon - if alt_num == 0: - alt_avg = gps.NaN - alt_ep = gps.NaN - else: - alt_avg = alt_sum / alt_num + # micro meters should be good enough + alt_data += "%.6f\n" % (alt) + + if 0 < len(alt_fixes): + # got altitude data + + # find min, max and mean of altitude + stat_alt.min_max_mean(alt_fixes, 0) + + # centroid in ECEF + self.centroid_ecef = wgs84_to_ecef([stat_lat.mean, + stat_lon.mean, + stat_alt.mean]) + + # once more through the data, looking for 3D max + for i in range(len(self.fixes)): + fix_lla = self.fixes[i][:3] + if not math.isnan(fix_lla[2]): + fix_ecef = wgs84_to_ecef(fix_lla[:3]) + dist3d = dist_3d(self.centroid_ecef, fix_ecef) + if dist_3d_max < dist3d: + dist_3d_max = dist3d + + # Convert fixes to offsets from avg in meters + alt_centered = [] + alt_data_centered = "" + alt_fixes.sort(key=lambda a: abs(a)) + for i in alt_fixes: + alt_centered.append(i - stat_alt.mean) + alt_data_centered += "%.6f\n" % (i - stat_alt.mean) + # Sort fixes by distance from average altitude - alt_fixes.sort(key=lambda a: abs(alt_avg - a)) - alt_ep = abs(alt_fixes[len(alt_fixes) // 2] - alt_avg) - if self.centroid[0] < 0: - latstring = "%.9fS" % -self.centroid[0] - elif self.centroid[0] == 0: - latstring = "0" + alt_centered.sort(key=lambda a: abs(a)) + # so we can rank fixes for EPs + alt_ep = abs(alt_centered[int(len(alt_centered) * 0.50)]) + alt_ep95 = abs(alt_centered[int(len(alt_centered) * 0.95)]) + alt_ep99 = abs(alt_centered[int(len(alt_centered) * 0.99)]) + + stat_alt.moments(alt_centered, 0) + + # HEP(50) calculated per RCC 261-00, Section 3.1.2 + calc_hep = 0.6745 * stat_alt.sigma + + # SEP(50) calculated per RCC 261-00, Section 3.1.3 (3) + calc_sep = 0.51 * (stat_lat.sigma + + stat_lon.sigma + + stat_alt.sigma) + + # MRSE calculated per RCC 261-00, Section 3.1.5 + calc_mrse = math.sqrt(stat_lat.sigma ** 2 + + stat_lon.sigma ** 2 + + stat_alt.sigma ** 2) + + fmt_lab11a = ('hep = %.3f meters\\n' + 'sep = %.3f meters\\n' + 'mrse = %.3f meters\\n' + ) % (calc_hep, calc_sep, calc_mrse) + + if stat_lat.mean < 0.0: + latstring = "%.9fS" % -stat_lat.mean + elif stat_lat.mean > 0.0: + latstring = "%.9fN" % stat_lat.mean else: - latstring = "%.9fN" % self.centroid[0] - if self.centroid[1] < 0: - lonstring = "%.9fW" % -self.centroid[1] - elif self.centroid[1] == 0: - lonstring = "0" + latstring = "0.0" + + if stat_lon.mean < 0.0: + lonstring = "%.9fW" % -stat_lon.mean + elif stat_lon.mean > 0.0: + lonstring = "%.9fE" % stat_lon.mean else: - lonstring = "%.9fE" % self.centroid[1] - fmt = "set autoscale\n" - fmt += "set format x \"%.3f\"\n" - fmt += "set format y \"%.3f\"\n" - fmt += 'set key below\n' - fmt += 'set key title "%s"\n' % gps.misc.isotime(int(time.time())) - fmt += 'set size ratio -1\n' - fmt += 'set style line 2 pt 1\n' - fmt += 'set style line 3 pt 2\n' - fmt += 'set xlabel "Meters east from %s"\n' % lonstring - fmt += 'set xtic rotate by -45\n' - fmt += 'set ylabel "Meters north from %s"\n' % latstring - fmt += 'set border 15\n' - if not gps.isnan(alt_avg): - fmt += 'set y2label "Meters Altitude from %.3f"\n' % alt_avg - fmt += 'set ytics nomirror\n' - fmt += 'set y2tics\n' - fmt += "set format y2 \"%.3f\"\n" - fmt += 'cep=%.9f\n' \ - % self.d((0, 0), self.recentered[len(self.fixes) // 2]) - fmt += 'cep95=%.9f\n' \ - % self.d((0, 0), self.recentered[int(len(self.fixes) * 0.95)]) - fmt += 'cep99=%.9f\n' \ - % self.d((0, 0), self.recentered[int(len(self.fixes) * 0.99)]) - fmt += 'set parametric\n' - fmt += 'set trange [0:2*pi]\n' - fmt += 'cx(t, r) = sin(t)*r\n' - fmt += 'cy(t, r) = cos(t)*r\n' - fmt += 'chlen = cep/20\n' - fmt += "set arrow from -chlen,0 to chlen,0 nohead\n" - fmt += "set arrow from 0,-chlen to 0,chlen nohead\n" + lonstring = "0.0" + + # oh, this is fun, mixing gnuplot and python string formatting + # Grrr, python implements %s max width or precision incorrectly... + # and the old and new styles also disagree... + fmt = ('set xlabel "Meters east from %s"\n' + 'set ylabel "Meters north from %s"\n' + 'cep=%.9f\n' + 'cep95=%.9f\n' + 'cep99=%.9f\n' + ) % (lonstring, latstring, + cep_meters, cep95_meters, cep99_meters) + fmt += ('set autoscale\n' + 'set multiplot\n' + # plot to use 95% of width + # set x and y scales to same distance + 'set size ratio -1 0.95,0.7\n' + # leave room at bottom for computed variables + 'set origin 0.025,0.30\n' + 'set format x "%.3f"\n' + 'set format y "%.3f"\n' + 'set key left at screen 0.6,0.30 vertical\n' + 'set key noautotitle\n' + 'set style line 2 pt 1\n' + 'set style line 3 pt 2\n' + 'set style line 5 pt 7 ps 1\n' + 'set xtic rotate by -45\n' + 'set border 15\n' + # now the CEP stuff + 'set parametric\n' + 'set trange [0:2*pi]\n' + 'cx(t, r) = sin(t)*r\n' + 'cy(t, r) = cos(t)*r\n' + 'chlen = cep/20\n' + # what do the next two lines do?? + 'set arrow from -chlen,0 to chlen,0 nohead\n' + 'set arrow from 0,-chlen to 0,chlen nohead\n') + + fmt += ('set label 11 at screen 0.01, screen 0.30 ' + '"RCC 261-00\\n' + 'cep = %.3f meters\\n' + '2drms = %.3f meters\\n%s' + '2d max = %.3f meters\\n' + '3d max = %.3f meters"\n' + ) % (calc_cep, calc_2drms, fmt_lab11a, cep100_meters, + dist_3d_max) + + # row labels + fmt += ('set label 12 at screen 0.01, screen 0.12 ' + '"RCC 261-00\\n' + '\\n' + 'Lat\\n' + 'Lon\\n' + 'Alt\\n' + 'Used"\n') + + # mean + fmt += ('set label 13 at screen 0.06, screen 0.12 ' + '"\\n' + ' mean\\n' + '%s\\n' + '%s\\n' + '%s\\n' + '%s"\n' + ) % ('{:>15}'.format(latstring), + '{:>15}'.format(lonstring), + '{:>15.3f}'.format(stat_alt.mean), + '{:>15.1f}'.format(stat_used.mean)) + + fmt += ('set label 14 at screen 0.23, screen 0.12 ' + '"\\n' + ' min max sigma ' + 'skewness kurtosis\\n' + '%s %s %s meters %s %s\\n' + '%s %s %s meters %s %s\\n' + '%s %s %s meters %s %s\\n' + '%12d %12d %s sats"\n' + ) % ('{:>10.3f}'.format(lat_min_o), + '{:>10.3f}'.format(lat_max_o), + '{:>10.3f}'.format(stat_lat.sigma), + '{:>10.1f}'.format(stat_lat.skewness), + '{:>10.1f}'.format(stat_lat.kurtosis), + '{:>10.3f}'.format(lon_min_o), + '{:>10.3f}'.format(lon_max_o), + '{:>10.3f}'.format(stat_lon.sigma), + '{:>10.1f}'.format(stat_lon.skewness), + '{:>10.1f}'.format(stat_lon.kurtosis), + '{:>10.3f}'.format(stat_alt.min - stat_alt.mean), + '{:>10.3f}'.format(stat_alt.max - stat_alt.mean), + '{:>10.3f}'.format(stat_alt.sigma), + '{:>10.1f}'.format(stat_alt.skewness), + '{:>10.1f}'.format(stat_alt.kurtosis), + stat_used.min, + stat_used.max, + '{:>10.1f}'.format(stat_used.sigma)) + + if debug: + fmt += ('set label 15 at screen 0.6, screen 0.12 ' + '"\\n' + ' min\\n' + '%s\\n' + '%s\\n' + '%s"\n' + ) % ('{:>15.9f}'.format(stat_lat.min), + '{:>15.9f}'.format(stat_lon.min), + '{:>15.3f}'.format(stat_alt.min)) + + fmt += ('set label 16 at screen 0.75, screen 0.12 ' + '"\\n' + ' max\\n' + '%s\\n' + '%s\\n' + '%s"\n' + ) % ('{:>15.9f}'.format(stat_lat.max), + '{:>15.9f}'.format(stat_lon.max), + '{:>15.3f}'.format(stat_alt.max)) + if len(self.fixes) > 1000: plot_style = 'dots' else: plot_style = 'points' - fmt += 'plot "-" using 1:2 with ' + plot_style \ - + ' ls 3 title "%d GPS fixes" ' % (len(self.fixes)) - if not gps.isnan(alt_avg): - fmt += ', "-" using ( %.3f ):($5 < 100000 ? $5 - %.3f : 1/0)' \ - 'axes x1y2 with %s ls 2 title ' \ - '" %d Altitude fixes, Average = %.3f, EP (50%%) = %.3f"' \ - % (lon_max + 1, alt_avg, plot_style, alt_num, alt_avg, - alt_ep) - fmt += ', cx(t, cep),cy(t, cep) ls 1 title "CEP (50%%) = %.3f meters"' \ - % (cep_meters) - fmt += ', cx(t, cep95),cy(t, cep95) title "CEP (95%%) = %.3f meters"' \ - % (cep95_meters) - fmt += ', cx(t, cep99),cy(t, cep99) title "CEP (99%%) = %.3f meters"' \ - % (cep99_meters) - fmt += "\n" - fmt += self.header() - fmt += self.data() - if not gps.isnan(alt_avg): - fmt += "e\n" + self.data() + + # got altitude data? + if not math.isnan(stat_alt.mean): + fmt += ('set ytics nomirror\n' + 'set y2tics\n' + 'set format y2 "%.3f"\n') + fmt += (('set y2label "Altitude from %.3f meters"\n') % + (stat_alt.mean)) + + # add ep(50)s + + altitude_x = cep100_meters * 1.2 + fmt += ('$EPData << EOD\n' + '%.3f %.3f\n' + '%.3f %.3f\n' + 'EOD\n' + ) % (altitude_x, alt_ep, + altitude_x, -alt_ep) + + fmt += ('$EP95Data << EOD\n' + '%.3f %.3f\n' + '%.3f %.3f\n' + 'EOD\n' + ) % (altitude_x, alt_ep95, + altitude_x, -alt_ep95) + + fmt += ('$EP99Data << EOD\n' + '%.3f %.3f\n' + '%.3f %.3f\n' + 'EOD\n' + ) % (altitude_x, alt_ep99, + altitude_x, -alt_ep99) + + # setup now done, plot it! + fmt += ('plot "-" using 1:2 with %s ls 3 title "%d GPS fixes" ' + ', cx(t,cep),cy(t,cep) ls 1 title "CEP (50%%) = %.3f meters"' + ', cx(t,cep95),cy(t,cep95) title "CEP (95%%) = %.3f meters"' + ', cx(t,cep99),cy(t,cep99) title "CEP (99%%) = %.3f meters"' + ) % (plot_style, len(self.fixes), + cep_meters, cep95_meters, cep99_meters) + + if not math.isnan(stat_alt.mean): + # add plot of altitude + fmt += (', "-" using ( %.3f ):( $1 - %.3f ) ' + 'axes x1y2 with points ls 2 lc "green"' + ' title " %d Altitude fixes"' + ) % (cep100_meters * 1.1, stat_alt.mean, len(alt_fixes)) + + # altitude EPs + fmt += (', $EPData using 1:2 ' + 'axes x1y2 with points ls 5 lc "dark-green"' + ' title " EP(50%%) = %.3f meters"' + ) % (alt_ep) + + fmt += (', $EP95Data using 1:2 ' + 'axes x1y2 with points ls 5 lc "blue"' + ' title " EP(95%%) = %.3f meters"' + ) % (alt_ep95) + + fmt += (', $EP99Data using 1:2 ' + 'axes x1y2 with points ls 5 lc "red"' + ' title " EP(99%%) = %.3f meters"' + ) % (alt_ep99) + + fmt += self.header() + latlon_data + if not math.isnan(stat_alt.mean): + # add altitude samples + fmt += 'e\n' + alt_data return fmt +class polarplot(plotter): + "Polar plot of signal strength" + name = "polar" + requires_time = False + + def __init__(self): + plotter.__init__(self) + self.watch = set(['SKY']) + + def sample(self): + if self.session.data["class"] == "SKY": + sats = self.session.data['satellites'] + seen = 0 + used = 0 + for sat in sats: + # u'ss': 42, u'el': 15, u'PRN': 18, u'az': 80, u'used': True + if (('polarunused' == self.name) and (sat['used'] is True)): + continue + if (('polarused' == self.name) and (sat['used'] is False)): + continue + self.fixes.append((sat['PRN'], sat['ss'], sat['az'], + sat['el'], sat['used'])) + + return True + + def header(self): + return "# Polar plot of signal strengths, %s\n" % self.whatami() + + def postprocess(self): + pass + + def data(self): + res = "" + for (prn, ss, az, el, used) in self.fixes: + res += "%d\t%d\t%d\t%d\t%s\n" % (prn, ss, az, el, used) + return res + + def plot(self): + ss_min = 99 # min satellite signal strength + ss_max = 0 # max satellite signal strength + for fix in self.fixes: + # collect min and max + if ss_max < fix[1]: + ss_max = fix[1] + if 0 < fix[1] and ss_min > fix[1]: + ss_min = fix[1] + + fmt = '''\ +unset border +set polar +set angles degrees # set gnuplot on degrees instead of radians + +set style line 10 lt 1 lc 0 lw 0.3 #redefine a new line style for the grid + +set grid polar 45 #set the grid to be displayed every 45 degrees +set grid ls 10 + +# x is angle, go from 0 to 360 degrees +# y is radius, go from 90 at middle to 0 at edge +set xrange [0:360] +set rrange [90:0] # 90 at center +set yrange [-90:90] + +# set xtics axis #display the xtics on the axis instead of on the border +# set ytics axis +set xtics axis nomirror; set ytics axis nomirror + +# "remove" the tics so that only the y tics are displayed +set xtics scale 0 +# set the xtics only go from 0 to 90 with increment of 30 +# but do not display anything. This has to be done otherwise the grid +# will not be displayed correctly. +set xtics ("" 90, "" 60, "" 30,) + +# make the ytics go from the center (0) to 360 with incrment of 90 +# set ytics 0, 45, 360 +set ytics scale 0 + +# set the ytics only go from 0 to 90 with increment of 30 +# but do not display anything. This has to be done otherwise the grid +# will not be displayed correctly. +set ytics ("" 90, "" 60, "" 30,) + +set size square + +set key lmargin + +# this places a compass label on the outside +set_label(x, text) = sprintf("set label '%s' at (93*cos(%f)), (93*sin(%f)) center", text, x, x) + +# here all labels are created +# we compute North (0) at top, East (90) at right +# bug gnuplot puts 0 at right, 90 at top +eval set_label(0, "E") +eval set_label(90, "N") +eval set_label(180, "W") +eval set_label(270, "S") + +set style line 11 pt 2 ps 2 #set the line style for the plot + +set style fill transparent solid 0.8 noborder + +set cbrange [10:60] +set palette defined (100 "blue", 200 "green", 300 "red") +''' + + count = len(self.fixes) + fmt += '''\ +set label 11 at screen 0.01, screen 0.15 "%s\\nmin SS = %d\\nmax SS = %d\\n\ +Samples %d" +''' % (self.name, ss_min, ss_max, count) + + fmt += '''\ +# and finally the plot +# flip azimuth to plot north up, east right +# plot "-" u (90 - $3):4 t "Sat" with points ls 11 +plot "-" u (90 - $3):4:(1):($2) t "Sat" w circles lc palette +''' + # return fmt + self.header() + self.data() + return self.header() + fmt + self.data() + + +class polarplotunused(polarplot): + "Polar plot of unused sats signal strength" + name = "polarunused" + + +class polarplotused(polarplot): + "Polar plot of used sats signal strength" + name = "polarused" + + class timeplot(plotter): "Time drift against PPS." name = "time" @@ -454,16 +968,40 @@ plot \\\n''' fmt = fmt[:-4] + "\n" return fmt + self.header() + (self.data() + "e\n") * len(legends) -formatters = (spaceplot, timeplot, uninstrumented, instrumented) + +formatters = (polarplot, polarplotunused, polarplotused, spaceplot, timeplot, + uninstrumented, instrumented) + + +def usage(): + "Print help, then exit" + sys.stderr.write('''\ +usage: gpsprof [OPTION]... [server[:port[:device]]] + [-D debuglevel] + [-d dumpfile] + [-f {%s}] + [-h] + [-l logfile] + [-m threshold] + [-n samplecount] + [-r] + [-s speed] + [-S subtitle] + [-T terminal] + [-t title] + +''' % ("|".join([x.name for x in formatters]))) + sys.exit(0) if __name__ == '__main__': try: (options, arguments) = getopt.getopt(sys.argv[1:], - "d:f:hl:m:n:rs:t:T:D:") + "d:f:hl:m:n:rs:S:t:T:D:") plotmode = "space" raw = False title = None + subtitle = None threshold = 0 await = 100 verbose = 0 @@ -472,8 +1010,19 @@ if __name__ == '__main__': logfp = None redo = False for (switch, val) in options: - if switch == '-f': + if switch == '-d': + dumpfile = val + elif switch == '-D': + # set debug level, 0 off, 1 medium, 2 high + verbose = int(val) + elif switch == '-f': plotmode = val + elif switch == '-h': + # print help, then exit + usage() + + elif switch == '-l': + logfp = open(val, "w") elif switch == '-m': threshold = int(val) elif switch == '-n': @@ -481,26 +1030,16 @@ if __name__ == '__main__': await = int(val[:-1]) * 360 else: await = int(val) + elif switch == '-r': + redo = True elif switch == '-t': + # replace title title = val + elif switch == '-S': + # add sub title + subtitle = val elif switch == '-T': terminal = val - elif switch == '-d': - dumpfile = val - elif switch == '-l': - logfp = open(val, "w") - elif switch == '-r': - redo = True - elif switch == '-D': - verbose = int(val) - elif switch == '-h': - sys.stderr.write( - "usage: gpsprof [-h] [-D debuglevel] [-m threshold] " - "[-n samplecount] [-d]\n" - + "\t[-f {" + "|".join([x.name for x in formatters]) - + "}] [-s speed] [-t title] [-T terminal] " - "[server[:port[:device]]]\n") - sys.exit(0) (host, port, device) = ("localhost", "2947", None) if len(arguments): @@ -536,9 +1075,16 @@ if __name__ == '__main__': # Ship the plot to standard output if not title: title = plot.whatami() + # escape " for gnuplot + title = title.replace('"', '\\"') + if subtitle: + title += '\\n' + subtitle if terminal: - sys.stdout.write("set terminal %s size 800,600\n" % terminal) - sys.stdout.write("set title \"%s\"\n" % title) + sys.stdout.write("set terminal %s truecolor enhanced size" + " 800,950\n" + % terminal) + # double quotes on title so \n is parsed by gnuplot + sys.stdout.write('set title noenhanced "%s"\n' % title) sys.stdout.write(plot.plot()) except KeyboardInterrupt: pass diff --git a/gpsprof.xml b/gpsprof.xml index d7151050..bdc22691 100644 --- a/gpsprof.xml +++ b/gpsprof.xml @@ -7,7 +7,7 @@ SPDX-License-Identifier: BSD-2-clause "-//OASIS//DTD DocBook XML V4.1.2//EN" "http://www.oasis-open.org/docbook/xml/4.1.2/docbookx.dtd"> -10 Feb 2005 +30 May 2018 gpsprof 1 @@ -22,16 +22,17 @@ SPDX-License-Identifier: BSD-2-clause gpsprof - -f plot_type - -m threshold - -n packetcount - -t title - -T terminal + -D debuglevel -d dumpfile + -f plot_type + -h -l logfile + -m threshold + -n samplecount -r - -D debuglevel - -h + -S subtitle + -T terminal + -t title [server[:port[:device]]] @@ -39,79 +40,94 @@ SPDX-License-Identifier: BSD-2-clause DESCRIPTION gpsprof performs accuracy, latency, -and time drift profiling on a GPS. It emits to standard output a -GNUPLOT program that draws one of several illustrative graphs. It can +skyview, and time drift profiling on a GPS. It emits to standard output +a GNUPLOT program that draws one of several illustrative graphs. It can also be told to emit the raw profile data. Information from the default spatial plot it provides can be -useful for establishing an upper bound on latency, and thus on -position accuracy of a GPS in motion. +useful for characterizing position accuracy of a GPS. gpsprof uses instrumentation built -into gpsd. - -To display the graph, use -gnuplot1. -Thus, for example, to display the default spatial scatter plot, do -this: - - -gpsprof | gnuplot -persist - - - -To generate an image file: - - -gpsprof -T png | gnuplot >image.png - +into gpsd. It can read data from a local +or remote running gpsd. Or it can read +data from a saved logfile. + +gpsprof is designed to be lightweight +and use minimal host resources. No graphics subsystem needs to be +installed on the host running gpsprof. Simply +copy the resultant plot file to another host to be rendered +with gnuplot. OPTIONS -The -f option sets the plot type. The X axis is samples (either -sentences with timestamps or PPS time drift messages). The Y axis is -normally latency in seconds, except for the spatial plot. Currently the -following plot types are defined: +The -f option sets the plot type. Currently the following plot +types are defined: space -Generate a scattergram of fixes and plot a probable-error -circle. This data is only meaningful if the GPS is held stationary -while gpsprof is running. -This is the default. +Generate a scatterplot of fixes and plot probable error circles. +This data is only meaningful if the GPS is held stationary while +gpsprof is running. Various statistics about +the fixes are listed at the bottom. This is the default plot type. -time +polar -Plot delta of system clock (NTP corrected time) against GPS time -as reported in PPS messages. +Generate a heat map of reported satellite Signal to Noise Ratio +(SNR) using polar coordinates. A colored dot is plotted for +each satellite seen by the GPS. The color of dot corresponds to the +SNR of the satellite. The dots are plotted by azimuth and +elevation. North, azimuth 0 degrees, is at the top of the plot. +Directly overhead, elevation of 90 degrees, is plotted at the center. +Useful for analyzing the quality of the skyview as seen by the GPS. + + -uninstrumented +polarunused -Plot total latency without instrumentation. Useful mainly as a -check that the instrumentation is not producing significant -distortion. It only plots times for reports that contain fixes; -staircase-like artifacts in the plot are created when elapsed time -from reports without fixes is lumped in. +Similar to the polar plot, but only unused satellites +are plotted. Useful for seeing which parts of the antenna skyview +are obstructed, degraded, below the GPS elevation mask, or otherwise +rejected. + -instrumented +polarused -Plot instrumented profile. -Plots various components of the total latency between the GPS's fix time -fix and when the client receives the fix. +Similar to the polar plot, but only satellites used to compute +fixs are plotted. Useful for seeing which parts of the antenna +skyview are being used in fixes. - + +time + +Plot delta of system clock (NTP corrected time) against GPS time +as reported in PPS messages. The X axis is sample time in seconds +from the start of the plot. The Y axis is the system clock delta from +GPS time. This plot only works if gpsd was +built with the timing (latency timing support) configure option +enabled. + + + + +instrumented + +Plot instrumented profile. Plots various components of the total +latency between the GPS's fix time and when the client receives the +fix. This plot only works if gpsd was built +with the timing (latency timing support) configuration option enabled. + For purposes of the description, below, start-of-reporting-cycle (SORC) is when a device's reporting cycle begins. This time is detected by watching to see when data availability follows a long @@ -161,34 +177,95 @@ beginning of a session. The -m option lets you set a latency threshold, in multiples of the cycle time, above which reports are discarded. -The -n option sets the number of packets to sample. The default -is 100. - -The -t option sets a text string to be included in the plot -title. + + + +uninstrumented + +Plot total latency without instrumentation. Useful mainly as a +check that the instrumentation is not producing significant distortion. +The X axis is sample time in seconds from the start of the plot. The Y +axs is latency in seconds.It only plots times for reports that contain +fixes; staircase-like artifacts in the plot are created when elapsed +time from reports without fixes is lumped in. + + + -The -T option generates a terminal type setting into the gnuplot code. -Typical usage is "-T png" telling gnuplot to write a PNG file. Without -this option gnuplot will call its X11 display code. The -d option dumps the plot data, without attached gnuplot code, to a specified file for post-analysis. +The -D sets debug level. + +The -h option makes gpsprof print +a usage message and exit. + The -l option dumps the raw JSON reports collected from the device to a specified file. +The -n option sets the number of packets to sample. The default +is 100. Most GPS are configured to emit one fix per second, so 100 +samples would then span 100 seconds. + The -r option replots from a JSON logfile (such as -l produces) on standard input. Both -n and -l options are ignored when this one is selected. -The -h option makes gpsprof print -a usage message and exit. +The -S option sets a text string to be included in the plot +as a subtitle. This will be below the title. -The -D sets debug level. +The -t option sets a text string to be the plot title. This +will replace the default title. +The -T option generates a terminal type setting into the gnuplot +code. Typical usage is "-T png", or "-T pngcairo" telling gnuplot to +write a PNG file. Without this option gnuplot will call its X11 display +code. + Different installations of gnuplot will +support different terminal types. Different terminal types may work better +for you than other ones. "-T png" will generate PNG images. Use "-T jpeg" +to generate JPEG images. "-T pngcairo" often works best, but is not +supported by some distributions. + + +SIGNALS Sending SIGUSR1 to a running instance causes it to write a completion message to standard error and resume processing. The first number in the startup message is the process ID to signal. + + +EXAMPLES +To display the graph, use +gnuplot1. +Thus, for example, to display the default spatial scatter plot, do +this: + + +gpsprof | gnuplot -persist + + + +To generate an image file: + + +gpsprof -T png | gnuplot > image.png + + + +To generate a polar plot, and save the GPS data for further plots: + +gpsprof -f polar -T jpeg -l polar.json | gnuplot > polar.png + +Then to make the matching polarused and polarunused plots and pngs from +the just saved the GPS data: + +gpsprof -f polarused -T jpeg -r < polar.json > polarused.plot +gnuplot < polarused.plot > polarused.png +gpsprof -f polarunused -T jpeg -r < polar.json > polarunused.plot +gnuplot < polarunused.plot > polarunused.png + + SEE ALSO -- cgit v1.2.1