# Copyright Trimble Inc 2021 # # You'll need something like this to setup the python path: # export PYTHONPATH=$HOME/GPSTools/pythonTools/:$HOME/GPSTools/RFPlaybackRegression/:$HOME/GPSTools/pythonTools/obstructions/ # import numpy as np #import io import os import json import sys import plotOverpass as obs import datetime import mutils as m import mutils.PosTypeConst as PosTypeConst from ProcessResults import cdf_vals, get_rounded_cdf_percents import matplotlib.pyplot as plt 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] if __name__ == "__main__": # 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 now = datetime.datetime.now() YYYYMMDD = str(now.year) + '-' + str(now.month).zfill(2) + '-' + str(now.day).zfill(2) if not os.path.exists(YYYYMMDD): os.makedirs(YYYYMMDD) if not os.path.exists("Trend"): os.makedirs("Trend") if(len(sys.argv) == 2): with open(sys.argv[1],'r') as f: config = json.load(f) else: print("Usage:\n python processOverpass.py config.json") sys.exit() plt.figure() summaryFig = plt.gcf().number firstRX = None processedDataSet = [] for RXRec in config['RX']: RX = RXRec['ID'] longLabel = RXRec['label'] # Too long for the title in many graphs titleModifier = 'RX' + str(RX) print("RX",RX,titleModifier) if(('rec' in RXRec) and (RXRec['rec'] == '-d35:16')): useRec16 = True titleModifier += '_StingerRec' else: useRec16 = False if(firstRX is None): firstRX = RX # 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/' PVT = [] for run in config['dataSets']: # The json file can include multiple runs - read and concat. print("Find run",run) # The data directory has a symlink to the latest data set. By expanding the symlink we can determine # the run number, we use this in the path to the results. latest = dataBase + 'latest-RX' + str(RX) + '-' + run + '/' # Now break out the run ID from the full path (expand the symlink) fullPath= os.path.realpath(latest) RunID = (fullPath.split('/')[-1]).split('-')[1] # Save this so we can provide a summary at the end if((firstRX == RX) and (useRec16 == False)): processedDataSet.append(RunID + ' ' + run + '.xml') filename = resBase + 'RX' + str(RX) + '-' + str(RunID) 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 += '/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 += '/Main_overpass_time_series.txt' print(filename) try: overpassData = np.loadtxt(filename,skiprows=1) numOverpasses = len(np.unique(overpassData[:,0])) if(numOverpasses > 0): if(len(PVT) == 0): PVT = overpassData else: passNum = np.max(PVT[:,0]) overpassData[:,0] += passNum + 1 PVT = np.concatenate((PVT,overpassData)) print("Got Data",np.shape(PVT)) except: print("Problem with run #",run) # 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 obs.plotTruth(PVT,titleModifier,numOverpasses,path=YYYYMMDD) obs.plotTotalEpochs(PVT,titleModifier,numOverpasses,path=YYYYMMDD) obs.plotErrvsRange(PVT,titleModifier,numOverpasses,path=YYYYMMDD) if(useRec16 == False): # No way we can get RTK data if we are using record 35:16s obs.plotErrvsRange(PVT,titleModifier,numOverpasses,RTKEngOnly=True,path=YYYYMMDD) obs.plotErrEstvsRange(PVT,titleModifier,numOverpasses,path=YYYYMMDD) if(useRec16 == False): # No way we can get RTK data if we are using record 35:16s obs.plotErrEstvsRange(PVT,titleModifier,numOverpasses,RTKEngOnly=True,path=YYYYMMDD) obs.plotStanford(PVT,titleModifier,numOverpasses,path=YYYYMMDD) if(useRec16 == False): # No way we can get RTK data if we are using record 35:16s obs.plotStanford(PVT,titleModifier,numOverpasses,RTKEngOnly=True,path=YYYYMMDD) # Now form some statistics # 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 ) print(stats) # 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) ) ) # Append the stats so we can trend over time. trendStr = 'Trend/RX' + str(RX) + '-' if(useRec16 == True): trendStr += 'Rec35_16.csv' else: trendStr += 'Rec35_2.csv' with open(trendStr,'a') as fid: fid.write('%d,%d,%d,%d,%d,%d,' % (now.year,now.month,now.day,len(twoD),len(stingerEpochs),numOverpasses)) for k in range(len(stats)): fid.write('%.4f' % (stats[k])) if(k < (len(stats)-1)): fid.write(',') else: fid.write('\n') # end of the data line plt.figure(summaryFig) legendStr = 'RX' + str(RX) + ' ' + longLabel + ' (' + str(len(twoD)) + ' Epochs)' plt.plot(cx2, cy2,label=legendStr) plt.grid(True) # Warning - assuming the num overpasses is the same for each receiver. # This should be correct, unless one receiver stops processing and the # overpass map is updated plt.title('Data +/-200m of ' + str(numOverpasses) + ' Overpasses') 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('CDF') plt.tight_layout() obs.savePlot(YYYYMMDD,YYYYMMDD + '_CDF.png') if(xRange[1] > 1.0): plt.xlim(0,1) plt.tight_layout() obs.savePlot(YYYYMMDD,YYYYMMDD + '_CDF-Max1m.png') if(xRange[1] > 0.2): plt.xlim(0,1) plt.tight_layout() obs.savePlot(YYYYMMDD,YYYYMMDD + '_CDF-Max20cm.png') fid = open(YYYYMMDD + '/' + YYYYMMDD + '-datasets.txt','w') for dataSet in processedDataSet: fid.write("%s\n" % dataSet) fid.close()