#!/usr/bin/env python """rowe-adcp-sink.py A simple utility to start the Rowe ADCP pinging and record the raw data it streams back over the serial port. TODO: Make this a more OO instrument interface, to simplify passing parameters and file handles around. """ import os import serial import struct import collections import time import crcmod from cStringIO import StringIO # TODO: migrate to io for python3 from scipy.io import loadmat, savemat import datetime import netCDF4 import numpy as np import binascii kermit = crcmod.predefined.mkPredefinedCrcFun('kermit') ccitt_false = crcmod.predefined.mkPredefinedCrcFun('crc-ccitt-false') ccitt_aug = crcmod.predefined.mkPredefinedCrcFun('crc-aug-ccitt') 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 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 vblist = [('b0', np.float32), ('b1', np.float32), ('b2', np.float32),('b3', np.float32)] vectorB_t = rg.createCompoundType(np.dtype(vblist),'vector_beam') vilist = [('u', np.float32), ('v', np.float32), ('w', np.float32),('e', np.float32)] vectorI_t = rg.createCompoundType(np.dtype(vilist),'vector_instrument') velist = [('N', np.float32), ('E', np.float32), ('D', np.float32),('e', np.float32)] vectorE_t = rg.createCompoundType(np.dtype(velist),'vector_earth') vbilist = [('b0', np.int16), ('b1', np.int16), ('b2', np.int16),('b3', np.int16)] vectorBint_t = rg.createCompoundType(np.dtype(vbilist),'int_vector_beam') viilist = [('u', np.int16), ('v', np.int16), ('w', np.int16),('e', np.int16)] vectorIint_t = rg.createCompoundType(np.dtype(viilist),'int_vector_instrument') veilist = [('N', np.int16), ('E', np.int16), ('D', np.int16),('e', np.int16)] vectorEint_t = rg.createCompoundType(np.dtype(veilist),'int_vector_earth') 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' ensemble_number_ = rg.createVariable('ensemble_number', np.int16, 'time') ensemble_number_.units = 'count' ed = rg.createGroup('ensemble_data') en = ed.createVariable('ensemble_number', np.int16, 'time') nbins = ed.createVariable('number_of_bins', np.int16, 'time') nbeams = ed.createVariable('number_of_beams', np.int16, 'time') ndp = ed.createVariable('number_of_desired_pings', np.int16, 'time') nap = ed.createVariable('number_of_acquired_pings', np.int16, 'time') for v in [ed, nbins, nbeams, ndp, nap]: v.units = 'count' status = ed.createVariable('status', np.int16, 'time') # TODO: this might be OK as a bool status.units = 'none' instrument_time = ed.createVariable('instrument_time', np.float64, 'time') instrument_time.units = 'seconds since 1970-00-00T00:00:00.0Z' instrument_time.calendar = 'standard' instrument_epoch = ed.createVariable('epoch', np.float64, 'time') instrument_epoch.units = 'seconds since 1970-00-00T00:00:00.0Z' bt = rg.createGroup('bottom_track') heading = bt.createVariable('heading', np.float32, 'time') pitch = bt.createVariable('pitch', np.float32, 'time') roll = bt.createVariable('roll', np.float32, 'time') for v in [heading, pitch, roll]: v.units = 'angular_degree' wt = bt.createVariable('water_temperature', np.float32, 'time') st = bt.createVariable('system_temperature', np.float32, 'time') for v in [bt, wt]: v.units = 'degree_celsius' dep = bt.createVariable('depth', np.float32, 'time') dep.units = 'meter' pres = bt.createVariable('pressure', np.float32, 'time') pres.units = 'bar' sal = bt.createVariable('salinity', np.float32, 'time') sal.units = 'part_per_thousand' ss = bt.createVariable('soundspeed', np.float32, 'time') ss.units = 'meter_per_second' so = bt.createVariable('status_or', np.int16, 'time') so.units = 'none' bt_nob = bt.createVariable('number_of_beams', np.int16, 'time') pc = bt.createVariable('ping_count', np.int16, 'time') for v in [bt_nob, pc]: v.units = 'count' vr = bt.createVariable('vertical_range', vectorB_t, 'time') vr.units = 'meter' snr = bt.createVariable('signal_to_noise_ratio', vectorB_t, 'time') snr.units = 'none' amp = bt.createVariable('amplitude', vectorB_t, 'time') amp.units = 'decibel' corr = bt.createVariable('correlation', vectorB_t, 'time') corr.units = 'none' bvel = bt.createVariable('beam_velocity', vectorB_t, 'time') pavg = bt.createVariable('pings_averaged', vectorBint_t, 'time') ivel = bt.createVariable('instrument_velocity', vectorI_t, 'time') isa = bt.createVariable('instrument_solutions_averaged', vectorIint_t, 'time') evel = bt.createVariable('earth_velocity', vectorE_t, 'time') esa = bt.createVariable('earth_solutions_averaged', vectorEint_t, 'time') for v in [bvel, ivel, evel]: v.units = 'meter_per_second' for v in [pavg, isa, esa]: v.units = 'count' return rg def write_ensemble_to_nc4(timestamp, ensemble_number, payload, nf, verbosity=1): 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) nf.variables['ensemble_number'][k] = ensemble_number ed = nf.groups['ensemble_data'] bt = nf.groups['bottom_track'] mat = loadmat(StringIO(payload)) edm = mat['E000008'].ravel() # only ravel 1D arrays btm = mat['E000010'].ravel() # only ravel 1D arrays ed.variables['ensemble_number'][k] = edm[0] ed.variables['number_of_bins'][k] = edm[1] ed.variables['number_of_beams'][k] = edm[2] ed.variables['number_of_desired_pings'][k] = edm[3] ed.variables['number_of_acquired_pings'][k] = edm[4] ed.variables['status'][k] = edm[5] its = datetime.datetime(edm[6], edm[7], edm[8], edm[9], edm[10], edm[11]) its += datetime.timedelta(seconds = 0.01 * edm[12]) ed.variables['instrument_time'][k] = netCDF4.date2num(its, units=ed.variables['instrument_time'].units, calendar=ed.variables['instrument_time'].calendar) ed.variables['epoch'][k] = unix_epoch_seconds(its) # TODO: what are in 0 & 1? bt.variables['heading'][k] = btm[2] bt.variables['pitch'][k] = btm[3] bt.variables['roll'][k] = btm[4] bt.variables['water_temperature'][k] = btm[5] bt.variables['system_temperature'][k] = btm[6] bt.variables['depth'][k] = btm[7] bt.variables['pressure'][k] = btm[8] bt.variables['salinity'][k] = btm[9] bt.variables['soundspeed'][k] = btm[10] bt.variables['status_or'][k] = btm[11] bt.variables['number_of_beams'][k] = btm[12] bt.variables['ping_count'][k] = btm[13] vr = bt.variables['vertical_range'] vr[k] = np.array((btm[14], btm[15], btm[16], btm[17]), vr.dtype) snr = bt.variables['signal_to_noise_ratio'] snr[k] = np.array((btm[18], btm[19], btm[20], btm[21]), snr.dtype) amp = bt.variables['amplitude'] amp[k] = np.array((btm[22], btm[23], btm[24], btm[25]), amp.dtype) corr = bt.variables['correlation'] corr[k] = np.array((btm[26], btm[27], btm[28], btm[29]), corr.dtype) bv = bt.variables['beam_velocity'] bv[k] = np.array((btm[30], btm[31], btm[32], btm[33]), bv.dtype) pa = bt.variables['pings_averaged'] pa[k] = np.array((btm[34], btm[35], btm[36], btm[37]), pa.dtype) iv = bt.variables['instrument_velocity'] iv[k] = np.array((btm[38], btm[39], btm[40], btm[41]), iv.dtype) isa = bt.variables['instrument_solutions_averaged'] isa[k] = np.array((btm[42], btm[43], btm[44], btm[45]), isa.dtype) ev = bt.variables['earth_velocity'] ev[k] = np.array((btm[46], btm[47], btm[48], btm[49]), ev.dtype) esa = bt.variables['earth_solutions_averaged'] esa[k] = np.array((btm[50], btm[51], btm[52], btm[53]), esa.dtype) if verbosity > 0: msg = '{0} vertical range: {1:0.1f}, {2:0.1f}, {3:0.1f}, {4:0.1f}; instrument velocity: {5:0.2f}, {6:0.2f}, {7:0.2f}, {8:0.2f}' print(msg.format(timestamp, vr[k]['b0'], vr[k]['b1'], vr[k]['b2'], vr[k]['b3'], iv[k]['u'], iv[k]['v'], iv[k]['w'], iv[k]['e'])) def configure(sio, config_filename, lf): with open(config_filename, 'r') as cf: for line in cf.readlines(): sio.write(line.rstrip() + '\r') sio.write('CSHOW\r') sio.write('SLEEPSECONDS\r') sio.write('CSAVE\r') time.sleep(0.1) # give it a bit to start responding while sio.inWaiting() > 0: response = sio.readline().rstrip() # grab the return from CSHOW, etc lf.write(response) print response time.sleep(1e-2) def set_time(sio, lf): now = datetime.datetime.utcnow() tstr = now.strftime('%Y/%m/%d,%H:%M:%S') sio.write('STIME ' + now.strftime('%Y/%m/%d,%H:%M:%S') + '\r') ack = sio.readline() # this is an acknowledgement and an echo of the STIME set command lf.write(ack) print 'host time: {0}'.format(datetime.datetime.utcnow()) sio.write('STIME\r') ack = sio.readline() # this is the acknowledgment and an echo of the STIME command lf.write(ack) ack = sio.readline() # this is the acknowledgment and an echo of the STIME command lf.write(ack) stime = sio.readline() lf.write(stime) print 'sensor time: {0}'.format(stime) def stop(sio, n, s): for i in range(n): sio.write('STOP\r') sio.flush() # do it now, don't wait with the buffer time.sleep(s) def reset(sio, n=2, s = 0.1): for i in range(n): sio.sendBreak() time.sleep(s) sio.flushInput() sio.flushOutput() sio.sendBreak() def read_starting_with(sio, blob, nbytes = 0): """Return nbytes, starting with specified blob.""" blob_length = len(blob) if nbytes < blob_length: nbytes = blob_length q = collections.deque(sio.read(nbytes), maxlen=len(blob)) while ''.join(q) != blob: q.append(sio.read(1)) return blob + sio.read(nbytes - blob_length) def run(sio, period=90, ef = None, nf = None, lf = None, verbosity = 1): if type(period) in (int, float): period = datetime.timedelta(seconds=period) try_number = 0 sio.write('STIME\r') sio.write('START\r') sio.flush() # do it now, don't wait with the buffer tstart = datetime.datetime.utcnow() stime = sio.readline() # grab the STIME, the timestamp tstr = sio.readline() start = sio.read(6) # grab the START response = [stime, tstr, start] for line in response: lf.write(line) print line.lstrip().rstrip() while datetime.datetime.utcnow() < tstart + period: if sio.inWaiting() > 31: chunk = read_starting_with(sio, 16 * "\x80", 32) ensemble_number, payload_size, preamble_position = parse_header(chunk) if (ensemble_number > 0) and (payload_size > 0): timestamp = datetime.datetime.utcnow() ef.write(chunk) payload = sio.read(payload_size) ef.write(payload) checksum = sio.read(4) ef.write(checksum) msg = '{0} read ensemble # {1} with {2}-byte payload; checksum: {3} {4}, kermit: {5}, crc-ccitt-false: {6}, crc-aug-ccitt: {7}\n' cksum_msg = msg.format(timestamp, ensemble_number, payload_size, binascii.hexlify(checksum), struct.unpack('hh', checksum), kermit(payload), ccitt_false(payload), ccitt_aug(payload), #binascii.hexlify(kermit(payload)), #binascii.hexlify(ccitt_false(payload)), #binascii.hexlify(ccitt_aug(payload)) ) if verbosity > 1: print cksum_msg.rstrip() lf.write(cksum_msg) write_ensemble_to_nc4(timestamp, ensemble_number, payload, nf, verbosity) if verbosity > 1: inspect_MATv4(payload) try_number = 0 elif preamble_position > 0: print 'try number {0} with preamble at position {1}'.format(try_number, preamble_position) else: try_number += 1 print '...finished recording...' stop(sio, 9, 0.25) sio.flushInput() sio.flushOutput() stop(sio, 3, 0.25) sio.write('STIME\r') while sio.inWaiting() > 0: response = sio.readline() print response.lstrip().rstrip() lf.write(response) def parse_header(blob, verbosity = 0): """Look for the Rowe header in the blob -- parse if found.""" preamble = 16 * "\x80" ensemble_number = -1 payload_size = -1 preamble_position = -1 preamble_position = blob.find(preamble) if preamble_position == 0: if verbosity > 1: print 'blob starts with preamble' i = struct.unpack(8*'i', blob) ensemble_number = i[4] payload_size = i[6] elif preamble_position != -1: print 'found preamble at position {0}'.format(preamble_position) else: print 'did not find preamble in blob', blob return ensemble_number, payload_size, preamble_position def inspect_MATv4(payload): for k, v in loadmat(StringIO(payload)).iteritems(): print k, v.shape def convert_ens_to_nc4(ens_filename=None, nc_filename=None, verbosity=0): if nc_filename is None: nc_filename = '.'.join((ens_filename, 'nc4')) nf = setup_netCDFv4(nc_filename) with open(ens_filename, 'rb') as ens: # TODO: probably a more efficient implementation with generator expressions raw = ens.read() # reads entire file. There's probably a smarter option. payloads = [ensemble[16:] for ensemble in raw.split('\x80'*16)[1:]] for ensemble_number, payload in enumerate(payloads): write_ensemble_to_nc4(datetime.datetime.utcnow(), ensemble_number, payload, nf, verbosity) # TODO: handle None as an arg for timestamp, use the device timestamp nf.close() def main(port='/dev/ttyUSB0', baudrate=115200, timeout=1.0, period = 10, config_filename='rowe-adcp.conf', output_path = '/tmp', ens_filename=None, nc_filename=None, log_filename=None, verbosity = 0): tstart = datetime.datetime.utcnow() if type(period) in (int, float): period = datetime.timedelta(seconds=period) if ens_filename is None: ens_filename = '.'.join((tstart.strftime('%Y%m%dT%H%M%S'), 'rowe-adcp', 'ens')) if nc_filename is None: nc_filename = '.'.join((tstart.strftime('%Y%m%dT%H%M%S'), 'rowe-adcp', 'nc4')) if log_filename is None: log_filename = '.'.join((tstart.strftime('%Y%m%dT%H%M%S'), 'rowe-adcp', 'log')) ef = open(os.path.join(output_path, ens_filename), 'wb') nf = setup_netCDFv4(os.path.join(output_path, nc_filename)) lf = open(os.path.join(output_path, log_filename), 'w', 1) # TODO: Use logger library instead. 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) reset(sio) banner = [sio.readline() for i in range(6)] # try reading the banner for line in banner: print (line.lstrip().rstrip()) lf.writelines(banner) configure(sio, config_filename, lf) set_time(sio, lf) run(sio, period, ef, nf, lf) sio.close() lf.close() nf.close() ef.close() if __name__ == "__main__": import argparse program_description = """read data from Rowe ADCP connected via 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 communicate on') parser.add_argument('-b', '--baudrate', default=115200, type=int, help='baud rate to communicate 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('-c', '--config-filename', default='rowe-adcp.conf', help='path to output directory') parser.add_argument('-o', '--output-path', default='/tmp', help='path to output directory') parser.add_argument('-e', '--ens-filename', default=None, help='name of output file') 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', action='count', help='increase verbosity of output (does not affect logs)') parser.add_argument('-C','--convert', action='store_true', default=False, help='convert an existing ens file into a netCDF4 file') # TODO: make v flag work like normal instead of requiring an arg args = parser.parse_args() if args.convert: if args.ens_filename is not None: convert_ens_to_nc4(args.ens_filename, args.nc_filename, args.verbosity) else: print('Cannot convert ens file to nc4 without ens_filename input.') else: main(args.port, args.baudrate, args.timeout, args.period, args.config_filename, args.output_path, args.ens_filename, args.nc_filename, args.log_filename, args.verbosity)