gpsprof: Add statistics to scatterplot. Add satellite heat maps.

Plus the man page.

This work supported by the kind generatosity of Virgin Orbit.

1 files changed, 665 insertions, 119 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