import matplotlib
# Allow running headless from the command line
matplotlib.use("agg")

import math
import os
import datetime
import argparse
import sys;
import json;

from pylab import *
import numpy as np
import mutils as m
import pandas as pd
 

# Add stats to the figure legend
def plotLabel(stats,sigma, label):
  labelOut  = label + '\n'
  labelOut += r'$\sigma$ = ' + str.format('{0:.1f}', sigma * 1000.0) + ' mm/s'
  labelOut = plotPercentileLabel(stats,labelOut)
  return(labelOut)

def plotPercentileLabel(stats, label):
  labelOut  = label + '\n'
  labelOut += '68% = ' + str.format('{0:.1f}',stats['68%'] * 1000.0) + ' mm/s\n'
  labelOut += '95% = ' + str.format('{0:.1f}',stats['95%'] * 1000.0) + ' mm/s'
  return(labelOut)

def plotVel(figHandle,subplot,X,Y,percentiles,sigma,labelStr):
  ax=figHandle.add_subplot(subplot)
  plot(X,Y,label=plotLabel(percentiles,sigma,labelStr))
  ylabel('Velocity [m/s]')
  grid()
  legend(loc='best')

today = True

######################################################################
# Parse arguments
parser = argparse.ArgumentParser(description='Extracts and analyzes Velocity from a T04 file')
parser.add_argument('-f','--file', help='Parse a single file - provide the filename --file myFile.T04')
parser.add_argument('-d','--date', help='Optional start date year,month,day e.g. --date YYYY,MM,DD')
args = parser.parse_args()
######################################################################

if(args.date):
  if(args.file):
    print("Can't add a date if you provide a file for processing")
    sys.exit(1)
  date = args.date.split(",")
  today = False

if(args.file):
  MaxRX = 1
else:
  # Load the JSON file 
  with open('receivers.json', 'r') as f:
    receivers = json.load(f)

  MaxRX = len(receivers['receivers'])

for rx in range(MaxRX):
  if(today == True):
    #start  = datetime.datetime.now()
    # The NoPi files aren't available until the following day so get
    # yesterday's date
    start = datetime.datetime.now() - datetime.timedelta(days=1)
    num = 1
  else:
    start = datetime.datetime(int(date[0]),int(date[1]),int(date[2]),0,0,0)
    stop  = datetime.datetime.now()
    delta = ((stop-start).total_seconds()/86400)
    num = int(math.ceil(delta))

  if(args.file):
    checkExt = False
    if(os.name == 'nt'):
      # Running on windows
      RXStr = ''
    else:
      RXStr = '.'
  else:
    checkExt = True
    sys  = receivers['receivers'][rx]['sys']
    RXData = receivers['receivers'][rx]['RXDir']
    RXStr  = receivers['receivers'][rx]['RXStr']
    Ext  = receivers['receivers'][rx]['ext']

    # Some data sets have changed format overtime (e.g. from T02 to T04).
    # If this has occurred there will be an "oldExt" and an "oldDate" 
    # which define what the old file type was and when the transition occurred.
    # Check for this and change the extension time if we are before
    # the transition date.
    if('oldExt' in receivers['receivers'][rx]):
      print("got OldExt")
      OldExt  = receivers['receivers'][rx]['oldExt']
      OldDate = receivers['receivers'][rx]['oldDate']

    if(not os.path.exists(RXStr)):
      os.makedirs(RXStr)

  for index in range(num):
    day   = start.day
    month = start.month
    year  = start.year

    if(checkExt == True):
      # Some data sets have changed format overtime (e.g. from T02 to T04).
      # If this has occurred there will be an "oldExt" and an "oldDate" 
      # which define what the old file type was and when the transition occurred.
      # Check for this and change the extension time if we are before
      # the transition date.
      if('oldExt' in receivers['receivers'][rx]):
        YearOld  = int(OldDate[0:4])
        MonthOld = int(OldDate[4:6])
        DayOld   = int(OldDate[6:8])
        pivot = datetime.datetime(YearOld,MonthOld,DayOld)
        if((start-pivot).days < 0):
          Ext = OldExt

    start = start + datetime.timedelta(days=1)

    YYYYMMDD =  str(year) + str(month).zfill(2) + str(day).zfill(2) 
    if(args.file):
      filename = args.file
      resPrefix = filename.split('.')[-2]
    else:
      resPrefix = YYYYMMDD
      if(sys == "NoPi"):
        # Get the data from a NoPi directory
        filename  = "/net/daffy/mnt/data_drive/GNSSResults/" + str(year) + "/"
        filename +=  YYYYMMDD + "/" + RXData + "_" + YYYYMMDD + "." + Ext
      else:
        # Path is absolute
        filename = RXData + "/" + YYYYMMDD + "." + Ext

    print(filename)
    if(os.path.isfile(filename) and (os.path.getsize(filename) > 0)):
      (ref,k)=m.vd2arr(filename,rec='-d29',opt=['--dec=1000'])

      # The time gets extended by the week to handle week rollovers
      minWeek = ref[:,k.WEEK].min()
      # Load the data - convert time to hours
      data = pd.DataFrame({'East Vel':ref[:,k.VLON], 
                           'North Vel':ref[:,k.VLAT], 
                           'Up Vel':ref[:,k.VHGT],
                           'East abs(Vel)':abs(ref[:,k.VLON]), 
                           'North abs(Vel)':abs(ref[:,k.VLAT]), 
                           'Up abs(Vel)':abs(ref[:,k.VHGT]),
                           '2D Vel':np.sqrt(ref[:,k.VLON]**2 + ref[:,k.VLAT]**2)},
                           index=((ref[:,k.TIME] + ((ref[:,k.WEEK]-minWeek)*7*86400))/3600))
      # Define the time label - if we change the time to seconds,
      # changing this string will update the x-axis label
      data.index.name = 'Time [Hour in GPS Week]'
      # Get the stats including defined percentiles
      Stats = data.describe(percentiles=[0.68,0.95,0.99,0.999]) 

      # Plot the data - for the percentiles use the "abs()" version
      fig = figure()
      plotVel(fig,311,data.index,data['East Vel'],Stats['East abs(Vel)'],Stats['East Vel']['std'],'East Velocity')
      plotVel(fig,312,data.index,data['North Vel'],Stats['North abs(Vel)'],Stats['North Vel']['std'],'North Velocity')
      plotVel(fig,313,data.index,data['Up Vel'],Stats['Up abs(Vel)'],Stats['Up Vel']['std'],'Up Velocity')
      xlabel(data.index.name)
      tight_layout()
      show()
      # Save the data as a PNG file
      if(os.name == 'nt'):
        savefig(resPrefix + "-velocity.png",dpi=150)
      else:
        savefig(RXStr + "/" + resPrefix + "-velocity.png",dpi=150) 
      close()

      # Plot the 2D data
      fig = figure()
      plot( data.index, data['2D Vel'],  label=plotPercentileLabel(Stats['2D Vel'],'Horz Velocity'))
      ylabel('Velocity [m/s]')
      grid()
      legend(loc='best')
      xlabel(data.index.name)
      tight_layout()
      show()
      # Save the data as a PNG file
      if(os.name == 'nt'):
        savefig(resPrefix + "-velocity2D.png",dpi=150)
      else:
        savefig(RXStr + "/" + resPrefix + "-velocity2D.png",dpi=150)
      close()

      # Now plot a CDF of E/N/U
      fig = figure()
      dataLen = len(data.index)
      Y = np.arange(0,dataLen,1) / (dataLen-1)
      plot( np.sort(data['East abs(Vel)']),  Y, label=plotPercentileLabel(Stats['East abs(Vel)'],'East Velocity'))
      plot( np.sort(data['North abs(Vel)']), Y, label=plotPercentileLabel(Stats['North abs(Vel)'],'North Velocity'))
      plot( np.sort(data['Up abs(Vel)']),    Y, label=plotPercentileLabel(Stats['Up abs(Vel)'],'Up Velocity'))
      plot( np.sort(data['2D Vel']),         Y, label=plotPercentileLabel(Stats['2D Vel'],'2D Velocity'))
      xlabel('Velocity [m/s]')
      ylabel('CDF')
      xStart,xStop = xlim()
      xlim(0.0,xStop)
      ylim(0.0,1.0)
      grid()
      legend(loc='best')
      tight_layout()
      show()
      # Save the data as a PNG file
      if(os.name == 'nt'):
        savefig(resPrefix + "-velocityCDF.png",dpi=150)
      else:
        savefig(RXStr + "/" + resPrefix + "-velocityCDF.png",dpi=150)
      close()

      if(not args.file):
        fid = open(RXStr + "/VelocityStats.txt","a")
        fid.write(  "%d %d %d " % (year,month,day))
        fid.write(  "%.0f %.5f %.5f %.5f %.5f %.5f %.5f " %
                   (Stats['East Vel']['count'],
                    Stats['East Vel']['mean'],
                    Stats['East Vel']['std'],
                    Stats['East abs(Vel)']['68%'],
                    Stats['East abs(Vel)']['95%'],
                    Stats['East abs(Vel)']['99%'],
                    Stats['East abs(Vel)']['max']))
        fid.write(  "%.5f %.5f %.5f %.5f %.5f %.5f " %
                   (Stats['North Vel']['mean'],
                    Stats['North Vel']['std'],
                    Stats['North abs(Vel)']['68%'],
                    Stats['North abs(Vel)']['95%'],
                    Stats['North abs(Vel)']['99%'],
                    Stats['North abs(Vel)']['max']))
        fid.write(  "%.5f %.5f %.5f %.5f %.5f %.5f " %
                   (Stats['Up Vel']['mean'],
                    Stats['Up Vel']['std'],
                    Stats['Up abs(Vel)']['68%'],
                    Stats['Up abs(Vel)']['95%'],
                    Stats['Up abs(Vel)']['99%'],
                    Stats['Up abs(Vel)']['max']))
        fid.write(  "%.5f %.5f %.5f %.5f\n" %
                   (Stats['2D Vel']['68%'],
                    Stats['2D Vel']['95%'],
                    Stats['2D Vel']['99%'],
                    Stats['2D Vel']['max']))
        fid.close()
      
      # Output summary statistics to stdout
      print(Stats)