#!/usr/bin/env python """ note: I considered using `pynmea`, but since we're only dealing with GGA sentences and `pynmea` doesn't actually handle the ddmm.mm to dd.ddd (i.e., degree minute without separator to decimal degree) conversion, it doesn't really add anything extra over a simple string parser. There is probably a simpler or cleaner way to do this. In the future, this could be extended to handle other NMEA strings and sort them into different groups, probably by string type. """ import os import serial import datetime import warnings import numpy as np import netCDF4 def unix_epoch_seconds(t): """Simple function to convert datetime to epoch seconds. Providing this to handle cases where the python and datetime versions are old, and as a workaround to get sub-second time resolution in the netCDF4 file. """ delta = t - datetime.datetime.utcfromtimestamp(0) try: return delta.total_seconds() except AttributeError: # probably an older version of python and datetime return delta.days * 86400 + delta.seconds + delta.microseconds * 1e-6 def float_or_nan(a): try: return float(a) except: return np.nan def parse_gpgga(gga): """Parse a single GPGGA string into numbers. Since we are not concerned with a series of GPGGA strings here, it is probably most straigtforward to use sting methods. """ ccs = nmea_checksum(gga) elements = gga.split(',') sentence_identifier = elements[0] if sentence_identifier != '$GPGGA': print 'string' + gga + 'is not $GPGGA' tstr = elements[1] time_of_day = datetime.time(int(tstr[:2]), int(tstr[2:4]), int(tstr[4:])) latitude = degrees_minutes_to_decimal_degrees(elements[2]) north_latitude = (elements[3] == 'N')*2 - 1 latitude *= north_latitude longitude = degrees_minutes_to_decimal_degrees(elements[4]) east_longitude = (elements[5] == 'E')*2 - 1 longitude *= east_longitude return dict(time_of_day = time_of_day, latitude = latitude, longitude=longitude, fix_quality = int(elements[6]), number_of_satellites = int(elements[7]), horizontal_dilution_of_precision = float_or_nan(elements[8]), altitude = float_or_nan(elements[9]), altitude_units = elements[10], height_above_geoid = float_or_nan(elements[11]), height_above_geoid_units = elements[12], time_since_last_DGPS_update = float_or_nan(elements[13]), checksum = int(elements[14][1:3],16), calculated_checksum = ccs) def nmea_checksum(s): #print s frame, cs = s[1:].split('*') ccs = 0 for c in frame: ccs ^= ord(c) if int(cs.rstrip(), 16) is not ccs: warnings.warn('NMEA checksum mismatch! read: {0}, calculated: {1}'.format(hex(cs), hex(ccs))) return ccs def degrees_minutes_to_decimal_degrees(dm): """Convert degrees, minutes, and hemisphere to decimal degrees. """ boundary = dm.find('.') - 2 return float(dm[:boundary]) + float(dm[boundary:]) / 60. def setup_netCDFv4(filename): """Set up groups, metadata, and dimensions for a netCDFv4 log. """ rg = netCDF4.Dataset(filename, 'w', clobber=False, format='NETCDF4') # TODO write nc4 metadata time = rg.createDimension('time', None) # set up time as an unlimited dimension ts = rg.createVariable('time', np.float64, 'time') ts.units = 'seconds since 1970-00-00T00:00:00.0Z' ts.calendar = 'standard' # e.g., ts[0] = netCDF4.date2num(t[0], units=ts.units, calendar=ts.calendar) epoch = rg.createVariable('epoch', np.float64, 'time') epoch.units = 'seconds since 1970-00-00T00:00:00.0Z' it = rg.createVariable('gps_time', np.float64, 'time') it.units = 'seconds since 1970-00-00T00:00:00.0Z' it.calendar = 'standard' ie = rg.createVariable('gps_epoch', np.float64, 'time') ie.units = 'seconds since 1970-00-00T00:00:00.0Z' lat = rg.createVariable('latitude', np.float64, 'time') lat.units = 'degrees_latitude' lon = rg.createVariable('longitude', np.float64, 'time') lon.units = 'degrees_longitude' fq = rg.createVariable('fix_quality', np.int32, 'time') nos = rg.createVariable('number_of_satellites', np.int32, 'time') for v in (fq, nos): v.units = 'count' hdop = rg.createVariable('horizontal_dilution_of_precision', np.float64, 'time') alt = rg.createVariable('altitude', np.float64, 'time') hag = rg.createVariable('height_above_geoid', np.float64, 'time') for v in (alt, hag): v.units = 'm' # guess tslDu = rg.createVariable('time_since_last_DGPS_update', np.float64, 'time') cs = rg.createVariable('checksum', np.uint16, 'time') ccs = rg.createVariable('calculated_checksum', np.uint16, 'time') gcs = rg.createVariable('good_checksum', np.uint16, 'time') return rg def write_fix_to_nc4(timestamp, gpgga, nf): k = nf.variables['time'].shape[0] # get the index of the next timestep nf.variables['time'][k] = netCDF4.date2num(timestamp, units=nf.variables['time'].units, calendar=nf.variables['time'].calendar) nf.variables['epoch'][k] = unix_epoch_seconds(timestamp) raw = parse_gpgga(gpgga) gps_timestamp = datetime.datetime.combine(timestamp.date(), raw['time_of_day']) # TODO: This probably doesn't roll over well at midnight. nf.variables['gps_time'][k] = netCDF4.date2num(gps_timestamp, units=nf.variables['gps_time'].units, calendar=nf.variables['gps_time'].calendar) nf.variables['gps_epoch'][k] = unix_epoch_seconds(gps_timestamp) nf.variables['latitude'][k] = raw['latitude'] nf.variables['longitude'][k] = raw['longitude'] nf.variables['fix_quality'][k] = raw['fix_quality'] nf.variables['number_of_satellites'][k] = raw['number_of_satellites'] nf.variables['horizontal_dilution_of_precision'][k] = raw['horizontal_dilution_of_precision'] nf.variables['altitude'][k] = raw['altitude'] nf.variables['height_above_geoid'][k] = raw['height_above_geoid'] nf.variables['time_since_last_DGPS_update'][k] = raw['time_since_last_DGPS_update'] nf.variables['checksum'][k] = raw['checksum'] nf.variables['calculated_checksum'][k] = raw['calculated_checksum'] nf.variables['good_checksum'][k] = (raw['checksum'] == raw['calculated_checksum']) def main(port = '/dev/ttyUSB0', baudrate = 9600, timeout = 0.1, period = 10, output_path='/tmp', nc_filename = None, log_filename = None, verbosity = 0): tstart = datetime.datetime.utcnow() if type(period) in (int, float): period = datetime.timedelta(seconds=period) if log_filename is None: log_filename = '.'.join((tstart.strftime('%Y%m%dT%H%M%S'), 'gpgga', 'log')) if nc_filename is None: nc_filename = '.'.join((tstart.strftime('%Y%m%dT%H%M%S'), 'gpgga', 'nc4')) lf = open(os.path.join(output_path,log_filename), 'w', 1) nf = setup_netCDFv4(os.path.join(output_path, nc_filename)) siokw = dict(bytesize=8, parity='N', stopbits=1, xonxoff=0, rtscts=0, writeTimeout=None, dsrdtr=None) sio = serial.Serial(port=port, baudrate=baudrate, timeout=timeout, **siokw) while datetime.datetime.utcnow() < tstart + period: line = sio.readline().rstrip() ts = datetime.datetime.utcnow() if line.startswith('$GPGGA'): lf.write('{0} {1}\n'.format(ts, line)) write_fix_to_nc4(ts, line, nf) else: if verbosity > 0: print 'received non-$GPGGA line', line print('{0} {1}'.format(ts, line)) sio.close() nf.close() lf.close() if __name__ == "__main__": import argparse program_description = """read $GPGGA strings from serial and log to nc4""" parser = argparse.ArgumentParser(description=program_description) parser.add_argument('-V', '--version', action='version', version='%(prog)s 0.0.1', help='display version information and exit') parser.add_argument('-p', '--port', default='/dev/ttyUSB0', help='serial port to listen on') parser.add_argument('-b', '--baudrate', default=9600, type=int, help='baud rate to listen with') parser.add_argument('-t', '--timeout', default=1, type=float, help='timeout for serial port') parser.add_argument('-P', '--period', default=10, type=float, help='number of seconds to listen') parser.add_argument('-o', '--output-path', default='/tmp', help='path to output directory') parser.add_argument('-n', '--nc_filename', default=None, help='name of output file') parser.add_argument('-l', '--log_filename', default=None, help='name of output file') parser.add_argument('-v','--verbosity', default=0, type=int, help='verbosity to use with the console (does not affect logs)') # TODO: make v flag work like normal instead of requiring an arg args = parser.parse_args() main(args.port, args.baudrate, args.timeout, args.period, args.output_path, args.nc_filename, args.log_filename, args.verbosity)