gpsd/clients/gpsprof.py.in

#!@PYSHEBANG@
#
# @GENERATED@
# This file is Copyright 2010 by the GPSD project
# SPDX-License-Identifier: BSD-2-clause
#
# 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"
#
# TODO: put date from data on plot, not time of replot.
# TODO: add lat/lon to polar plots
#
# This code runs compatibly under Python 2 and 3.x for x >= 2.
# Preserve this property!
# Codacy D203 and D211 conflict, I choose D203
# Codacy D212 and D213 conflict, I choose D212
"""Collect and plot latency-profiling data from a running gpsd.

Requires gnuplot, but gnuplot can be on another host.
"""

from __future__ import absolute_import, print_function, division

import argparse
import copy
import math
import os
import signal
import socket
import sys
import time

# pylint wants local modules last
try:
    import gps
    import gps.clienthelpers
except ImportError as e:
    sys.stderr.write(
        "gpsprof: can't load Python gps libraries -- check PYTHONPATH.\n")
    sys.stderr.write("%s\n" % e)
    sys.exit(1)

gps_version = '@VERSION@'
if gps.__version__ != gps_version:
    sys.stderr.write("gpsprof: ERROR: need gps module version %s, got %s\n" %
                     (gps_version, gps.__version__))
    sys.exit(1)


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."""

    def __init__(self, prompt, endmsg=None):
        """Init class baton."""
        self.stream = sys.stderr
        self.stream.write(prompt + "...")
        if os.isatty(self.stream.fileno()):
            self.stream.write(" \b")
        self.stream.flush()
        self.count = 0
        self.endmsg = endmsg
        self.time = time.time()

    def twirl(self, ch=None):
        """Twirl the baton."""
        if self.stream is None:
            return
        if ch:
            self.stream.write(ch)
        elif os.isatty(self.stream.fileno()):
            self.stream.write("-/|\\"[self.count % 4])
            self.stream.write("\b")
        self.count = self.count + 1
        self.stream.flush()
        return

    def end(self, msg=None):
        """Write the end message."""
        if msg is None:
            msg = self.endmsg
        if self.stream:
            self.stream.write("...(%2.2f sec) %s.\n"
                              % (time.time() - self.time, msg))


class stats(object):
    """Class for 1D stats: min, max, mean, sigma, skewness, kurtosis."""

    def __init__(self):
        """Init class stats."""
        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 __str__(self):
        """Return a nice string, for debug."""
        return ("min %f, max %f, mean %f, median %f, sigma %f, skewedness %f, "
                "kurtosis %f" %
                (self.min, self.max, self.mean, self.median,
                 self.sigma, self.skewness, self.kurtosis))

    def min_max_mean(self, fixes, index):
        """Find min, max, and mean of fixes[index]."""
        if not fixes:
            return

        # might be fast to go through list once?
        if isinstance(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(fixes) / len(fixes)
            self.min = min(fixes)
            self.max = max(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.
        self.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".
        self.kurtosis = float('nan')

        if not fixes:
            return

        m3 = 0.0
        m4 = 0.0
        if isinstance(fixes[0], tuple):
            sum_squares = [(x[index] - self.mean) ** 2 for x in fixes]
            sigma = math.sqrt(sum(sum_squares) / (len(fixes) - 1))
            for fix in fixes:
                m3 += pow(fix[index] - sigma, 3)
                m4 += pow(fix[index] - sigma, 4)
        else:
            # must be float
            sum_squares = [(x - self.mean) ** 2 for x in fixes]
            sigma = math.sqrt(sum(sum_squares) / (len(fixes) - 1))
            for fix in fixes:
                m3 += pow(fix - sigma, 3)
                m4 += pow(fix - sigma, 4)

        self.sigma = sigma
        if sigma > 0.0001:
            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."""
    requires_time = False  # Default

    def __init__(self):
        """Init class plotter."""
        self.device = []
        self.fixes = []
        self.in_replot = False
        self.session = None
        self.start_time = int(time.time())
        self.watch = set(['TPV'])

    def whatami(self):
        """How do we identify this plotting run?"""
        desc = "%s" % gps.misc.isotime(self.start_time)
        if 'driver' in self.device:
            desc += ", %s" % self.device['driver']
        if 'bps' in self.device:
            desc += ", %d %dN%d, cycle %.3gs" % \
                (self.device['bps'], 9 - self.device['stopbits'],
                 self.device['stopbits'], self.device['cycle'])
        if 'path' in self.device:
            desc += ", %s" % self.device['path']
        if 'subtype' in self.device:
            desc += "\\n%s" % self.device['subtype']

        return desc

    def collect(self, verb, log_fp=None):
        """Collect data from the GPS."""
        try:
            self.session = gps.gps(host=options.host, port=options.port,
                                   verbose=verb)
        except socket.error:
            sys.stderr.write("gpsprof: gpsd unreachable.\n")
            sys.exit(1)
        # Initialize
        self.session.read()
        if self.session.version is None:
            sys.stderr.write("gpsprof: requires gpsd to speak new protocol.\n")
            sys.exit(1)
        # Set parameters
        flags = gps.WATCH_ENABLE | gps.WATCH_JSON
        if self.requires_time:
            flags |= gps.WATCH_TIMING
        if options.device:
            flags |= gps.WATCH_DEVICE
        try:
            signal.signal(signal.SIGUSR1,
                          lambda empty, unused: sys.stderr.write(
                              "%d of %d (%d%%)..."
                              % (options.wait - countdown, options.wait,
                                 ((options.wait - countdown) * 100.0 /
                                  options.wait))))
            signal.siginterrupt(signal.SIGUSR1, False)
            self.session.stream(flags, options.device)
            baton = Baton("gpsprof: %d looking for fix" % os.getpid(), "done")
            countdown = options.wait
            basetime = time.time()
            while countdown > 0:
                if self.session.read() == -1:
                    sys.stderr.write("gpsprof: gpsd has vanished.\n")
                    sys.exit(1)
                baton.twirl()
                if self.session.data["class"] == "ERROR":
                    sys.stderr.write(" ERROR: %s.\n"
                                     % self.session.data["message"])
                    sys.exit(1)
                if self.session.data["class"] == "DEVICES":
                    if ((len(self.session.data["devices"]) != 1 and
                         not options.device)):
                        sys.stderr.write("ERROR: multiple devices connected, "
                                         "you must explicitly specify the "
                                         "device.\n")
                        sys.exit(1)
                    for i in range(len(self.session.data["devices"])):
                        self.device = copy.copy(
                            self.session.data["devices"][i])
                        if self.device['path'] == options.device:
                            break
                if self.session.data["class"] == "WATCH":
                    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.
                if log_fp:
                    log_fp.write(self.session.response)
                # Ignore everything but what we're told to
                if self.session.data["class"] not in self.watch:
                    continue
                # We can get some funky artifacts at start of self.session
                # apparently due to RS232 buffering effects. Ignore
                # them.
                if ((options.threshold and
                     (time.time() - basetime <
                      self.session.cycle * options.threshold))):
                    continue
                if self.session.fix.mode <= gps.MODE_NO_FIX:
                    continue
                if self.sample():
                    if countdown == options.wait:
                        sys.stderr.write("first fix in %.2fsec, gathering %d "
                                         "samples..."
                                         % (time.time() - basetime,
                                            options.wait))
                    countdown -= 1
            baton.end()
        finally:
            self.session.stream(gps.WATCH_DISABLE | gps.WATCH_TIMING)
            signal.signal(signal.SIGUSR1, signal.SIG_DFL)

    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:
            baton.twirl()
            self.session.unpack(line)
            if self.session.data["class"] == "DEVICES":
                self.device = copy.copy(self.session.data["devices"][0])
            elif self.session.data["class"] not in self.watch:
                continue
            self.sample()
        baton.end()

    def dump(self):
        """Dump the raw data for post-analysis."""
        return self.header() + self.data()


class spaceplot(plotter):
    """Spatial scattergram of fixes."""

    name = "space"
    requires_time = False

    def __init__(self):
        """Initialize class spaceplot."""
        plotter.__init__(self)
        self.centroid = None
        self.centroid_ecef = None
        self.recentered = []

    def sample(self):
        """Grab samples."""
        # Watch out for the NaN value from gps.py.
        if (((self.in_replot or self.session.valid) and
             self.session.data["class"] == "TPV" and
             # Skip TPV without an actual fix.
             self.session.data["mode"] >= gps.MODE_2D)):
            # get sat used count
            sats_used = 0
            for sat in self.session.satellites:
                if sat.used:
                    sats_used += 1

            if 'altHAE' not in self.session.data:
                self.session.data['altHAE'] = gps.NaN

            self.fixes.append((self.session.data['lat'],
                               self.session.data['lon'],
                               self.session.data['altHAE'], sats_used))
        return True

    def header(self):
        """Return header."""
        return "\n# Position uncertainty, %s\n" % self.whatami()

    def postprocess(self):
        """Postprocess the sample data."""
        return

    def data(self):
        """Format data for dump."""
        res = ""
        for i in range(len(self.recentered)):
            (lat, lon) = self.recentered[i][:2]
            (raw1, raw2, alt) = self.fixes[i]
            res += "%.9f\t%.9f\t%.9f\t%.9f\t%.9f\n" \
                % (lat, lon, raw1, raw2, alt)
        return res

    def plot(self):
        """Plot the data."""
        stat_lat = stats()
        stat_lon = stats()
        stat_alt = stats()
        stat_used = stats()
        # recentered stats
        stat_lat_r = stats()
        stat_lon_r = stats()
        stat_alt_r = stats()

        if not self.fixes:
            print("plot(): fixes array is empty\n")
            raise Exception('When plotting, there are no fixes.')

        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_used.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 positive 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]

        stat_lat_r.min_max_mean(self.recentered, 0)
        stat_lon_r.min_max_mean(self.recentered, 1)
        # compute sigma, skewness and kurtosis of lat/lon
        stat_lat_r.moments(self.recentered, 0)
        stat_lon_r.moments(self.recentered, 1)

        # CEP(50) calculated per RCC 261-00, Section 3.1.1
        # Note that this is not the distance such that 50% of the
        # points are within that distance; it is a distance calculated
        # from the standard deviation as specified by a standard, and
        # is likely a true CEP(50) if the distribution meets
        # expectations, such as normality.  However, it seems in
        # general to be close to the actual population CEP(50).
        calc_cep = 0.5887 * (stat_lat_r.sigma + stat_lon_r.sigma)

        # 2DRMS calculated per RCC 261-00, Section 3.1.4
        calc_2drms = 2 * math.sqrt(stat_lat_r.sigma ** 2 +
                                   stat_lon_r.sigma ** 2)

        # Note that the fencepost error situation in the following is
        # not entirely clear.  Issues include CEP(50) = x being
        # defined as 50% of the points being < x distance vs <= x
        # distance, and the interaction of truncation of len*0.5 vs
        # zero-based arrays.  (If you care about this because you are
        # concerned about your results, then you almost certainly are
        # not using enough samples.)

        # Compute actual CEP(50%) for our input set.  This makes no
        # assumptions about the distribution being normal, etc. -- it
        # is literally finding the distance such that 50% of our
        # samples are within that distance.  We simply find the
        # distance of the point whose index is half the number of
        # points, after having sorted the points by distance.
        cep_meters = gps.misc.EarthDistance(
            self.centroid[:2], self.fixes[int(len(self.fixes) * 0.50)][:2])

        # Compute actual CEP(95%) for our input set.  As above, but
        # index at 0.95 in the sorted points, so that 95% of points
        # are closer, 5% farther.
        cep95_meters = gps.misc.EarthDistance(
            self.centroid[:2], self.fixes[int(len(self.fixes) * 0.95)][:2])

        # Compute actual CEP(99%) for our input set.  As above, but
        # index at 0.99 in the sorted points, so that 99% of points
        # are closer, 1% farther.
        cep99_meters = gps.misc.EarthDistance(
            self.centroid[:2], self.fixes[int(len(self.fixes) * 0.99)][:2])

        # Compute actual CEP(100%) for our input set.  This is the
        # maximum distance from the centroid for any point in the
        # input set, and hence the distance of last point in sorted order.
        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
        dist_3d_max = 0.0
        alt_fixes = []
        alt_fixes_r = []
        latlon_data = ""
        alt_data = ""
        # init calculated hep, sep and mrse
        calc_hep = gps.NaN
        calc_sep = gps.NaN
        calc_mrse = gps.NaN

        # grab and format the fixes as gnuplot will use them
        for i in range(len(self.recentered)):
            # 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)
                # micro meters should be good enough
                alt_data += "%.6f\n" % (alt)

        if alt_fixes:
            # got altitude data

            # Convert fixes to offsets from avg in meters
            alt_data_centered = ""

            # find min, max and mean of altitude
            stat_alt.min_max_mean(alt_fixes, 0)
            for alt in alt_fixes:
                alt_fixes_r.append(alt - stat_alt.mean)
                alt_data_centered += "%.6f\n" % (alt - stat_alt.mean)

            stat_alt_r.min_max_mean(alt_fixes_r, 0)
            stat_alt_r.moments(alt_fixes_r, 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 fix_lla in self.fixes:
                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

            # Sort fixes by distance from average altitude
            alt_fixes_r.sort(key=abs)
            # so we can rank fixes for EPs
            alt_ep = abs(alt_fixes_r[int(len(alt_fixes_r) * 0.50)])
            alt_ep95 = abs(alt_fixes_r[int(len(alt_fixes_r) * 0.95)])
            alt_ep99 = abs(alt_fixes_r[int(len(alt_fixes_r) * 0.99)])

            # HEP(50) calculated per RCC 261-00, Section 3.1.2
            calc_hep = 0.6745 * stat_alt_r.sigma

            # SEP(50) calculated per RCC 261-00, Section 3.1.3 (3)
            calc_sep = 0.51 * (stat_lat_r.sigma +
                               stat_lon_r.sigma +
                               stat_alt_r.sigma)

            # MRSE calculated per RCC 261-00, Section 3.1.5
            calc_mrse = math.sqrt(stat_lat_r.sigma ** 2 +
                                  stat_lon_r.sigma ** 2 +
                                  stat_alt_r.sigma ** 2)

        fmt_lab11a = ('hep = %.3f meters\\n'
                      'sep = %.3f meters\\n'
                      'mrse = %.3f meters\\n'
                      ) % (calc_hep, calc_sep, calc_mrse)

        if self.centroid[0] < 0.0:
            latstring = "%.9fS" % -self.centroid[0]
        elif stat_lat.mean > 0.0:
            latstring = "%.9fN" % self.centroid[0]
        else:
            latstring = "0.0"

        if self.centroid[1] < 0.0:
            lonstring = "%.9fW" % -self.centroid[1]
        elif stat_lon.mean > 0.0:
            lonstring = "%.9fE" % self.centroid[1]
        else:
            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'
                'AltHAE\\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'
                ) % ('{0:>15}'.format(latstring),
                     '{0:>15}'.format(lonstring),
                     '{0:>15.3f}'.format(stat_alt.mean),
                     '{0:>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'
                ) % ('{0:>10.3f}'.format(lat_min_o),
                     '{0:>10.3f}'.format(lat_max_o),
                     '{0:>10.3f}'.format(stat_lat_r.sigma),
                     '{0:>10.1f}'.format(stat_lat_r.skewness),
                     '{0:>10.1f}'.format(stat_lat_r.kurtosis),
                     '{0:>10.3f}'.format(lon_min_o),
                     '{0:>10.3f}'.format(lon_max_o),
                     '{0:>10.3f}'.format(stat_lon_r.sigma),
                     '{0:>10.1f}'.format(stat_lon_r.skewness),
                     '{0:>10.1f}'.format(stat_lon_r.kurtosis),
                     '{0:>10.3f}'.format(stat_alt_r.min),
                     '{0:>10.3f}'.format(stat_alt_r.max),
                     '{0:>10.3f}'.format(stat_alt_r.sigma),
                     '{0:>10.1f}'.format(stat_alt_r.skewness),
                     '{0:>10.1f}'.format(stat_alt_r.kurtosis),
                     stat_used.min,
                     stat_used.max,
                     '{0:>10.1f}'.format(stat_used.sigma))

        if 1 < options.debug:
            fmt += ('set label 15 at screen 0.6, screen 0.12 '
                    '"\\n'
                    '      min\\n'
                    '%s\\n'
                    '%s\\n'
                    '%s"\n'
                    ) % ('{0:>15.9f}'.format(stat_lat_r.min),
                         '{0:>15.9f}'.format(stat_lon_r.min),
                         '{0:>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'
                    ) % ('{0:>15.9f}'.format(stat_lat_r.max),
                         '{0:>15.9f}'.format(stat_lon_r.max),
                         '{0:>15.3f}'.format(stat_alt.max))

        if len(self.fixes) > 1000:
            plot_style = 'dots'
        else:
            plot_style = 'points'

        # 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 "AltHAE 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 AltHAE 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
    seen_used = []          # count of seen and used in each SKY

    def __init__(self):
        """Init class polarplot."""
        plotter.__init__(self)
        self.watch = set(['SKY'])

    def sample(self):
        """Grab samples."""
        if self.session.data["class"] == "SKY":
            sats = self.session.data['satellites']
            seen = 0
            used = 0
            for sat in sats:
                seen += 1
                # u'ss': 42, u'el': 15, u'PRN': 18, u'az': 80, u'used': True
                if ((('az' not in sat) or
                     ('el' not in sat))):
                    continue
                if sat['used'] is True:
                    used += 1
                    if 'polarunused' == self.name:
                        continue
                if (('polarused' == self.name) and (sat['used'] is False)):
                    continue
                self.fixes.append((sat['PRN'], sat['ss'], sat['az'],
                                  sat['el'], sat['used']))
            self.seen_used.append((seen, used))

        return True

    def header(self):
        """Return header."""
        return "# Polar plot of signal strengths, %s\n" % self.whatami()

    def postprocess(self):
        """Postprocess the sample data."""
        return

    def data(self):
        """Format data for dump."""
        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):
        """Format data for dump."""
        # calc SNR: mean, min, max, sigma
        stat_ss = stats()
        stat_ss.min_max_mean(self.fixes, 1)
        stat_ss.moments(self.fixes, 1)

        # calc sats seen data: mean, min, max, sigma
        stat_seen = stats()
        stat_seen.min_max_mean(self.seen_used, 0)
        stat_seen.moments(self.seen_used, 0)

        # calc sats used data: mean, min, max, sigma
        stat_used = stats()
        stat_used.min_max_mean(self.seen_used, 1)
        stat_used.moments(self.seen_used, 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 increment 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 rmargin then put colorbox in the margin.
set lmargin at screen .08
set rmargin at screen .85
set cbrange [10:60]
set palette defined (100 "blue", 200 "green", 300 "red")
set colorbox user origin .92, .15 size .03, .6
set label 10 at screen 0.89, screen 0.13 "SNR dBHz"
'''

        count = len(self.fixes)
        fmt += '''\
set label 11 at screen 0.01, screen 0.15 "%s plot, samples %d"
set label 12 at screen 0.01, screen 0.10 "\\nSS\\nSeen\\nUsed"
''' % (self.name, count)

        fmt += '''\
set label 13 at screen 0.11, screen 0.10 "min\\n%d\\n%d\\n%d" right
''' % (stat_ss.min, stat_seen.min, stat_used.min)

        fmt += '''\
set label 14 at screen 0.21, screen 0.10 "max\\n%d\\n%d\\n%d" right
''' % (stat_ss.max, stat_seen.max, stat_used.max)

        fmt += '''\
set label 15 at screen 0.31, screen 0.10 "mean\\n%.1f\\n%.1f\\n%.1f" right
''' % (stat_ss.mean, stat_seen.mean, stat_used.mean)

        fmt += '''\
set label 16 at screen 0.41, screen 0.10 "sigma\\n%.1f\\n%.1f\\n%.1f" right
''' % (stat_ss.sigma, stat_seen.sigma, stat_used.sigma)

        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"
    requires_time = True

    def __init__(self):
        """Init class timeplot."""
        plotter.__init__(self)
        self.watch = set(['PPS'])

    def sample(self):
        """Grab samples."""
        if self.session.data["class"] == "PPS":
            self.fixes.append((self.session.data['real_sec'],
                               self.session.data['real_nsec'],
                               self.session.data['clock_sec'],
                               self.session.data['clock_nsec']))
        return True

    def header(self):
        """Return header."""
        return "# Time drift against PPS, %s\n" % self.whatami()

    def postprocess(self):
        """Postprocess the sample data."""
        return

    def data(self):
        """Format data for dump."""
        res = ""
        for (real_sec, real_nsec, clock_sec, clock_nsec) in self.fixes:
            res += "%d\t%d\t%d\t%d\n" % (real_sec, real_nsec, clock_sec,
                                         clock_nsec)
        return res

    def plot(self):
        """Format data for dump."""
        fmt = '''\
set autoscale
set key below
set ylabel "System clock delta from GPS time (nsec)"
plot "-" using 0:((column(1)-column(3))*1e9 + (column(2)-column(4))) \
 title "Delta" with impulses
'''
        return fmt + self.header() + self.data()


class uninstrumented(plotter):
    """Total times without instrumentation."""
    name = "uninstrumented"
    requires_time = False

    def __init__(self):
        """Init class uninstrumentd."""
        plotter.__init__(self)

    def sample(self):
        """Grab samples."""
        if self.session.fix.time:
            seconds = time.time() - gps.misc.isotime(self.session.data.time)
            self.fixes.append(seconds)
            return True

        return False

    def header(self):
        """Return header."""
        return "# Uninstrumented total latency, " + self.whatami() + "\n"

    def postprocess(self):
        """Postprocess the sample data."""
        return

    def data(self):
        """Format data for dump."""
        res = ""
        for seconds in self.fixes:
            res += "%2.6lf\n" % seconds
        return res

    def plot(self):
        """Plot the data."""
        fmt = '''\
set autoscale
set key below
set key title "Uninstrumented total latency"
plot "-" using 0:1 title "Total time" with impulses
'''
        return fmt + self.header() + self.data()


class instrumented(plotter):
    """Latency as analyzed by instrumentation."""

    name = "instrumented"
    requires_time = True

    def __init__(self):
        """Initialize class instrumented()."""
        plotter.__init__(self)

    def sample(self):
        """Grab the samples."""
        if 'rtime' in self.session.data:
            self.fixes.append((gps.misc.isotime(self.session.data['time']),
                               self.session.data["chars"],
                               self.session.data['sats'],
                               self.session.data['sor'],
                               self.session.data['rtime'],
                               time.time()))
            return True

        return False

    def header(self):
        """Return the header."""
        res = "# Analyzed latency, " + self.whatami() + "\n"
        res += "#-- Fix time --  - Chars -  --   Latency  - RS232-  " \
               "Analysis  - Recv -\n"
        return res

    def postprocess(self):
        """Postprocess the sample data."""
        return

    def data(self):
        """Format data for dump."""
        res = ""
        for (fix_time, chars, sats, start, xmit, recv) in self.fixes:
            rs232_time = (chars * 10.0) / self.device['bps']
            res += "%.3f  %9u  %2u  %.6f  %.6f  %.6f  %.6f\n" \
                % (fix_time, chars, sats, start - fix_time,
                   (start - fix_time) + rs232_time, xmit - fix_time,
                   recv - fix_time)
        return res

    def plot(self):
        """Do the plot."""
        legends = (
            "Reception delta",
            "Analysis time",
            "RS232 time",
            "Fix latency",
        )
        fmt = '''\
set autoscale
set key title "Analyzed latency"
set key below
plot \\\n'''
        for (i, legend) in enumerate(legends):
            j = len(legends) - i + 4
            fmt += '    "-" using 0:%d title "%s" with impulses, \\\n' \
                % (j, legend)
        fmt = fmt[:-4] + "\n"
        return fmt + self.header() + (self.data() + "e\n") * len(legends)


if __name__ == '__main__':
    usage = '%(prog)s [OPTIONS] [host[:port[:device]]]'
    epilog = ('BSD terms apply: see the file COPYING in the distribution '
              'root for details.')

    # get default units from the environment
    # GPSD_UNITS, LC_MEASUREMENT and LANG
    # FIXME: use default_units
    default_units = gps.clienthelpers.unit_adjustments()

    parser = argparse.ArgumentParser(usage=usage, epilog=epilog)
    parser.add_argument(
        '-?',
        action="help",
        help='show this help message and exit'
    )
    parser.add_argument(
        '-D',
        '--debug',
        dest='debug',
        default=0,
        type=int,
        help='Set level of debug. Must be integer. [Default %(default)s]',
    )
    parser.add_argument(
        '-d', '--dumpfile',
        dest='dumpfile',
        default=None,
        help='Set dumpfile. [Default %(default)s]',
    )
    parser.add_argument(
        '--device',
        dest='device',
        default='',
        help='The device to connect. [Default %(default)s]',
    )
    parser.add_argument(
        '-f', '--formatter',
        choices=('instrumented',
                 'polar',
                 'polarunused',
                 'polarused',
                 'space',
                 'time',
                 'uninstrumented',
                 ),
        dest='plotmode',
        default='space',
        help='Formatter to use. [Default %(default)s]',
    )
    parser.add_argument(
        '--host',
        dest='host',
        default='localhost',
        help='The host to connect. [Default %(default)s]',
    )
    parser.add_argument(
        '-l',
        '--logfile',
        dest='logfile',
        default=None,
        help='Select logfile for JSON. [Default %(default)s]',
    )
    parser.add_argument(
        '-m', '--threshold',
        dest='threshold',
        default=None,
        help='Select logfile. [Default %(default)s]',
        type=int,
    )
    parser.add_argument(
        '-n', '--wait',
        dest='wait',
        default='100',
        help=('Select wait time in seconds before plotting.\n'
              '[Default %(default)s]'),
        type=str,
    )
    parser.add_argument(
        '--port',
        dest='port',
        default=gps.GPSD_PORT,
        help='The port to connect. [Default %(default)s]',
    )
    parser.add_argument(
        '-r', '--redo',
        action='store_true',
        dest='redo',
        default=False,
        help=('Redo from a previous JSON logfile on stdin.\n'
              '[Default %(default)s]'),
    )
    parser.add_argument(
        '-S', '--subtitle',
        dest='subtitle',
        default=None,
        help='Plot subtitle. [Default %(default)s]',
    )
    parser.add_argument(
        '-t', '--title',
        dest='title',
        default=None,
        help='Plot title. [Default %(default)s]',
    )
    parser.add_argument(
        '-T', '--terminal',
        dest='terminal',
        default='x11',
        help='gnuplot terminal type. [Default %(default)s]',
    )
    parser.add_argument(
        '-V', '--version',
        action='version',
        version="%(prog)s: Version " + gps_version + "\n",
        help='Output version to stderr, then exit',
    )
    parser.add_argument(
        'target',
        nargs='?',
        help='[host[:port[:device]]]',
    )

    options = parser.parse_args()

    if options.logfile:
        options.logfp = open(options.logfile, "w")
    else:
        options.logfp = None

    # the options host, port, device are set by the defaults
    if options.target:
        # override with target
        arg = options.target.split(':')
        len_arg = len(arg)
        if len_arg == 1:
            (options.host,) = arg
        elif len_arg == 2:
            (options.host, options.port) = arg
        elif len_arg == 3:
            (options.host, options.port, options.device) = arg
        else:
            parser.print_help()
            sys.exit(0)

    if not options.port:
        options.port = gps.GPSD_PORT

    formatters = (polarplot, polarplotunused, polarplotused,
                  spaceplot, timeplot, uninstrumented, instrumented)
    if options.plotmode:
        for formatter in formatters:
            if formatter.name == options.plotmode:
                plot = formatter()
                break
        else:
            sys.stderr.write("gpsprof: no such formatter.\n")
            sys.exit(1)

    if options.wait[-1] == 'h':
        options.wait = int(options.wait[:-1]) * 360
    else:
        options.wait = int(options.wait)

    try:
        # Get fix data from the GPS
        if options.redo:
            plot.replot(sys.stdin)
        else:
            plot.collect(options.debug, options.logfp)
        plot.postprocess()
        # Save the timing data (only) for post-analysis if required.
        if options.dumpfile:
            with open(options.dumpfile, "w") as fp:
                fp.write(plot.dump())
        if options.logfp:
            options.logfp.close()
        # Ship the plot to standard output
        if not options.title:
            options.title = plot.whatami()
            # escape " for gnuplot
            options.title = options.title.replace('"', '\\"')
        if options.subtitle:
            options.title += '\\n' + options.subtitle
        term_opts = ""
        truecolor_terms = ['png', 'sixelgd', 'wxt']
        if options.terminal in truecolor_terms:
            term_opts = 'truecolor'
        sys.stdout.write("set terminal %s size 800,950 %s\n"
                         "set termoption enhanced\n"
                         % (options.terminal, term_opts))
        # double quotes on title so \n is parsed by gnuplot
        sys.stdout.write('set title noenhanced "%s\\n\\n"\n' % options.title)
        sys.stdout.write(plot.plot())
    except KeyboardInterrupt:
        pass

# The following sets edit modes for GNU EMACS
# Local Variables:
# mode:python
# End: