# GPS accuracy

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

‘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")

with open(fp_name, "wb") as fp:

#---------------------------------------------------------------------------------
# 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"
```

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