#!/usr/bin/perl -w
#errorplot.pl: processes NMEA data to produce GNUPlot-ready error plot data
#Mon Feb 27 08:50:13 BRT 2006 - Durval Menezes.
#$Id: errorplot.pl,v 1.4 2006/02/28 00:18:35 durval Exp durval $

use strict;	#don't leave home without it

my ($lat, $lng, $hdop, @lat, @lng, @hdop, $hdop_class);
my ($lat_sum, $lng_sum, $nsamples, $lat_avg, $lng_avg);
my ($lat_delta_m, $lng_delta_m, $ix, $set);

my $DEG2M = 105520.917000265180;	#degrees/meters, obtained empirically

MAIN:
{

#we are interested in GGA records, example:
#$GPGGA,112931.000,2255.2405,S,04303.0489,W,1,06,3.7,232.4,M,-5.7,M,,0000*42
#$1     $2         $3       $4 $5        $6 $7 $8 $9
#$1: record type ('$GPGGA')
#$2: current time (GMT+0)
#$3: latitude (DDMM.MMMM)
#$4: north/south indicator ('N' or 'S')
#$5: longitude (DDDMM.MMMM)
#$6: east/west indicator ('E' or 'W')
#$7: fix flag ('0': invalid, '1': valid SPS, '2': valid DGPS, '3': valid PPS)
#$8: number of satellites being used
#$9: HDOP (Horizontal Dilution Of Precision)
#The other fields do not interest us
#source: http://www.commlinx.com.au/NMEA_sentences.htm

#run thru the file, converting/loading lat/lng/hdop into arrays and
#calculating averages for lat/lng
while(<>)
{
	chomp;
	#only process valid records
	if (/^(\$GPGGA),([0-9.]+),([0-9.]+),([SN]),([0-9.]+),([EW]),([0123]),(\d\d),([0-9.]+),[-0-9.]+,M,[-0-9.]+,M,,0000\*[0-9A-F]{2}\r{0,1}$/ && $7 != 0)
	{
		#extract lat/lng/hdop, converting the lat/lng to decimal degrees
		$lat = substr($3,0,2)+substr($3,2)/60.; $lat *= -1 if $4 eq 'S';
		$lng = substr($5,0,3)+substr($5,3)/60.; $lng *= -1 if $6 eq 'W';
		$hdop = $9;

		#push them at the end of the respective arrays
		push @lat, $lat;
		push @lng, $lng;
		push @hdop, $hdop;

		#accumulate sums for average calculations
		$lat_sum += $lat;
		$lng_sum += $lng;
		++$nsamples;
	}
}

#calculate averages
$lat_avg = $lat_sum / $nsamples;
$lng_avg = $lng_sum / $nsamples;

#now generate output
#print header
print <<EOM;
#GPS error plot data for gnuplot
#Fields are:
#  lat_delta  lng_delta  tot_delta  lat  lng  hdop  HDOP_class
#  lat_delta, lng_deltas: delta from average lat/lng in meters
#	obtained by multiplying delta degrees by $DEG2M);
#  tot_delta: total (straight line) delta, calculated by classic
#	pythagoras: sqrt(lat_delta^2 + lng_delta^2);
#  lat,lng,hdop: original GPS-generated latitude, longitude and HDOP
#  HDOP_class: A if HDOP <=1.0, B <=1.0, C <= 2.0, D <=4.0, E >4.0;
#deltas are relative to from average lat/lng of the whole dataset:
#  lat_avg : $lat_avg
#  lng_avg : $lng_avg
#  nsamples: $nsamples
#now let's plot
set title 'GPS error plot'
set xlabel 'Longitude (East/West) error (m)'
set ylabel 'Latitude (North/South) error (m)'
set size square
#plot first the worst precision sets (more dots) so they won't cover
#the best precision ones (fewer dots)
plot '-' using 2:1 title 'HDOP >4.0' with points 1,	\\
     '-' using 2:1 title 'HDOP <=4.0' with points 15, \\
	 '-' using 2:1 title 'HDOP <=2.0' with points 9,	\\
	 '-' using 2:1 title 'HDOP <=1.5' with points 2,	\\
	 '-' using 2:1 title 'HDOP <=1.0' with points 3;
EOM
#print records
for ($set=4; $set >= 0; --$set)
{
	for ($ix=0; $ix <$nsamples; ++$ix)
	{
		if 	($hdop[$ix] <= 1.0) {
			$hdop_class = 'A';
		}
		elsif	($hdop[$ix] <= 1.5) {
			$hdop_class = 'B';
		}
		elsif	($hdop[$ix] <= 2.0) {
			$hdop_class = 'C'
		}
		elsif	($hdop[$ix] <= 4.0) {
			$hdop_class = 'D'
		}
		else { 	#>4.0
			$hdop_class = 'E';
		}

		if ((ord($hdop_class) - ord('A')) == $set)
		{
			$lat_delta_m = ($lat[$ix] - $lat_avg) * $DEG2M;
			$lng_delta_m = ($lng[$ix] - $lng_avg) * $DEG2M;

			printf "%.2f\t%.2f\t%.2f\t%.9f\t%.9f\t%.2f\t%s\n",
				$lat_delta_m, $lng_delta_m,
				sqrt($lat_delta_m*$lat_delta_m + $lng_delta_m*$lng_delta_m),
				$lat[$ix], $lng[$ix], $hdop[$ix],
				$hdop_class;
		}

	}
	print "e\n";	#end of plot set
}

}#MAIN

#eof errorplot.pl
# vim:ts=4:sw=4:ai