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 ++++++++++++++++++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 665 insertions(+), 119 deletions(-) (limited to 'gpsprof') 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 -- cgit v1.2.1