Surveyors use differential GPS systems to get accurate locations. A stationary GPS and a rover work together to cancel out some errors. I thought I could build a poor-man's differential GPS with two inexpensive GPSs by collecting the data from both and subtracting the motion vectors of the stationary device from the rover. But this does not work with consumer grade GPSs. If I'd read up first, I might have gone directly to...
Averaging several GPS outputs from adjacent devices does produce a better answer than a single device. I wound up with three on a window sill, feeding a Linux machine via USB ports.
USB GPSs in a windowUSB GPSs in a windowUSB GPSs in a windowUSB GPSs in a windowUSB GPSs in a window
I made another mistake when I wrote the Python code that reads the devices, near the end:
### edited for presentation, may not be functional without some TLC
# dmesg after plug in gives usb port info
import serial
import math
import re
import datetime
import sys
# https://github.com/Knio/pynmea2/blob/master/pynmea2/nmea_utils.py
def timestamp(s):
'''
Converts a timestamp given in "hhmmss[.ss]" ASCII format to a
datetime.time object
'''
ms_s = s[6:]
ms = ms_s and int(float(ms_s) * 1000000) or 0
t = datetime.time(
hour=int(s[0:2]),
minute=int(s[2:4]),
second=int(s[4:6]),
microsecond=ms)
return t
def datestamp(s):
'''
Converts a datestamp given in "DDMMYY" ASCII format to a
datetime.datetime object
'''
return datetime.datetime.strptime(s, '%d%m%y').date()
def dm_to_sd(dm):
'''
Converts a geographic coordiante given in "degres/minutes" dddmm.mmmm
format (ie, "12319.943281" = 123 degrees, 19.953281 minutes) to a signed
decimal (python float) format
'''
# '12319.943281'
if not dm or dm == '0':
return 0.
d, m = re.match(r'^(\d+)(\d\d\.\d+)$', dm).groups()
return float(d) + float(m) / 60
ser0 = serial.Serial('/dev/ttyACM0', 2*9600)
ser1 = serial.Serial('/dev/ttyACM1', 2*9600)
ser2 = serial.Serial('/dev/ttyACM2', 2*9600)
lat0=0
lon0=0
lat1=0
lon1=0
lat2=0
lon2=0
time = ""
while True :
state=ser0.readline().split(",")
if state[0]=="$GPRMC":
lat0 = dm_to_sd(state[3])
lon0 = -dm_to_sd(state[5])
time = datetime.datetime.combine(datestamp(state[9]), timestamp(state[1]))
state=ser1.readline().split(",")
if state[0]=="$GPRMC":
lat1 = dm_to_sd(state[3])
lon1 = -dm_to_sd(state[5])
time = datetime.datetime.combine(datestamp(state[9]), timestamp(state[1]))
state=ser2.readline().split(",")
if state[0]=="$GPRMC":
lat2 = dm_to_sd(state[3])
lon2 = -dm_to_sd(state[5])
time = datetime.datetime.combine(datestamp(state[9]), timestamp(state[1]))
print time.strftime("%Y-%m-%dT%H:%M:%S")+","+str(lat0)+","+str(lon0)+","+str(lat1)+","+str(lon1)+","+str(lat2)+","+str(lon2)
sys.stdout.flush()
The code loops forever, reading each of three devices, looking for the GPRMC message that has the interesting data, and printing a formatted message to the output file. In retrospect that simple print statement should have formatted the data to be the shape I eventually wanted. What I have is a file like
2017-06-29T13:59:59, 35.111,-78.222,35.112,-78.223,35.113,-78.224
where 35 and -78 are latitude and longitude, approximately. But what I'm going to want is quite a bit rearranged:
2017-06-29T13:59:59, 35.111,-78.222,1
2017-06-29T13:59:59, 35.112,-78.223,2
2017-06-29T13:59:59, 35.113,-78.224,3
where the trailing 1,2,3 identifies the sensor. Since I might import the data again, I wrote this JSL to do the conversion. (Actually most of the code is to produce the images, notes at the end.)
colors={{255,0,0},{0,255,0},{0,0,255},{0,0,0}};
ipicture=1e7; // sequence number with leading digit
path="H:\GPSWander\images\";
takePicture = function({},
// {filename=path||"gps"||char(ipicture++)||".png"},
// Report( gb )[framebox( 1 )]<<SAVEPICTURE>= .003, trans = trans * .80,
Report( gb )[framebox( 1 )] << Transparency( trans );
takePicture();
);
Report( gb )[framebox( 1 )] << Transparency( 1.0 );
Report( gb )[framebox( 1 )] << Marker Size( 2 );
dt << selectrows( 5 :: N Rows( dt ) );
dt << hide( 1 );
dt << clearselect;
step = N Rows( dt ) / 10000; // frames in video
For( i = 1, i < N Rows( dt ), i += step,
dt << selectrows( i :: i + step + step + step + step ); // trail
dt << hide( 0 );
takePicture();
dt << hide( 1 );
dt << clearselect;
);
</SAVEPICTURE>
Look for the open function for "gpsLog3.txt". It has the script I copied from a data table the import preview helped me with. The next few lines delete some rows with invalid data. There were only 6 bad ones; I'm not sure why the GPS on channel 2 reported 0,0.
The next few lines create lat9,lon9 as the average of lat1,2,3 and lon1,2,3. A <
The next step makes subsets with just the columns to be concatenated, renames them, adds the id variable. Then the tables are concatenated, unneeded tables closed.
The summary table that gets generated next is used to set reflines and axis limits in the GraphBuilder that follows. Then there are some loops to drive an animation sequence in graph builder (and the odd function at the beginning to take pictures.)
- I borrowed some github code to handle some of the NMEA GPS issues, thanks Knio!
- If you try to do this, you'll probably need "sudo adduser $USER dialout" too, look it up.
Explaining the density ellipse in the videoExplaining the density ellipse in the videoExplaining the density ellipse in the videoExplaining the density ellipse in the videoExplaining the density ellipse in the video
The GPS Wander video.
