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.
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(filename,"png");
Report( gb )<<updatewindow;
);
dt = Open( "H:\gpsLog3.txt",
columns(
Column( "Time", Numeric, "Continuous", Format( "yyyy-mm-ddThh:mm:ss", 19, 0 ), Input Format( "yyyy-mm-ddThh:mm:ss", 0 ) ),
Column( "lat1", Numeric, "Continuous", Format( "Best", 12 ) ),
Column( "lon1", Numeric, "Continuous", Format( "Best", 12 ) ),
Column( "lat2", Numeric, "Continuous", Format( "Best", 12 ) ),
Column( "lon2", Numeric, "Continuous", Format( "Best", 12 ) ),
Column( "lat3", Numeric, "Continuous", Format( "Best", 12 ) ),
Column( "lon3", Numeric, "Continuous", Format( "Best", 12 ) )
),
Import Settings(
End Of Line( CRLF, CR, LF ), End Of Field( Comma, CSV( 0 ) ), Strip Quotes( 1 ),
Use Apostrophe as Quotation Mark( 0 ), Use Regional Settings( 0 ), Scan Whole File( 1 ),
Treat empty columns as numeric( 0 ), CompressNumericColumns( 0 ),
CompressCharacterColumns( 0 ), CompressAllowListCheck( 0 ), Labels( 0 ),
Column Names Start( 1 ), Data Starts( 1 ), Lines To Read( "All" ), Year Rule( "20xx" )
)
);
dt<<selectwhere(lat1==0|lon1==0|lat2==0|lon2==0|lat3==0|lon3==0);
dt<<deleterows;
dt<<newcolumn("lat9",formula((lat1+lat2+lat3)/3));
dt<<newcolumn("lon9",formula((lon1+lon2+lon3)/3));
dt<<runformulas;
dt:lat9<<deleteformula();
dt:lon9<<deleteformula();
dt9=dt << Subset( All rows, columns( :Time, :lat9, :lon9 ) );
dt9:lat9<<setname("lat");
dt9:lon9<<setname("lon");
(dt9<<newcolumn("id"))<<setEachValue(9);
dt1=dt << Subset( All rows, columns( :Time, :lat1, :lon1 ) );
dt1:lat1<<setname("lat");
dt1:lon1<<setname("lon");
(dt1<<newcolumn("id"))<<setEachValue(1);
dt2=dt << Subset( All rows, columns( :Time, :lat2, :lon2 ) );
dt2:lat2<<setname("lat");
dt2:lon2<<setname("lon");
(dt2<<newcolumn("id"))<<setEachValue(2);
dt3=dt << Subset( All rows, columns( :Time, :lat3, :lon3 ) );
dt3:lat3<<setname("lat");
dt3:lon3<<setname("lon");
(dt3<<newcolumn("id"))<<setEachValue(3);
close(dt,"nosave");
dt4 = dt1 << Concatenate(
dt2,
dt3,
dt9
);
close(dt1,"nosave");
close(dt2,"nosave");
close(dt3,"nosave");
close(dt9,"nosave");
dt = dt4 << Sort( By( :Time ), Order( Ascending ) );
close(dt4,"nosave");
dtsummary = dt << Summary(
Group( :id ),
Mean( :lat ),
Std Dev( :lat ),
Min( :lat ),
Max( :lat ),
Mean( :lon ),
Std Dev( :lon ),
Min( :lon ),
Max( :lon ),
Freq( "None" ),
Weight( "None" )
);
lonTarget = dtsummary:name("Mean(lon)")[4]; // row 4 is id==9 for the average
latTarget = dtsummary:name("Mean(lat)")[4]; // row 4 is id==9 for the average
lonMin = min(dtsummary:name("Min(lon)")[1::3]); // rows 1 thru 3
lonMax = max(dtsummary:name("Max(lon)")[1::3]); // need the min/max
latMin = min(dtsummary:name("Min(lat)")[1::3]); // across sensors
latMax = max(dtsummary:name("Max(lat)")[1::3]);
//dummy=vlistbox(
gb = dt << Graph Builder(
Size( 1924, 1084 ), Show Control Panel( 0 ), Show Legend( 0 ), lock scales(1),
Variables( X( :lon ), Y( :lat ), Overlay( :id ) ), Elements( Points( X, Y, Legend( 4 ) ) ),
SendToReport( Dispatch( {}, "lon", ScaleBox,
{Scale( "Geodesic" ), Min(lonMin ), Max( lonMax), Inc( 0.0005 ), Minor Ticks( 0 ),
Add Ref Line( lonTarget, "Dotted", "Medium Dark Blue", "", 2 ), Label Row(
{Inside Ticks( 1 ), Show Major Labels( 0 ), Show Major Ticks( 0 ), Show Minor Ticks( 0 )}
)}
),
Dispatch( {}, "lat", ScaleBox,
{Scale( "Geodesic" ), Min( latMin ), Max( (latMax+latTarget)/2 /*adjust viewport-artistic license*/), Inc( 0.0005 ), Minor Ticks( 0 ),
Add Ref Line( latTarget, "Dotted", "Medium Dark Blue", "", 2 ), Label Row(
{Inside Ticks( 1 ), Show Major Labels( 0 ), Show Major Ticks( 0 ), Show Minor Ticks( 0 )}
)}
),
Dispatch( {}, "400", ScaleBox, {Legend Model( 4,
Properties( 0, {Line Color( colors[1] )} ),
Properties( 1, {Line Color( colors[2] )} ),
Properties( 2, {Line Color( colors[3] )} ),
Properties( 3, {Line Color( colors[4] )} )
)}
),
Dispatch( {}, "graph title", TextEditBox, {Set Text( "GPS tracker" )} ),
Dispatch( {}, "X title", TextEditBox, {Set Text( "" )} ),
Dispatch( {}, "Y title", TextEditBox, {Set Text( "" )} ),
Dispatch( {}, "Graph Builder", FrameBox,
{Background Map( Images( "Street Map Service" ) ), Transparency( 1 ),
Add Graphics Script( 8, Description( "Overlay" ),
Pen Color( white ); Pen Size( 4 ); // placeholder Circle( {lonTarget, latTarget}, .0001 );
),
Grid Line Order( 2 ), Reference Line Order( 7 )}
)
)
);
//);
takePicture();
// graph builder can add shaded ellipses without doing this...
For( i = 1, i <= 4, i += 1,
sensor = [1, 2, 3, 9][i];
biv = dt << Bivariate( Y( :lat ), X( :lon ), where( id == sensor ), Show Points( 0 ), Density Ellipse( 0.99, {Line Color( colors[i] )} ) );
Report( biv )[framebox( 1 )] << copyFrameContents;
Report( gb )[framebox( 1 )] << pasteFrameContents;
biv << closewindow;
takePicture();
);
Report( gb )[framebox( 1 )] << Marker Size( 0 );
For( trans = 1, trans >= .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;
);
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 <<runformula message makes sure they are evaluated, and then the formula is removed, because it won't make sense after the next step.
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.
The GPS Wander video.
You must be a registered user to add a comment. If you've already registered, sign in. Otherwise, register and sign in.