#!/usr/bin/env python
#
__author__    = 'Mike McCann'
__version__ = '1.0'
__date__ = '2013-07-18'
__copyright__ = '2013'
__license__   = 'GPL v3'
__contact__   = 'mccann at mbari.org'

__doc__ = '''

Parse logged data from BEDs .EVT or .WAT file to compute translations in x, y, z and rotations
about these axes.  Output as NetCDF with the intention of of loading into a STOQS database.
The resulting Netcdf file will be of featureType timeSeries for stationary events (the default).
If the --trajectory option is given then a NetCDF file with featureType of trajectory will be
created using the path in the file specified with the --trajectory option.
'''

import os
import sys
import csv
import math
import coards
import datetime
import numpy as np
from BEDS import BEDS
from seawater import eos80
from netCDF4 import Dataset


class BEDS_NetCDF(BEDS):

    def __init__(self):
        '''
        Initialize with options
        '''

        return super(BEDS_NetCDF, self).__init__()

    def createNetCDFfromFile(self):
        '''
        Read data from EVT or WAT file and apply operations to convert it to data with units that
        are then written to a NetCDF file.
        '''

        # Suppress 'util.py:70: RuntimeWarning: invalid value encountered in sqrt'
        np.seterr(all='ignore')

        for fileName in self.inputFileNames:
            # Make sure input file is openable
            print 'Input fileName = ', fileName
            try:
                with open(fileName): 
                    pass
            except IOError:
                raise Exception('Cannot open input file %s' % fileName)

            if self.sensorType == 'Invensense':
                self.readBEDsFile(fileName)
                self.processAccelerations()
                self.processRotations()
            else:
                raise Exception("No handler for sensorType = %s" % self.sensorType)

        if self.args.output:
            self.outFile = self.args.output
        elif len(self.inputFileNames) == 1:
            self.outFile = '.'.join(self.inputFileNames[0].split('.')[:-1]) + '.nc'
        else:
            raise Exception("Must specify --output if more than one input file.")

        # Interpolate data to regularly spaced time values - may need to do this to improve accuracy
        # (See http://www.freescale.com/files/sensors/doc/app_note/AN3397.pdf)
        ##si = linspace(self.s2013[0], self.s2013[-1], len(self.s2013))
        ##axi = interp(si, self.s2013, self.ax)

        # TODO: Review this calculation - may need to rotate these to the absolute (not rotating) frame
        # Double integrate accelerations to get position and construct X3D position values string
        # (May need to high-pass filter the data to remove noise that can give unreasonably large positions.)
        ##x = self.cumtrapz(self.s2013, self.cumtrapz(self.s2013, self.ax))
        ##y = self.cumtrapz(self.s2013, self.cumtrapz(self.s2013, self.ay))
        ##z = self.cumtrapz(self.s2013, self.cumtrapz(self.s2013, self.az))

        dateCreated = datetime.datetime.now().strftime("%d %B %Y %H:%M:%S")
        yearCreated = datetime.datetime.now().strftime("%Y")

        # Create the NetCDF file
        self.ncFile = Dataset(self.outFile, 'w')

        # Time dimensions for both trajectory and timeSeries datasets - IMU and pressure have different times
        self.ncFile.createDimension('time', len(self.s2013))
        self.time = self.ncFile.createVariable('time', 'float64', ('time',))
        self.time.standard_name = 'time'
        self.time.units = 'seconds since 2013-01-01'
        self.time[:] = self.s2013

        if self.featureType == 'timeseries':
            # Save with COARDS compliant station coordinates
            self.ncFile.createDimension('ptime', len(self.ps2013))
            self.ptime = self.ncFile.createVariable('ptime', 'float64', ('ptime',))
            self.ptime.standard_name = 'time'
            self.ptime.units = 'seconds since 2013-01-01'
            self.ptime[:] = self.ps2013

            self.ncFile.createDimension('latitude', 1)
            self.latitude = self.ncFile.createVariable('latitude', 'float64', ('latitude',))
            self.latitude.long_name = 'LATITUDE'
            self.latitude.standard_name = 'latitude'
            self.latitude.units = 'degree_north'
            self.latitude[0] = self.lat
    
            self.ncFile.createDimension('longitude', 1)
            self.longitude = self.ncFile.createVariable('longitude', 'float64', ('longitude',))
            self.longitude.long_name = 'LONGITUDE'
            self.longitude.standard_name = 'longitude'
            self.longitude.units = 'degree_east'
            self.longitude[0] = self.lon
    
            self.ncFile.createDimension('depth', 1)
            self.depth = self.ncFile.createVariable('depth', 'float64', ('depth',))
            self.depth.long_name = 'depth'
            self.depth.standard_name = 'depth'
            self.depth.units = 'm'
            self.depth[0] = self.dep

            # Record Variable - Pressure and Depth
            pr = self.ncFile.createVariable('PRESS', 'float64', ('ptime', 'depth', 'latitude', 'longitude'))
            pr.long_name = 'External Instrument Pressure'
            pr.coordinates = 'ptime depth latitude longitude'
            pr.units = 'bar'
            pr[:] = self.pr.reshape(len(self.pr), 1, 1, 1)

            bd = self.ncFile.createVariable('BED_DEPTH', 'float64', ('ptime', 'depth', 'latitude', 'longitude'))
            bd.long_name = 'Instrument Depth'
            bd.coordinates = 'ptime depth latitude longitude'
            bd.units = 'm'
            bd[:] = eos80.dpth(self.pr * 10.0, self.lat).reshape(len(self.pr), 1, 1, 1)    # Convert pressure to depth

            bdi = self.ncFile.createVariable('BED_DEPTH_INTERP', 'float64', ('time', 'depth', 'latitude', 'longitude'))
            bdi.long_name = 'Instrument Depth'
            bdi.coordinates = 'ptime depth latitude longitude'
            bdi.units = 'm'
            bdi[:] = np.interp(self.s2013, self.ps2013, bd[:].reshape(len(self.pr))).reshape(len(self.s2013), 1, 1, 1)

            # Record Variables - Accelerations
            xa = self.ncFile.createVariable('XA', 'float64', ('time', 'depth', 'latitude', 'longitude'))
            xa.long_name = 'Acceleration along X-axis'
            xa.coordinates = 'time depth latitude longitude'
            xa.units = 'g'
            xa[:] = self.ax.reshape(len(self.ax), 1, 1, 1)

            ya = self.ncFile.createVariable('YA', 'float64', ('time', 'depth', 'latitude', 'longitude'))
            ya.long_name = 'Acceleration along Y-axis'
            ya.coordinates = 'time depth latitude longitude'
            ya.units = 'g'
            ya[:] = self.ay.reshape(len(self.ay), 1, 1, 1)

            za = self.ncFile.createVariable('ZA', 'float64', ('time', 'depth', 'latitude', 'longitude'))
            za.long_name = 'Acceleration along X-axis'
            za.coordinates = 'time depth latitude longitude'
            za.units = 'g'
            za[:] = self.az.reshape(len(self.az), 1, 1, 1)

            a = self.ncFile.createVariable('A', 'float64', ('time', 'depth', 'latitude', 'longitude'))
            a.long_name = 'Acceleration Magnitude'
            a.coordinates = 'time depth latitude longitude'
            a.units = 'g'
            a[:] = self.a.reshape(len(self.a), 1, 1, 1)

            # Record Variables - Rotations
            xr = self.ncFile.createVariable('XR', 'float64', ('time', 'depth', 'latitude', 'longitude'))
            xr.long_name = 'Rotation about X-axis'
            xr.coordinates = 'time depth latitude longitude'
            xr.units = 'degree'
            xr[:] = (self.rx * 180 / np.pi).reshape(len(self.rx), 1, 1, 1)

            yr = self.ncFile.createVariable('YR', 'float64', ('time', 'depth', 'latitude', 'longitude'))
            yr.long_name = 'Rotation about Y-axis'
            yr.coordinates = 'time depth latitude longitude'
            yr.units = 'degree'
            yr[:] = (self.ry * 180 / np.pi).reshape(len(self.ry), 1, 1, 1)
    
            zr = self.ncFile.createVariable('ZR', 'float64', ('time', 'depth', 'latitude', 'longitude'))
            zr.long_name = 'Rotation about Z-axis'
            zr.coordinates = 'time depth latitude longitude'
            zr.units = 'degree'
            zr[:] = (self.rz * 180 / np.pi).reshape(len(self.rz), 1, 1, 1)

            mx = self.ncFile.createVariable('MX', 'float64', ('time', 'depth', 'latitude', 'longitude'))
            mx.long_name = 'X-component of rotation vector'
            mx.coordinates = 'time depth latitude longitude'
            mx.units = ''
            mx[:] = self.mx.reshape(len(self.mx), 1, 1, 1)

            my = self.ncFile.createVariable('MY', 'float64', ('time', 'depth', 'latitude', 'longitude'))
            my.long_name = 'Y-component of rotation vector'
            my.coordinates = 'time depth latitude longitude'
            my.units = ''
            my[:] = self.my.reshape(len(self.my), 1, 1, 1)

            mz = self.ncFile.createVariable('MZ', 'float64', ('time', 'depth', 'latitude', 'longitude'))
            mz.long_name = 'Z-component of rotation vector'
            mz.coordinates = 'time depth latitude longitude'
            mz.units = ''
            mz[:] = self.mz.reshape(len(self.mz), 1, 1, 1)

            rotrate = self.ncFile.createVariable('ROTRATE', 'float64', ('time', 'depth', 'latitude', 'longitude'))
            rotrate.long_name = 'Absolute rotation rate about rotation vector'
            rotrate.coordinates = 'time depth latitude longitude'
            rotrate.units = 'degree/second'
            rotrate[:] = self.rotrate.reshape(len(self.rotrate), 1, 1, 1)

            rotcount = self.ncFile.createVariable('ROTCOUNT', 'float64', ('time', 'depth', 'latitude', 'longitude'))
            rotcount.long_name = 'Rotation Count - Cumulative Sum of ROTRATE * dt / 360 deg'
            rotcount.coordinates = 'time depth latitude longitude'
            rotcount[:] = (self.rotcount).reshape(len(self.rotcount), 1, 1, 1)

            # Pressure sensor data interpolated to IMU samples
            p = self.ncFile.createVariable('P', 'float64', ('time','depth', 'latitude', 'longitude'))
            p.long_name = 'Pressure'
            p.coordinates = 'time depth latitude longitude'
            p.units = 'dbar'
            pres = np.interp(self.s2013, self.ps2013, self.pr * 10.0)
            p[:] = pres.reshape(len(pres), 1, 1, 1)


        elif self.featureType == 'trajectory':
            print "Writing trajectory data"
            meanLat = np.mean(self.traj_lat)
            print "Using meanLat =", meanLat, "for depth calculation"
            # Coordinate variables for trajectory 
            # Interpolate trajectory lat and lon onto the times of the data from the event - two time axes time and ptime mapped to lat,lon,depth
            self.latitude = self.ncFile.createVariable('latitude', 'float64', ('time',))
            self.latitude.long_name = 'LATITUDE'
            self.latitude.standard_name = 'latitude'
            self.latitude.units = 'degree_north'
            self.latitude[:] = np.interp(np.linspace(0,1,len(self.s2013)), np.linspace(0,1,len(self.traj_lat)), self.traj_lat)
    
            self.longitude = self.ncFile.createVariable('longitude', 'float64', ('time',))
            self.longitude.long_name = 'LONGITUDE'
            self.longitude.standard_name = 'longitude'
            self.longitude.units = 'degree_east'
            self.longitude[:] = np.interp(np.linspace(0,1,len(self.s2013)), np.linspace(0,1,len(self.traj_lon)), self.traj_lon)
    
            self.depth = self.ncFile.createVariable('depth', 'float64', ('time',))
            self.depth.long_name = 'DEPTH'
            self.depth.standard_name = 'depth'
            self.depth.units = 'm'
            self.depth[:] = np.interp(self.s2013, self.ps2013, eos80.dpth(self.pr * 10.0, meanLat))    # Convert pressure to depth

            # Record Variables - Accelerations
            xa = self.ncFile.createVariable('XA', 'float64', ('time',))
            xa.long_name = 'Acceleration along X-axis'
            xa.coordinates = 'time depth latitude longitude'
            xa.units = 'g'
            xa[:] = self.ax

            ya = self.ncFile.createVariable('YA', 'float64', ('time',))
            ya.long_name = 'Acceleration along Y-axis'
            ya.coordinates = 'time depth latitude longitude'
            ya.units = 'g'
            ya[:] = self.ay

            za = self.ncFile.createVariable('ZA', 'float64', ('time',))
            za.long_name = 'Acceleration along X-axis'
            za.coordinates = 'time depth latitude longitude'
            za.units = 'g'
            za[:] = self.az

            a = self.ncFile.createVariable('A', 'float64', ('time',))
            a.long_name = 'Acceleration Magnitude'
            a.coordinates = 'time depth latitude longitude'
            a.units = 'g'
            a[:] = self.a

            # Record Variables - Rotations
            xr = self.ncFile.createVariable('XR', 'float64', ('time',))
            xr.long_name = 'Rotation about X-axis'
            xr.standard_name = 'platform_roll_angle'
            xr.coordinates = 'time depth latitude longitude'
            xr.units = 'degree'
            xr[:] = (self.rx * 180 / np.pi)

            yr = self.ncFile.createVariable('YR', 'float64', ('time',))
            yr.long_name = 'Rotation about Y-axis'
            yr.standard_name = 'platform_pitch_angle'
            yr.coordinates = 'time depth latitude longitude'
            yr.units = 'degree'
            yr[:] = (self.ry * 180 / np.pi)
    
            zr = self.ncFile.createVariable('ZR', 'float64', ('time',))
            zr.long_name = 'Rotation about Z-axis'
            zr.standard_name = 'platform_yaw_angle'
            zr.coordinates = 'time depth latitude longitude'
            zr.units = ''
            zr[:] = (self.rz * 180 / np.pi)

            # Rotation about rotation vector - X3D format
            mx = self.ncFile.createVariable('MX', 'float64', ('time',))
            mx.long_name = 'X-component of rotation vector'
            mx.coordinates = 'time depth latitude longitude'
            mx.units = ''
            mx[:] = self.mx

            my = self.ncFile.createVariable('MY', 'float64', ('time',))
            my.long_name = 'Y-component of rotation vector'
            my.coordinates = 'time depth latitude longitude'
            my.units = ''
            my[:] = self.my

            mz = self.ncFile.createVariable('MZ', 'float64', ('time',))
            mz.long_name = 'Z-component of rotation vector'
            mz.coordinates = 'time depth latitude longitude'
            mz.units = ''
            mz[:] = self.mz

            rotrate = self.ncFile.createVariable('ROTRATE', 'float64', ('time',))
            rotrate.long_name = 'Absolute rotation rate about rotation vector'
            rotrate.coordinates = 'time depth latitude longitude'
            rotrate.units = 'degree/second'
            rotrate[:] = self.rotrate

            rotcount = self.ncFile.createVariable('ROTCOUNT', 'float64', ('time', ))
            rotcount.long_name = 'Rotation Count - Cumulative Sum of ROTRATE * dt / 360 deg'
            rotcount.coordinates = 'time depth latitude longitude'
            rotcount[:] = (self.rotcount)

            # Pressure sensor data interpolated to IMU samples
            p = self.ncFile.createVariable('P', 'float64', ('time',))
            p.long_name = 'Pressure'
            p.coordinates = 'time depth latitude longitude'
            p.units = 'dbar'
            p[:] = np.interp(self.s2013, self.ps2013, self.pr * 10.0)

        self.add_global_metadata()

        self.ncFile.close()

    def process_command_line(self):

        import argparse
        from argparse import RawTextHelpFormatter

        examples = 'Examples:' + '\n\n'
        examples += '  For 12 April 2013 BED01 deployment:\n'
        examples += '    ' + sys.argv[0] + " --input BED00048.EVT --output BED00048.nc --lat 36.793458 --lon -121.845703 --depth 295\n"
        examples += '  For 1 June 2013 BED01 Canyon Event:\n'
        examples += '    ' + sys.argv[0] + " --input BED00038.EVT --lat 36.793458 --lon -121.845703 --depth 340\n"
        examples += '    ' + sys.argv[0] + " --input BED00039.EVT --lat 36.785428 --lon -121.903602 --depth 530\n"
        examples += '    ' + sys.argv[0] + " --input BED00038.EVT BED00039.EVT --output BED01_1_June_2013.nc --trajectory BEDSLocation1_ThalwegTrace.csv \n"
        examples += '  For 18 February 2014 BED03 Canyon Event:\n'
        examples += '    ' + sys.argv[0] + " --input 30100046_partial_decimated10.EVT --lat 36.793367 --lon -121.8456035 --depth 292\n"

        parser = argparse.ArgumentParser(formatter_class=RawTextHelpFormatter,
                                         description='Convert BED event file(s) to a NetCDF file',
                                         epilog=examples)

        parser.add_argument('-i', '--input', action='store', nargs='*', required=True, help="Specify input event file name(s)")
        parser.add_argument('-o', '--output', action='store', help="Specify output NetCDF file name, if different from <base>.nc of input file name")
        parser.add_argument('-t', '--trajectory', action='store', help="csv file with columns of latitude and longitude where first and lat row corresponds to first and last records of input event files")
        parser.add_argument('--lat', type=float, action='store', help="latitude of BED device")
        parser.add_argument('--lon', type=float, action='store', help="longitude of BED device")
        parser.add_argument('--depth', type=float, action='store', help="depth of BED device")
        parser.add_argument('-v', '--verbose', type=int, choices=range(3), action='store', default=0, help="Specify verbosity level, values greater than 1 give more details ")

        self.args = parser.parse_args()

        if not self.args.input:
            parser.error("Must specify --input\n")

        if self.args.trajectory:
            pass
        else:
            if (self.args.lat and self.args.lon and self.args.depth):
                pass
            else:
                parser.error("If no --trajectory specified then must specify --lat, --lon, and --depth")

        self.commandline = ' '.join(sys.argv)
        self.inputFileNames = self.args.input
        if self.args.lat and self.args.lon and self.args.depth:
            self.lat = self.args.lat
            self.lon = self.args.lon
            self.dep = self.args.depth
            self.featureType = 'timeseries'
            if self.args.verbose > 0:
                print "self.lat = %f, self.lon = %f, self.dep = %f" % (self.lat, self.lon, self.dep)
        elif self.args.trajectory:
            self.featureType = 'trajectory'
            self.readTrajectory()
        else:
            raise Exception("Unknown featureType - must be timeseries or trajectory")

        for fileName in self.inputFileNames:
            if fileName.endswith('.EVT') or fileName.endswith('.WAT'):
                self.sensorType = 'Invensense'
            else:
                raise Exception("Unknown file: %s.  Input file must end with '.EVT' or '.WAT'." % fileName)
    
if __name__ == '__main__':

    beds_netcdf = BEDS_NetCDF()

    beds_netcdf.process_command_line()

    beds_netcdf.createNetCDFfromFile()

    print "Wrote file %s" % beds_netcdf.outFile