GPS accuracy

This is a plot of the GPS signal from my office window over the course of a few minutes.

'Stationary' GPS

‘Stationary’ GPS

I’m seriously not impressed – the units are meters – typically this is from 3 or 4 satellites.  I’ll try a run in the garden to see what the results are like in a more open space.  The GPS output also includes a measure of predicted accuracy of each sample – I need to work out a way get get excel to mark a centre dot, with a surrounding error circle.

The GPS receiver is by Crius running a ublox chip and Sarantel geohelix antenna I got from ebay a few years ago when Sarantel went bust.

The graph is based upon the change in distance compared to the initial GPS input. Here’s the code – one bug fixed to get longitude and latitude the right way round!

from __future__ import division
import gps
import os
import math

# Listen on port 2947 (gpsd) of localhost
session = gps.gps("localhost", "2947")
session.stream(gps.WATCH_ENABLE | gps.WATCH_NEWSTYLE)

num_sats = 0
latitude = 0.0
longitude = 0.0
time = ""
epx = 0.0
epy = 0.0
epv = 0.0
ept = 0.0
eps = 0.0
climb = 0.0
altitude = 0.0
speed = 0.0
direction = 0.0

lat = 0.0
lon = 0.0
new_lat = False
new_lon = False
base_lon = 0.0
base_lat = 0.0

dx = 0.0
dy = 0.0

R = 6371000 # radius of the earth in meters

fp_name = "gpstrack.csv"
header = "time, latitude, longitude, satellites, climb, altitude, speed, direction, dx, dy"

os.system("clear")

print header

with open(fp_name, "wb") as fp:
    fp.write(header + "\n")

    #---------------------------------------------------------------------------------
    # With a based level longitude and latitude in degrees, we can be the current X and Y coordinates
    # relative to the takeoff position thus:
    # psi = latitude => p below
    # lambda = longitude => l below
    # Using equirectangular approximation: 
    #
    # x = (l2 - l1) * cos ((p1 + p2) / 2)
    # y = (p2 - p1)
    # d = R * (x*x + y*y) ^ 0.5
    # 
    # More at http://www.movable-type.co.uk/scripts/latlong.html
    #---------------------------------------------------------------------------------

    while True:
        try:
            report = session.next()
#            print report
#            os.system("clear")
            if report['class'] == 'TPV':
                if hasattr(report, 'time'):  # Time 
                    time = report.time

                if hasattr(report, 'ept'):   # Estimated timestamp error - seconds 
                    ept = report.ept

                if hasattr(report, 'lon'):   # Longitude in degrees
                    longitude = report.lon
                    new_lon = True

                if hasattr(report, 'epx'):   # Estimated longitude error - meters 
                    epx = report.epx

                if hasattr(report, 'lat'):   # Latitude in degrees
                    latitude = report.lat
                    new_lat = True

                if hasattr(report, 'epy'):   # Estimated latitude error - meters 
                    epy = report.epy

                if hasattr(report, 'alt'):   # Altitude - meters
                    altitude = report.alt
                if hasattr(report, 'epv'):   # Estimated altitude error - meters
                    epv = report.epv

                if hasattr(report, 'track'): # Direction - degrees from true north
                    direction = report.track     
                if hasattr(report, 'epd'):   # Estimated direction error - degrees 
                    epd = report.epd

                if hasattr(report, 'climb'): # Climb velocity - meters per second
                    climb = report.climb
                if hasattr(report, 'epc'):   # Estimated climb error - meters per seconds
                    epc = report.epc

                if hasattr(report, 'speed'): # Speed over ground - meters per second
                    speed = report.speed     
                if hasattr(report, 'eps'):   # Estimated speed error - meters per second
                    eps = report.eps


            if report['class'] == 'SKY':
                if hasattr(report, 'satellites'):
                    num_sats = 0
                    for satellite in report.satellites:
                        if hasattr(satellite, 'used') and satellite.used:
                            num_sats += 1

            #-----------------------------------------------------------------------------
            # Calculate the X,Y coordinates in meters
            #-----------------------------------------------------------------------------
            if new_lon and new_lat:

                new_lon = False
                new_lat = False
 
                lat = latitude * math.pi / 180
                lon = longitude * math.pi / 180

                if base_lat == 0.0 and base_lon == 0.0:
                    base_lat = lat
                    base_lon = lon
                    continue

                dx = (lon - base_lon) * math.cos((lat + base_lat) / 2) * R
                dy = (lat - base_lat) * R
     
            else:
                continue

            
            output = "%s, %f, %f, %d, %f, %f, %f, %f, %f, %f" % (time, 
                                                                 latitude,
                                                                 longitude,
                                                                 num_sats,
                                                                 climb, 
                                                                 altitude, 
                                                                 speed,
                                                                 direction,
                                                                 dx,
                                                                 dy)




            print output
            fp.write(output + "\n")
        except KeyError:
            pass
        except KeyboardInterrupt:
            quit()
        except StopIteration:
            session = None
            print "GPSD has terminated"

 

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.