# Copyright Trimble Inc 2021

import numpy as np
import matplotlib.cm as cm
import matplotlib.pyplot as plt
import matplotlib.colors as col
from matplotlib.patches import Polygon
from PIL import Image
import seaborn as sns
import io
import os
import json
import requests
import sys
import mutils as m
import mutils.PosTypeConst as PosTypeConst
from ProcessResults import cdf_vals, get_rounded_cdf_percents

Headless = True
if(Headless == True):
  from matplotlib import use
  # Allow running headless from the command line
  use("agg")

# Override pylab find() to avoid deprecation warning
# pylint: disable=function-redefined
def find(x):
  return np.where(x)[0]

# Save a plot as a PNG file. Matplotlib isn't very good at 
# compressing PNG image files. Therefore, use PIL functions
# to compress a memory version of the file. The compressed 
# images are close to half the size of the Matplotlib generated
# version
def savePlot(directory,name):
  # Save the data as a PNG file
  if not os.path.exists(directory):
    os.makedirs(directory)

  # Save the png to memory
  ram = io.BytesIO()
  plt.savefig(ram,format='png',dpi=600)
  ram.seek(0)
  # Compress the png data before saving to the file system, the
  # matplotlib png's are not well compressed
  im = Image.open(ram)
  im2 = im.convert('RGB').convert('P', palette=Image.ADAPTIVE)
  im2.save(directory + "/" + name,format='PNG')

# Create a stacked bar chart showing the number of epochs as a function
# of distance to/from the obstruction. 10m bins are used.
def plotTotalEpochs(PVT,titleModifier,numOverpasses,path='.'):
  # Range to overpass in meters
  startBin = -195
  stopBin  = 205
  stepBin  = 10

  RTKData = []
  RTXData = []
  INSData = []
  StingerData = []
  AllData = []
  for dist in range(startBin, stopBin, stepBin):
    i = find( (PVT[:,2] > (dist-5.0)) & (PVT[:,2] <= (dist+5.0)) )
    filt = PVT[i,:]

    i = find( (filt[:,3] == PosTypeConst.dRec29_fixType_RTK_float) | (filt[:,3] == PosTypeConst.dRec29_fixType_RTK_fix) )
    RTKData.append(len(i))

    i = find(   (filt[:,3] == PosTypeConst.dRec29_fixType_RTX) 
              | ( (filt[:,3] >= PosTypeConst.dRec29_fixType_RTX_fast) & (filt[:,3] <= PosTypeConst.dRec29_fixType_RTX_fast_RAM) ) ) 
    RTXData.append(len(i))
    
    i = find( (filt[:,3] >= PosTypeConst.dRec29_fixType_INS_auto) & (filt[:,3] <= PosTypeConst.dRec29_fixType_INS_dead_reckoning) ) 
    INSData.append(len(i))

    ##postype = np.unique(filt[:,3])
    ##print("PosTypes Range = ",dist,"Total PVT",len(filt),"Types = ",postype)

    # Test - note Novatel report some code solutions as dead reckoning!
    #          | (filt[:,3] == PosTypeConst.dRec29_fixType_INS_dead_reckoning)  #### Temp - Novatel reports this
    
    i = find(   (filt[:,3] <= PosTypeConst.dRec29_fixType_SBAS)
              | (filt[:,3] >= PosTypeConst.dRec29_fixType_QZSS_SLAS)
              | (filt[:,3] == PosTypeConst.dRec29_fixType_GVBS)
              | ( (filt[:,3] >= PosTypeConst.dRec29_fixType_KF_auto) & (filt[:,3] <= PosTypeConst.dRec29_fixType_KF_SBASplus) ) ) 
    StingerData.append(len(i))
    AllData.append(len(filt))

  x_ = np.arange(startBin,stopBin,stepBin)
  
  # Create the figure
  fig, ax = plt.subplots()

  if(np.sum(RTKData) > 0):
    plt.bar(x_,RTKData,8,label='RTK')
  bottom = RTKData

  if(np.sum(RTXData) > 0):
    plt.bar(x_,RTXData,8,bottom=bottom,label='RTX')
    bottom = np.add(bottom,RTXData)
  
  if(np.sum(INSData) > 0):
    plt.bar(x_,INSData,8,bottom=bottom,label='INS')
    bottom = np.add(bottom,INSData)
  
  if(np.sum(StingerData) > 0):
    plt.bar(x_,StingerData,8,bottom=bottom,label='Code')
    bottom = np.add(bottom,StingerData)

  # See if we missed anything
  delta = np.subtract(AllData,bottom)
  if(np.sum(delta) > 0):
    plt.bar(x_,delta,8,bottom=bottom,label='Other')

  plt.xlim([-200,200]) # Range to/from the overpass
  plt.xlabel('Distance to Overpass (' + str(stepBin) + 'm bins) [m]')
  plt.ylabel('Number of Epochs')
  plt.legend()

  plt.title('Number of Epochs - ' + titleModifier + " (" + str(numOverpasses) + " Obstructions)")

  plt.grid(True)
  plt.tight_layout()
  savePlot(path,titleModifier + '-NumEpochs.png')

  if(not Headless):
    plt.show()

  if(plt):
    plt.close()


# Create a plot that shows the percentage of positions <= a set of thresholds. Data
# for each threshold is a different color; the data is plotted as a set of polygons.
# The X-axis is the distance to/from the obstruction in 10m steps
def plotErrvsRange(PVT,titleModifier,numOverpasses,RTKEngOnly=False,path='.'):
  if(RTKEngOnly == True):
    print("RTK Only")
  else:
    print("All Data")

  # PVT threshold in cm
  Max = 50
  StepSize = -5

  # Range to overpass in meters
  startBin = -195
  stopBin  = 205
  stepBin  = 10

  # We'll use this to get the RGB value to plot each polygon
  #cmap = matplotlib.cm.get_cmap('rainbow')
  cmap = cm.get_cmap('rainbow')

  fig, ax = plt.subplots()

  twoD = np.sqrt( PVT[:,8]**2 + PVT[:,9]**2 )
  for cmThres in range(Max,0,StepSize):
    thres = cmThres / 100.0 # Convert to meters
    percent = []
    for dist in range(startBin, stopBin, stepBin):
      i = find( (PVT[:,2] > (dist-5.0)) & (PVT[:,2] <= (dist+5.0)))
      if(RTKEngOnly == True):
        j = find( ( (PVT[:,3] == 3) | (PVT[:,3] == 4) ) & (PVT[:,2] > (dist-5.0)) & (PVT[:,2] <= (dist+5.0)) & (twoD <= thres) )
      else:
        j = find( (PVT[:,2] > (dist-5.0)) & (PVT[:,2] <= (dist+5.0)) & (twoD <= thres) )
      percent.append(100*len(j) / len(i))
   
    if(len(percent) > 0):
      x_ = np.arange(startBin,stopBin,stepBin)
      y_ = np.array(percent)

      # Now close the polygon
      x_ = np.append(x_,[x_[-1], x_[0], x_[0]])
      y_ = np.append(y_,[     0,     0, y_[0]])

      # Add the closed polygon to the plot setting the color appropriately
      # (based on threshold). There are a number of ways to show a colored
      # contour plot (e.g. contourf). However, most (all?) require a regular
      # x/y grid. We could manipulate the data, but as the data will change
      # monotonically (each polygon will overlap or be inside the last one)
      # we can simply plot a set of colored polygons.
      data = np.transpose([x_,y_])
      patches = Polygon(data,closed=True,color=cmap( (cmThres + StepSize)/(Max+StepSize)))
      ax.add_patch(patches)
      print("Percent for error <= ",thres,'m |',startBin,'m thru',stopBin-stepBin,'m in ',stepBin,'m steps')
      for val in percent:
        print("%.1f " % val,end="")
      print("") # Force a line return

  ax.set_xlim([-200,200]) # Range to/from the overpass
  ax.set_ylim([0,100])    # 0-100%

  bounds  = np.arange(0,Max-StepSize,-1*StepSize)
  # Location of the tick
  tickLoc = bounds[:-1] - (StepSize/2)
  # Value of the tick - we offset the value as we want
  # to output discrete values that match the plot at the 
  # midpoint of the color segment in the colorbar
  tickVal = bounds[:-1] - StepSize

  norm = col.BoundaryNorm(bounds, cmap.N)
  cb = fig.colorbar(cm.ScalarMappable(cmap='rainbow',norm=norm), ax=ax, ticks=tickLoc)

  numLabels = len(cb.ax.get_yticklabels())
  # Adjust the labels - we want to label the midpoint of each color segment
  for i in range(numLabels):
    handle = cb.ax.get_yticklabels()[i]
    val = tickVal[i]
    handle.set_text(str.format("{0:d}",int(val)))

  cb.ax.set_yticklabels(cb.ax.get_yticklabels())
  cb.set_label('PVT Threshold [cm]')

  plt.xlabel('Distance to Overpass (' + str(stepBin) + 'm bins) [m]')
  plt.ylabel('Percent of positions <= Threshold [%]')

  titleStr = '2D position error - ' + titleModifier
  if(RTKEngOnly == True):
    titleStr += ' (Fixed or Float only)'
  titleStr += " (" + str(numOverpasses) + " Obstructions)"
  plt.title(titleStr)

  plt.grid(True)
  plt.tight_layout()
  if(RTKEngOnly == True):
    savePlot(path,titleModifier + '-ErrvsRange-FixedFloat.png')
  else:
    savePlot(path,titleModifier + '-ErrvsRange.png')

  if(not Headless):
    plt.show()

  if(plt):
    plt.close()

# Create a plot that shows the error relative to the error estimates. The
# percentage of actual errors less than N x Estimated error is plotted as
# a set of polygons. The X-axis is the distance to/from the obstruction.
def plotErrEstvsRange(PVT,titleModifier,numOverpasses,RTKEngOnly=False,path='.'):

  if(RTKEngOnly == True):
    print("RTK Only")
    i = find( (PVT[:,3] == PosTypeConst.dRec29_fixType_RTK_float) | (PVT[:,3] == PosTypeConst.dRec29_fixType_RTK_fix) )
    PVT = PVT[i,:]
  else:
    print("All Data")

  # Max scalefactor * 100
  Max = 600
  StepSize = -50

  # Range to overpass in meters
  startBin = -195
  stopBin  = 205
  stepBin  = 10

  # We'll use this to get the RGB value to plot each polygon
  #cmap = matplotlib.cm.get_cmap('rainbow')
  cmap = cm.get_cmap('rainbow')

  fig, ax = plt.subplots()

  twoD = np.sqrt( PVT[:,8]**2 + PVT[:,9]**2 )
  twoDEst = np.sqrt( PVT[:,11]**2 + PVT[:,12]**2 )

  for sfInt in range(Max,0,StepSize):

    # sfInt is the scale factor multipled by 100 (need ints for range())
    sf = sfInt / 100.0
    percent = []
    for dist in range(startBin, stopBin, stepBin):
      testMetric = twoD < (twoDEst * sf)
      # All data in the bin
      i = find( (PVT[:,2] > (dist-5.0)) & (PVT[:,2] <= (dist+5.0)))
      # Just the data that meets the test metric
      j = find( (PVT[:,2] > (dist-5.0)) & (PVT[:,2] <= (dist+5.0)) & testMetric)
      if(len(i) > 0):
        percent.append(100*len(j) / len(i))
      else:
        percent.append(0)
   
    if(len(percent) > 0):
      x_ = np.arange(startBin,stopBin,stepBin)
      y_ = np.array(percent)

      # Now close the polygon
      x_ = np.append(x_,[x_[-1], x_[0], x_[0]])
      y_ = np.append(y_,[     0,     0, y_[0]])

      # Add the closed polygon to the plot setting the color appropriately
      # (based on threshold). There are a number of ways to show a colored
      # control plot (e.g. contourf). However, most (all?) require a regular
      # x/y grid. We could manipulate the data, but as the data will change
      # monotonically (each polygon will overlap or be inside the last one)
      # we can simply plot a set of colored polygons.
      data = np.transpose([x_,y_])
      #patches = Polygon(data,closed=True,color=cmap( (sf - 1)* 100.0/(Max-100)))
      patches = Polygon(data,closed=True,color=cmap( (sf - 0.5)* 100.0/(Max-50)))
      ax.add_patch(patches)

      print(percent)
      print(sf)

  ax.set_xlim([-200,200]) # Range to/from the overpass
  ax.set_ylim([0,100])    # 0-100%

  # Now add a colorbar and adjust the labels
  bounds  = np.arange(0,Max,-1*StepSize)
  # Location of the tick
  tickLoc = bounds[0::2] - (StepSize/2)
  # Value of the tick - we offset the value as we want
  # to output discrete values that match the plot at the 
  # midpoint of the color segment in the colorbar
  tickVal = (bounds[0::2] - 2*StepSize) / 100.0

  bounds  = np.arange(0,Max-StepSize,-1*StepSize)
  tickLoc = bounds[:-1] - (StepSize/2)
  tickVal = (bounds[:-1] - StepSize) / 100.0

  norm = col.BoundaryNorm(bounds, cmap.N)
  cb = fig.colorbar(cm.ScalarMappable(cmap='rainbow',norm=norm), ax=ax, ticks=tickLoc)
  numLabels = len(cb.ax.get_yticklabels())
  # Adjust the labels - we want to label the midpoint of each color segment
  for i in range(numLabels):
    handle = cb.ax.get_yticklabels()[i]
    val = tickVal[i]
    handle.set_text(str.format("{0:.2f}",val))

  cb.ax.set_yticklabels(cb.ax.get_yticklabels())
  cb.set_label('Estimated Error Scale Factor')

  plt.xlabel('Distance to Overpass (' + str(stepBin) + 'm bins) [m]')
  plt.ylabel('Percent Actual Error < Estimated Error * Scale Factor')
  
  titleStr = '2D error vs. est. error - ' + titleModifier
  if(RTKEngOnly == True):
    titleStr += ' (RTK)'
  titleStr += " (" + str(numOverpasses) + " Obstructions)"
  plt.title(titleStr)

  plt.grid(True)
  plt.tight_layout()
  if(RTKEngOnly == True):
    savePlot(path,titleModifier + '-FixedFloat-ErrEstvsErr.png')
  else:
    savePlot(path,titleModifier + '-ErrEstvsErr.png')

  if(not Headless):
    plt.show()

  if(plt):
    plt.close()

# Plot the estimated error from the "truth" (likely POSPac) solution
def plotTruth(PVT,titleModifier,numOverpasses,path='.'):
  fig, ax = plt.subplots()
  
  truth2D = np.sqrt(PVT[:,17]**2 + PVT[:,18]**2)
  plt.plot(PVT[:,2],truth2D,'.')
  plt.xlim([-200,200]) # Range to/from the overpass
  plt.xlabel('Distance to Overpass [m]')
  plt.ylabel('Truth (POSPac) reported 2D 1-sigma estimate [m]')
  plt.title('Estimated truth accuracy - ' + titleModifier + " (" + str(numOverpasses) + " Obstructions)")
  plt.grid(True)
  plt.tight_layout()
  savePlot(path,titleModifier + '-TruthErr.png')

  if(not Headless):
    plt.show()

  if(plt):
    plt.close()

# Create a Stanford plot that shows the error against the reported estimated error. All data is 
# processed in a single plot (irrespective of the distance to/from the overpass). The data is 
# partioned into 10cm bins. Each bin is colored based on the number of epochs (log10 scale) that
# fall into that cell.
def plotStanford(PVT,titleModifier,numOverpasses,RTKEngOnly=False,path='.'):
  ##################
  # Now form a stanford diagram showing the error against estimated error
  ##################

  if(RTKEngOnly == True):
    i = find( (PVT[:,3] == 3) | (PVT[:,3] == 4) )
    PVT = PVT[i,:]

  # 2D error estimate and 2D error in cm
  twoDEst = 100.0 * np.sqrt( PVT[:,11]**2 + PVT[:,12]**2 )
  twoD = 100.0 * np.sqrt( PVT[:,8]**2 + PVT[:,9]**2 )

  errStep = 10
  errMin = 0
  errMax = 200
  NumErr = int((errMax - errMin) / errStep )

  errData = np.zeros( (NumErr,NumErr) )
  mask    = np.zeros_like(errData)
  errIndex = 0
  for x in range(errMin,errMax,errStep):
    estIndex = 0
    for y in range(errMin,errMax,errStep):
      # Find data in the current bin, make sure that any data
      # off the plot is captured in the last cell on the edge of
      # the plot
      if( (x == (errMax - errStep)) and (y == (errMax - errStep)) ):
        i = find( (twoDEst >= y) & (twoD >= x))
      elif(x == (errMax - errStep)):
        i = find( (twoDEst >= y) & (twoDEst < (y+errStep)) & (twoD >= x) )
      elif(y == (errMax - errStep)):
        i = find( (twoDEst >= y) & (twoD >= x) & (twoD < (x+errStep)))
      else:
        i = find( (twoDEst >= y) & (twoDEst < (y+errStep)) & (twoD >= x) & (twoD < (x+errStep)))

      # the heatmap expects a matrix that is ordered [row,column]
      if(len(i) > 0):
        # errData 
        #   row is an index to the actual error (yaxis)
        #   col is an index to the estimated error (xaxis)
        errData[errIndex,estIndex] = np.log10(len(i))
      else:
        errData[errIndex,estIndex] = 0
        mask[errIndex,estIndex] = True


      estIndex += 1
    errIndex += 1

  print(errData)

  # Set the figure size, otherwise the x and y axis becomes crowded when saved to a png
  fig, ax = plt.subplots(figsize=(8,6.4))
  sns.heatmap(errData,ax=ax,xticklabels=2,yticklabels=2,mask=mask,vmin=0,vmax=4,cmap="rainbow")
  # the upper left is 0,0 in heatmap, invert the y-axis so the lower left is 0,0
  ax.invert_yaxis()

  # Adjust the axis labels from m into cm
  for handle in ax.get_xticklabels():
    val = int(handle.get_text())
    handle.set_text( str.format("{0:.3f}",0.5*errStep/100.0 + val * errStep/100.0) )
    ax.set_xticklabels(ax.get_xticklabels(),rotation=-45)

  for handle in ax.get_yticklabels():
    val = int(handle.get_text())
    handle.set_text( str.format("{0:.3f}",0.5*errStep/100.0 + val * errStep/100.0) )
    ax.set_yticklabels(ax.get_yticklabels(), rotation=0)

  # Adjust the color bar (convert from log to linear)
  cbar = ax.collections[0].colorbar
  cbar.set_ticks([0,1,2,3,4])
  cbar.set_ticklabels(['1','10','100','1,000','>= 10,000'])
  cbar.set_label('Number of epochs')


  plt.xlabel('Estimated 2D Error [m]')
  plt.ylabel('Actual 2D Error [m]')
  plt.grid(True)
  print(NumErr)
  plt.plot([0,NumErr],[0, NumErr],'--',color='blue', label='Actual Err. = Err. Estimate')
  plt.plot([0,NumErr],[0, 2*NumErr],'--',color='orange', label='Actual Err. = 2 x Err. Estimate')
  plt.plot([0,NumErr],[0, 3*NumErr],'--',color='red', label='Actual Err. = 3 x Err. Estimate')

  titleStr = '2D Estimated vs Actual position error - ' + titleModifier
  if(RTKEngOnly == True):
    titleStr += ' (Fixed or Float only)'
  titleStr += " (" + str(numOverpasses) + " Obstructions)"
  plt.title(titleStr)

  plt.legend()
  plt.tight_layout()
  if(RTKEngOnly == True):
    savePlot(path,titleModifier + '-vsErrFixedFloat-StanfordPlot.png')
  else:
    savePlot(path,titleModifier + '-vsErr-StanfordPlot.png')

  if(not Headless):
    plt.show()

  if(plt):
    plt.close()


# Gets overpass position info from SNPC directory.
# In: config = JSON config file which points to SNPC info
# Returns: (PVT,titleModifier)
#          PVT = overpass time-series info
#          titleModifier = text for plot titles
def getOverpassDataFromSNPCConfig( config ):
  import mutils
  from overpasses import find_overpasses
  from ProcessResults import gen_overpass_text_time_series
  pp = mutils.doload(config['SNPC_pospac'])
  d = mutils.vd2cls(config['SNPC_t04'],'-d35:2')
  titleModifier = config['SNPC_title']
  _, events = find_overpasses( pp, False )
  with open("temp_overpasses.txt","w") as f_out:
    gen_overpass_text_time_series( f_out, d, pp, events )
  PVT = np.loadtxt("temp_overpasses.txt",skiprows=1)
  os.remove("temp_overpasses.txt")
  return (PVT,titleModifier)

# Gets overpass position info from RF regression system.
# In: config = JSON config file which points to RF playback info
# Returns: (PVT,titleModifier)
#          PVT = overpass time-series info
#          titleModifier = text for plot titles
def getOverpassDataFromRFConfig( config ):
  RX = int(config['RX'])
  titleModifier = config['label']

  if(('Record16' in config) and (config['Record16'] == 'Y')):
    useRec16 = True
  else:
    useRec16 = False

  # See if we have direct access to the data or whether we need to use HTTP
  dataBase = '/net/fermion/mnt/data_drive/SpirentTest/GPSTools/RFPlaybackRegression/DataDir/'
  resBase  = '/net/fermion/mnt/data_drive/SpirentTest/GPSTools/RFPlaybackRegression/OutputResults/plots/'
  if( os.path.isdir(resBase) and os.path.isdir(dataBase)):
    fileRead = True
  else:
    fileRead = False # We'll get the data via http

  PVT = []
  if(fileRead == True): # the results directory is mounted so perform a direct read
    for run in config['runNum']: # The json file can include multiple runs - read and concat.
      if(useRec16):
        # "Under the hood" Stinger record 16 - if the system is operating in RTK mode,
        # the majority of these positions will be DGNSS.
        filename = resBase + 'RX' + str(RX) + '-' + str(run) + '/Stinger_overpass_time_series.txt'
      else:
        # Primary record - may contain RTK and Stinger solutions. If RTK is inactive this
        # may be identical to the record 16 data
        filename = resBase + 'RX' + str(RX) + '-' + str(run) + '/Main_overpass_time_series.txt'
      try:
        overpassData = np.loadtxt(filename,skiprows=1)
        numOverpasses = len(np.unique(overpassData[:,0]))
        if(len(PVT) == 0):
          PVT = overpassData
        else:
          passNum = np.max(PVT[:,0])
          overpassData[:,0] += passNum + 1
          PVT = np.concatenate((PVT,overpassData))
        print(np.shape(PVT))

        # The following demonstrates how we can get the configuration informtion
        # for this run. We may want to extend this to get other summary information.
        dataDir = dataBase + 'RX' + str(RX) + '-' + str(run)
        dirList = os.listdir(dataDir)
        for thisDir in dirList:
          if( (thisDir.endswith('.xml')) and (not thisDir.endswith('ALM.xml')) and (not thisDir.endswith('config.xml'))):
            fid = open(dataDir + '/' + thisDir,'r')
            soup = BeautifulSoup(fid, 'html.parser')
            fid.close()
            desc = soup.desc.text
            print(thisDir,"num overpasses = ",numOverpasses,desc)
            break

      except:
        print("Problem with run #",run)

  else: # Get the data via HTTP
    baseURL = 'http://fermion.eng.trimble.com:5000/RFPlaybackRegression/OutputResults/plots/RX'
    NumDataSets = len(config['runNum'])
    truth = []
    for run in config['runNum']:
      data = requests.get(baseURL + str(RX) + '-' + str(run) + '/Main_overpass_time_series.txt')
      try:
        if(data.status_code == 200):
          # Quick hack - find a way to read this from the buffer
          with open('tmp.txt','wb') as fid:
            fid.write(data.content)
            overpassData = np.loadtxt('tmp.txt',skiprows=1)
            numOverpasses = len(np.unique(overpassData[:,0]))
            if(len(PVT) == 0):
              PVT = overpassData
            else:
              passNum = np.max(PVT[:,0])
              overpassData[:,0] += passNum + 1
              PVT = np.concatenate((PVT,overpassData))
            print(np.shape(PVT))

          # Now find the config XML data file - this provides us with extra
          # description about the test. For HTTP access we get an HTML file with
          # the directory listing that we need to parse
          dataURL = 'http://fermion.eng.trimble.com:5000/RFPlaybackRegression/DataDir/RX'
          data = requests.get(dataURL + str(RX) + '-' + str(run) )
          soup = BeautifulSoup(data.text, 'html.parser')
          dirList = [node.get('href') for node in soup.find_all('a') if node.get('href').endswith('.xml')]
          for thisDir in dirList:
            if( (thisDir.endswith('.xml')) and (not thisDir.endswith('ALM.xml')) and (not thisDir.endswith('config.xml'))):
              data = requests.get(thisDir)
              soup = BeautifulSoup(data.text, 'html.parser')
              desc = soup.desc.text
              print(thisDir,"num overpasses = ",numOverpasses,desc)
              break
      except:
        print("Problem with run #",run)
  return (PVT,titleModifier)

# Load the data file, if a JSON file is provided multiple data sets are combined. Once loaded generate a set
# of graphs that are saved as PNG files.
if __name__ == "__main__":
  from bs4 import BeautifulSoup

  # data format:
  # pass_num time range_to_overpass fixtype lat lon hgt pdop dE dN dU sigmaE sigmaN sigmaU true_heading true_roll true_pitch

  if(len(sys.argv) == 1):
    print("Usage:")
    print(" python plotOverPass.py datafile.txt PlotTitle")
    print(" python plotOverPass.py config.json\n")
    print(" datafile.txt - overpass analysis file from the RF regression system")
    print(" PlotTitle    - string added to the titles and output (png) plot file names\n")
    print(" config.json  - configuration file - allows merging of multple input files")
    print("   For SNPC use, needs keys:")
    print("     SNPC_pospac : /path/to/POSPac_truth.txt")
    print("     SNPC_t04 : path/to/navAPP.t04")
    print("     SNPC_title : string to put on plot titles and filenames")
    print("   For RF regression system use, needs keys:")
    print("     RX : integer - receiver #")
    print("     titleModifier : string to put on plot titles and filenames")
    print("     Record16 : [optional] 'Y' to use rec16 data")
    print("     runNum : multiple(?) run numbers to analyze")
    sys.exit()
  elif(len(sys.argv) == 3):
    PVT = np.loadtxt(sys.argv[1],skiprows=1)
    titleModifier = sys.argv[2]
  elif( (len(sys.argv) == 2) and sys.argv[1].endswith('json')):
    # We have a json config file, this allows us to read in and combine
    # multiple results files
     
    with open(sys.argv[1],'r') as f:
      config = json.load(f)

    if 'SNPC_pospac' in config:
      PVT,titleModifier = getOverpassDataFromSNPCConfig( config )
    elif 'RX' in config:
      PVT,titleModifier = getOverpassDataFromRFConfig( config )
    else:
      raise RuntimeError("Unknown JSON config format")

  # Start by plotting the truth data.
  (row,cols) = np.shape(PVT)
  numOverpasses = len(np.unique(PVT[:,0]))
  print('Total Data Sets',numOverpasses)
  if(cols == 20): # Early versions of the RF Regression didn't include the POSPac truth
    plotTruth(PVT,titleModifier,numOverpasses)

  plotTotalEpochs(PVT,titleModifier,numOverpasses)

  plotErrvsRange(PVT,titleModifier,numOverpasses)
  plotErrvsRange(PVT,titleModifier,numOverpasses,RTKEngOnly=True)

  plotErrEstvsRange(PVT,titleModifier,numOverpasses)
  plotErrEstvsRange(PVT,titleModifier,numOverpasses,RTKEngOnly=True)

  plotStanford(PVT,titleModifier,numOverpasses)
  plotStanford(PVT,titleModifier,numOverpasses,RTKEngOnly=True)


  # Range to overpass in meters
  startBin = -195
  stopBin  = 205
  stepBin  = 10

  twoD = np.sqrt( PVT[:,8]**2 + PVT[:,9]**2 )
  cx2,cy2 = m.docdf( twoD )
  stats = get_rounded_cdf_percents( cx2, cy2 )
    
  # Find all the Stinger epochs - assume the rest of Titan/Astra
  stingerEpochs = find(   (PVT[:,3] <= PosTypeConst.dRec29_fixType_SBAS)
                        | (PVT[:,3] >= PosTypeConst.dRec29_fixType_QZSS_SLAS)
                        | (PVT[:,3] == PosTypeConst.dRec29_fixType_GVBS)
                        | (   (PVT[:,3] >= PosTypeConst.dRec29_fixType_KF_auto) 
                            & (PVT[:,3] <= PosTypeConst.dRec29_fixType_KF_SBASplus) ) ) 

  print('Epochs %d\nStinger %d\nOverpasses %d' % (len(twoD),len(stingerEpochs),numOverpasses))
  for k in range(len(stats)):
    print('%s%%= %.4f' % (cdf_vals[k],stats[k]))
        
  plt.figure()
  legendStr =  titleModifier + ' (' + str(len(twoD)) + ' Epochs)'
  plt.plot(cx2, cy2,label=legendStr)
  plt.legend(fontsize=8)
  ax=plt.gca()
  xRange = ax.get_xlim()
  plt.xlim(0,xRange[1])
  plt.ylim(0,1)
  plt.xlabel('Error [m]')
  plt.ylabel('2D CDF')
  plt.grid(True)
  plt.title('Data +/-200m of ' + str(numOverpasses) + ' Overpasses')
  plt.tight_layout()
  savePlot('.',titleModifier + '_2D_CDF.png')
  if(xRange[1] > 1.0):
    plt.xlim(0,1)
    plt.tight_layout()
    savePlot('.',titleModifier + '_2D_CDF-Max1m.png')
  
  if(xRange[1] > 0.2):
    plt.xlim(0,1)
    plt.tight_layout()
    savePlot('.',titleModifier + '_2D_CDF-Max20cm.png')


  cdfData = {}
  cdfData['title'] = titleModifier
  cdfData['numOverpasses'] = numOverpasses
  cdfData['numTotalEpochs'] = len(twoD)
  cdfData['numCodeEpochs'] = len(stingerEpochs)
  cdfData['x'] = list(cx2)
  cdfData['y'] = list(cy2)

  with open(titleModifier + '_cdf.json', 'w') as fid:
    json.dump(cdfData, fid)