This video covers a first pass at converted a NMEA log from an Airmar weather station into a plot of the wander of the GPS position. The python code is not what I would want to have in the long run, but it is similar to what I might write as a first pass to make sure I know how parse the data and that my results make sense.
The notes for this video are not very polished. You might want to take a look at the bottom section that is a dump of the history from ipython after I produced the video.
wget http://vislab-ccom.unh.edu/~schwehr/Classes/2011/esci895-researchtools/examples/ccom-airmar-2011-10-11.bz2
bunzip2 ccom-airmar-2011-10-11.bz2wc -l ccom-airmar-2011-10-11 881046 ccom-airmar-2011-10-11
head ccom-airmar-2011-10-11 # START LOGGING UTC seconds since the epoch: 1318291200.14 # SPEED: 4800 # PORT: /dev/ttyS0 # TIMEOUT: 300.0 # STATIONID: ccom-airmar # DAEMON MODE: False $PNTZNT,1318291200.14,127.0.0.1,148.167.132.201,3,1318289568.2,0.000034,-20,0.086365,0.052994*12,ccom-airmar,1318291200.14 $GPVTG,277.9,T,293.3,M,0.1,N,0.1,K,D*26,rccom-airmar,1318291200.14 $WIMWV,307.3,R,0.7,N,A*23,rccom-airmar,1318291200.17 $GPZDA,000000,11,10,2011,00,00*4B,rccom-airmar,1318291200.27
What types of NMEA strings do we have?
grep '^\$' ccom-airmar-2011-10-11 | cut -d, -f1 | sort -u $GPGGA $GPVTG $GPZDA $HCHDT $PNTZNT $WIMDA $WIMWD $WIMWV
We can look up the NMEA sentences here:
http://gpsd.berlios.de/NMEA.html
=== GGA - Global Positioning System Fix Data ===
Time, Position and fix related data for a GPS receiver.
------------------------------------------------------------------------------
11
1 2 3 4 5 6 7 8 9 10 | 12 13 14 15
| | | | | | | | | | | | | | |
$--GGA,hhmmss.ss,llll.ll,a,yyyyy.yy,a,x,xx,x.x,x.x,M,x.x,M,x.x,xxxx*hh<CR><LF>
------------------------------------------------------------------------------
Field Number:
1. Universal Time Coordinated (UTC)
2. Latitude
3. N or S (North or South)
4. Longitude
5. E or W (East or West)
6. GPS Quality Indicator,
- 0 - fix not available,
- 1 - GPS fix,
- 2 - Differential GPS fix
(values above 2 are 2.3 features)
- 3 = PPS fix
- 4 = Real Time Kinematic
- 5 = Float RTK
- 6 = estimated (dead reckoning)
- 7 = Manual input mode
- 8 = Simulation mode
7. Number of satellites in view, 00 - 12
8. Horizontal Dilution of precision (meters)
9. Antenna Altitude above/below mean-sea-level (geoid) (in meters)
10. Units of antenna altitude, meters
11. Geoidal separation, the difference between the WGS-84 earth
ellipsoid and mean-sea-level (geoid), "-" means mean-sea-level
below ellipsoid
12. Units of geoidal separation, meters
13. Age of differential GPS data, time in seconds since last SC104
type 1 or 9 update, null field when DGPS is not used
14. Differential reference station ID, 0000-1023
15. Checksum
One NMEA GPS location message looks like this:
$GPGGA,000000,4308.1250,N,07056.3750,W,2,9,1.1,35.7,M,,,,*04
Start ipython with ipython –pylab
!head -20000 ccom-airmar-2011-10-11 | tail
$HCHDT,26.2,T*1F,rccom-airmar,1318293159.77
$WIMWV,214.5,T,0.3,N,A*24,rccom-airmar,1318293159.84
$WIMWV,218.9,R,0.3,N,A*22,rccom-airmar,1318293159.89
$HCHDT,26.1,T*1C,rccom-airmar,1318293160.17
$GPVTG,268.5,T,283.9,M,0.1,N,0.1,K,D*2F,rccom-airmar,1318293160.25
$WIMWV,211.9,R,0.4,N,A*2C,rccom-airmar,1318293160.3
$GPZDA,003240,11,10,2011,00,00*4E,rccom-airmar,1318293160.39
$WIMDA,30.1915,I,1.0224,B,17.6,C,,,,,,,229.0,T,244.4,M,0.3,N,0.2,M*2E,rccom-airmar,1318293160.52
$GPGGA,003240,4308.1263,N,07056.3750,W,2,9,1.1,41.3,M,,,,*06,rccom-airmar,1318293160.65
$WIMWD,228.9,T,244.3,M,0.3,N,0.2,M*5B,rccom-airmar,1318293160.74
line = '$GPGGA,003240,4308.1263,N,07056.3750,W,2,9,1.1,41.3,M,,,,*06,rccom-airmar,1318293160.65'
'GGA' in line
'VTG' in line#!/usr/bin/env python
num_gga = 0
for line in open('ccom-airmar-2011-10-11'):
if 'GGA' not in line:
# Skip all lines that are not GPS position messages
continue
num_gga += 1
print num_ggaline = '$GPGGA,003240,4308.1263,N,07056.3750,W,2,9,1.1,41.3,M,,,,*06,rccom-airmar,1318293160.65'
fields = line.split(',')
y = int(fields[2][:2]) + float(fields[2][2:])/ 100.
if fields[3] == 'S':
y = -y
x = int(fields[4][:3]) + float(fields[4][3:]) / 100.
if fields[5] == 'W':
x = -x
return '%f, %f' % (x,y)#!/usr/bin/env python
num_gga = 0
X = []
Y = []
for line_num, line in enumerate(open('ccom-airmar-2011-10-11')):
if line_num > 100:
break
if 'GGA' not in line:
# Skip all lines that are not GPS position messages
continue
fields = line.split(',')
y = int(fields[2][:2]) + float(fields[2][2:])/ 100.
if fields[3] == 'S':
y = -y
x = int(fields[4][:3]) + float(fields[4][3:]) / 100.
if fields[5] == 'W':
x = -x
#print x,y
X.append(x)
Y.append(y)#!/usr/bin/env python
def load_gga(filename):
num_gga = 0
X = []
Y = []
for line_num, line in enumerate(open(filename)):
#if line_num > 100:
# break
if 'GGA' not in line:
# Skip all lines that are not GPS position messages
continue
fields = line.split(',')
y = int(fields[2][:2]) + float(fields[2][2:])/ 100.
if fields[3] == 'S':
y = -y
x = int(fields[4][:3]) + float(fields[4][3:]) / 100.
if fields[5] == 'W':
x = -x
X.append(x)
Y.append(y)
return X, Yx,y = wx2.load_gga('ccom-airmar-2011-10-11')
plotls -F -l
run wx1
head ccom-airmar-2011-10-11
head -15 ccom-airmar-2011-10-11
line = '$GPGGA,000000,4308.1250,N,07056.3750,W,2,9,1.1,35.7,M,,,,*04,rccom-airmar,1318291200.53'
line
'GGA' in line
'MWD' in line
junk = '$WIMWD,338.0,T,353.4,M,0.6,N,0.3,M*56,rccom-airmar,1318291200.61'
'GGA' in junk
'MWD' in junk
'GGA' in line
'GGA' not in line
'GGA' not in junk
line
line = '$GPGGA,000000,4308.1250,N,07056.3750,W,2,9,1.1,35.7,M,,,,*04,rccom-airmar,1318291200.53'
line.split(',')
fields[2]
fields = line.split(',')
whos
# Latitude
fields[2]
fields[2][:2]
int(fields[2][:2])
fields[2][2:]
float(fields[2][2:])
float(fields[2][2:])/60.
int(fields[2][:2]) + float(fields[2][2:])/60.
fields
run wx1.py
# Longitude
fields
fields[4]
fields[4][:3]
int(fields[4][:3])
fields[4][3:]
float(fields[4][3:])
float(fields[4][3:]) / 60.
int(fields[4][:3]) + float(fields[4][3:]) / 60.
run wx1.py
import wx2
?wx2.load_gga
wx2.load_gga('ccom-airmar-2011-10-11')
_ip.magic("whos")
reload(wx2)
wx2.load_gga('ccom-airmar-2011-10-11')
x,y = wx2.load_gga('ccom-airmar-2011-10-11')
len(x)
print x
plot(x,y)
average(x)
average(y)file:./figures/v13-gps-wander.png