import sys import os import time import signal import datetime import matplotlib # Allow running headless from the command line matplotlib.use("agg") from numpy import * from pylab import * RXData = [ ["BD935_1_ROVER", "BD935_1_ROVER", "-a0"], ["BD935_2_ROVER", "BD935_2_ROVER", "-a0"], ["BD935_3_ROVER", "BD935_3_ROVER", "-a0"], ["BD940_1_ROVER", "BD940_1_ROVER", "-a0"], ["BD940_2_ROVER", "BD940_2_ROVER", "-a0"], ["BD940_3_ROVER", "BD940_3_ROVER", "-a0"], ["BD992_ROVER", "BD992_ROVER_A0","-a0"], ["BD992_ROVER", "BD992_ROVER_A1","-a1"] ] MaxRX = len(RXData) def signal_handler(signal, frame): print('\nInterrupted with CTRL-C downloading stopped....\n\n') sys.exit(0) signal.signal(signal.SIGINT, signal_handler) ############################################################# # Main code below ############################################################# ElevModel = [] now = time.time() dict = {'MP12' : 0, 'MP21' : 1, 'MP15' : 2, 'MP51' : 3, 'MP16' : 4, 'MP61' : 5, 'MP17' : 6, 'MP71' : 7, 'MP18' : 8, 'MP81' : 9 } islip = {'i12' : 0, 'i21' : 1, 'i15' : 2, 'i51' : 3, 'i16' : 4, 'i61' : 5, 'i17' : 6, 'i71' : 7, 'i18' : 8, 'i81' : 9 } mslip = {'m12' : 0, 'm21' : 1, 'm15' : 2, 'm51' : 3, 'm16' : 4, 'm61' : 5, 'm17' : 6, 'm71' : 7, 'm18' : 8, 'm81' : 9 } LUT = {0 : 'MP12', 1 : 'MP21', 2 : 'MP15', 3 : 'MP51', 4 : 'MP16', 5 : 'MP61', 6 : 'MP17', 7 : 'MP71', 8 : 'MP18', 9 : 'MP81' } def plotMPxx(filename): Day = [] MPXX = [] Slips = [] with open(filename,'r') as f: for line in f: # Col0 = Year # Col1 = Month # Col2 = Day # col3 = MPxx (summary MPxx for all data of this type) # Col4 = Num mslip epochs < 10 deg # Col5 = Num islip epochs < 10 deg # Col6 = Num mslip epochs >= 10 deg && < 20 # Col7 = Num islip epochs >= 10 deg && < 20 # Col8 = Num mslip epochs >= 20 deg # Col9 = Num islip epochs >= 20 deg data = line.rstrip().split() if(float(data[3]) < 1.0): # Remove any corrupt data if(int(data[1]) == 1): Day.append(31 + int(data[2])) else: Day.append(int(data[2])) MPXX.append(float(data[3])) Slips.append(int(data[7]) + int(data[9])) # First process the MPXX fig=figure() ax=fig.add_subplot(111) plot( array(Day), array(MPXX)) xlabel(r'Time [Days since 2017-12-01]') path, file_tmp = os.path.split(filename) file_start, file_extension = os.path.splitext(file_tmp) ylabel(file_start + ' [m]') grid(True) tight_layout() # Prevent the axis numers having an offset ax.get_xaxis().get_major_formatter().set_useOffset(False) ax.get_yaxis().get_major_formatter().set_useOffset(False) show() # Save the data as a PNG file savefig(filename + ".png",dpi=150) close() # Now process the cycle slips fig=figure() ax=fig.add_subplot(111) plot( array(Day), array(Slips)) xlabel(r'Time [Days since 2017-12-01]') path, file_tmp = os.path.split(filename) file_start, file_extension = os.path.splitext(file_tmp) ylabel(file_start + '-Slips (>10 deg) [number/day]') grid(True) tight_layout() # Prevent the axis numers having an offset ax.get_xaxis().get_major_formatter().set_useOffset(False) ax.get_yaxis().get_major_formatter().set_useOffset(False) show() # Save the data as a PNG file savefig(filename + "-Slips.png",dpi=150) close() def getElevData(data,dataType): if((len(data) >= 6) and (data[0].isdigit()) and (int(data[0])/5 <= 17)): index = int(1.0 + int(data[0])/5) ElevModel[0][index] = int(data[3]) ElevModel[1][index] = int(data[4]) ElevModel[2][index] = float(data[5]) # Determine whether it is time to terminate the parsing if((len(data) >= 1) and (data[0] == "<")): # Complete with the data from below the horizon ElevModel[0][0] = int(data[2]) ElevModel[1][0] = int(data[3]) ElevModel[2][0] = float(data[4]) #print dataType, ElevModel return {'processData':False, 'ElevModel':ElevModel} else: return {'processData':True} def timeStr2Secs(string): data = string.split(':') time = 3600 * int(data[0]) time = time + 60 * int(data[1]) time = time + float(data[1]) return time def parse17Sfile(fileroot,resroot): islipNum = [] mslipNum = [] print(fileroot,resroot) for i in range(3): islipNum.append([0] * 10) mslipNum.append([0] * 10) filename = "%s.Slips" % (fileroot) mp12Slip = [] mp12Time = [] with open(filename,'r') as f: for line in f: data = line.rstrip().split() if(float(data[1]) < 10.0): if(islip.get(data[3]) is not None): islipNum[0][islip.get(data[3])] = islipNum[0][islip.get(data[3])] + 1 if(mslip.get(data[3]) is not None): mslipNum[0][mslip.get(data[3])] = mslipNum[0][mslip.get(data[3])] + 1 elif(float(data[1]) < 20.0): if(data[3] == 'i12'): mp12Slip.append(float(data[1])) mp12Time.append(timeStr2Secs(data[9])) if(islip.get(data[3]) is not None): islipNum[1][islip.get(data[3])] = islipNum[1][islip.get(data[3])] + 1 if(mslip.get(data[3]) is not None): mslipNum[1][mslip.get(data[3])] = mslipNum[1][mslip.get(data[3])] + 1 else: if(islip.get(data[3]) is not None): islipNum[2][islip.get(data[3])] = islipNum[2][islip.get(data[3])] + 1 if(mslip.get(data[3]) is not None): mslipNum[2][mslip.get(data[3])] = mslipNum[2][mslip.get(data[3])] + 1 # Now output a graph of MP12 cycle slips Elev Vs Time fig=figure() ax=fig.add_subplot(111) plot( array(mp12Time), array(mp12Slip),'x') xlabel(r'Time [Seconds into Day]') ylabel('Elev Angle [deg]') grid(True) tight_layout() # Prevent the axis numers having an offset ax.get_xaxis().get_major_formatter().set_useOffset(False) ax.get_yaxis().get_major_formatter().set_useOffset(False) show() # Save the data as a PNG file savefig(fileroot + ".Slips.png",dpi=150) close() filename = "%s.MPX" % (fileroot) foundMPxx = False MPXXType = None gotObsSummary = False for i in range(3): ElevModel.append([None] * 19) with open(filename,'r') as f: for line in f: data = line.rstrip().split() if(len(data) >= 6 and data[0] == "Moving" and data[1] == "average"): if( dict.get(data[2]) is not None ): MPXX[ dict.get(data[2]) ] = float(data[4]) elif(len(data) >= 12 and data[0] == "first" and data[4] == "hrs"): gotObsSummary = True elif((len(data) >= 17) and (gotObsSummary == True)): gotObsSummary = False if( (data[11].isdigit()) and (data[12].isdigit()) and (data[16].isdigit())): expected = int(data[11]) have = int(data[12]) obsPerSlip = int(data[16]) gotDigits = True else: gotDigits = False # This is a legacy summary from Teqc, it only provides info # for the MP12/MP21 combinations. Provides obs per slip # info that some customers use if( (gotDigits == True) and (obsPerSlip > 0) ): fid = open("%s-summary.txt" % (resroot),'a') fid.write("%4d %d %d " % (year, month, day) ) numSlips = expected/obsPerSlip fid.write("%d %d %d %d\n" % (expected, have, obsPerSlip,numSlips) ) fid.close() elif( (len(data) >= 3) and (data[0] == "mean") and (data[2] == "rms") and (dict.get(data[1]) is not None) ): # Now look for the data as a function of elevation angle MPxxType = dict.get(data[1]) foundMPxx = True elif(foundMPxx == True): ret =getElevData(data, MPxxType) foundMPxx = ret['processData'] if(foundMPxx == False): # We have the elevation dependent MPxx and epoch/slip # information so save it. # # File information is as follows: # # Col1 = Year # Col2 = Month # Col3 = Day # col4 = MPxx (summary MPxx for all data of this type) # Col5 = Num mslip epochs < 10 deg # Col6 = Num islip epochs < 10 deg # Col7 = Num mslip epochs >= 10 deg && < 20 # Col8 = Num islip epochs >= 10 deg && < 20 # Col9 = Num mslip epochs >= 20 deg # Col10 = Num islip epochs >= 20 deg # Cols 11 - 29 = <0,0-5,5-10,10-15..,85-90 - Total epochs # this elev range # Cols 30 - 48 = <0,0-5,5-10,10-15..,85-90 - Total slips # this elev range # Cols 49 - 67 = <0,0-5,5-10,10-15..,85-90 - MPxx this elev # range fid = open("%s-%s.txt" % (resroot, LUT.get(MPxxType)),'a') fid.write("%4d %d %d " % (year, month, day) ) fid.write("%.4f " % MPXX[MPxxType]) # Output the information we parsed from the slip file fid.write("%d %d %d %d %d %d " % (mslipNum[0][MPxxType], islipNum[0][MPxxType], mslipNum[1][MPxxType], islipNum[1][MPxxType], mslipNum[2][MPxxType], islipNum[2][MPxxType]) ) for i in range(19): fid.write(str(ret['ElevModel'][0][i]) + " ") for i in range(19): fid.write(str(ret['ElevModel'][1][i]) + " ") for i in range(19): fid.write(str(ret['ElevModel'][2][i]) + " ") fid.write("\n") fid.close() RXLabel = [ "BD935_1", "BD935_2", "BD935_3", "BD940_1", "BD940_2", "BD940_3", "BD992-Pos", "BD992-Vec"] plotSymb = [ "bx", "rx", "cx", "gx", "kx", "yx", "b*", "r*"] file2Str = {'-noGLN-L2E-MP12.txt' : 'MP12 GPS (L1 CA / L2E)', '-noGLN-L2C-MP12.txt' : 'MP12 GPS (L1 CA / L2C)', '-GLN_CA-MP12.txt' : 'MP12 GLN (L1 CA / L2 CA)', '-GLN_P-MP12.txt' : 'MP12 GLN (L1 P / L2 P)', '-noGLN-L2E-MP15.txt' : 'MP15 GPS (L1 CA / L5) | GAL (E1 / E5A)', '-noGLN-L2E-MP16.txt' : 'MP16 GAL (E1 / E6) | BDS (B1 / B3)', '-noGLN-L2E-MP17.txt' : 'MP17 GAL (E1 / E5B) | BDS (B1 / B2)', '-noGLN-L2E-MP18.txt' : 'MP18 GAL E1 / E5AltBOC', '-noGLN-L2E-MP21.txt' : 'MP21 GPS (L2E / L1 CA)', '-noGLN-L2C-MP21.txt' : 'MP21 GPS (L2C / L1 CA)', '-GLN_CA-MP21.txt' : 'MP21 GLN (L2 CA / L1 CA)', '-GLN_P-MP21.txt' : 'MP21 GLN (L2 P / L1 P)', '-noGLN-L2E-MP51.txt' : 'MP51 GPS (L5 / L1 CA) | GAL (E5A / E1)', '-noGLN-L2E-MP61.txt' : 'MP61 GAL (E6 / E1) | BDS (B3 / B2)', '-noGLN-L2E-MP71.txt' : 'MP71 GAL (E5B / E51) | BDS (B2 / B1)', '-noGLN-L2E-MP81.txt' : 'MP81 GAL (E5AltBOC / E1)' } file2name = {'-noGLN-L2E-MP12.txt' : 'MP12-GPS_L1CA_L2E', '-noGLN-L2C-MP12.txt' : 'MP12-GPS_L1CA_L2C', '-GLN_CA-MP12.txt' : 'MP12-GLN_L1CA_L2CA', '-GLN_P-MP12.txt' : 'MP12-GLN_L1P_L2P', '-noGLN-L2E-MP15.txt' : 'MP15-GPS_L1CA_L5__GAL_E1_E5A', '-noGLN-L2E-MP16.txt' : 'MP16-GAL_E1_E6__BDS_B1_B3', '-noGLN-L2E-MP17.txt' : 'MP17-GAL_E1_E5B__BDS_B1_B2', '-noGLN-L2E-MP18.txt' : 'MP18-GAL_E1_E5AltBOC', '-noGLN-L2E-MP21.txt' : 'MP21-GPS_L2E_L1CA', '-noGLN-L2C-MP21.txt' : 'MP21-GPS_L2C_L1CA', '-GLN_CA-MP21.txt' : 'MP21-GLN_L2CA_L1CA', '-GLN_P-MP21.txt' : 'MP21-GLN_L2P_L1P', '-noGLN-L2E-MP51.txt' : 'MP51-GPS_L5_L1CA__GAL_E5A_E1', '-noGLN-L2E-MP61.txt' : 'MP61-GAL_E6_E1__BDS_B3_B1', '-noGLN-L2E-MP71.txt' : 'MP71-GAL_E5B_E1__BDS_B2_B1', '-noGLN-L2E-MP81.txt' : 'MP81-GAL_E5AltBOC_E1' } file2type = {'-noGLN-L2E-MP12.txt' : 'MP12', '-noGLN-L2C-MP12.txt' : 'MP12', '-GLN_CA-MP12.txt' : 'MP12', '-GLN_P-MP12.txt' : 'MP12', '-noGLN-L2E-MP15.txt' : 'MP15', '-noGLN-L2E-MP16.txt' : 'MP16', '-noGLN-L2E-MP17.txt' : 'MP17', '-noGLN-L2E-MP18.txt' : 'MP18', '-noGLN-L2E-MP21.txt' : 'MP21', '-noGLN-L2C-MP21.txt' : 'MP21', '-GLN_CA-MP21.txt' : 'MP21', '-GLN_P-MP21.txt' : 'MP21', '-noGLN-L2E-MP51.txt' : 'MP51', '-noGLN-L2E-MP61.txt' : 'MP61', '-noGLN-L2E-MP71.txt' : 'MP71', '-noGLN-L2E-MP81.txt' : 'MP81' } def plotMPXTrend(fileSuffix): havePlot = False fig=figure() ax=fig.add_subplot(111) for i in range(MaxRX): DOY = [] MPX = [] gotData = False stop = datetime.datetime.now() endYear = stop.year for j in range(endYear - 2018 + 1): filename = str(2018 + j) + '/' + RXData[i][1] + fileSuffix print(filename) if(os.path.isfile(filename)): with open(filename,'r') as f: for line in f: # Col0 = Year # Col1 = Month # Col2 = Day # col3 = MPxx (summary MPxx for all data of this type) # Col4 = Num mslip epochs < 10 deg # Col5 = Num islip epochs < 10 deg # Col6 = Num mslip epochs >= 10 deg && < 20 # Col7 = Num islip epochs >= 10 deg && < 20 # Col8 = Num mslip epochs >= 20 deg # Col9 = Num islip epochs >= 20 deg data = line.rstrip().split() if(float(data[3]) < 1.0): # Remove any corrupt data DOY.append( ( datetime.datetime((int)(data[0]), (int)(data[1]), (int)(data[2]), 0, 0) - datetime.datetime(2018,1,1)).days) MPX.append(float(data[3])) gotData = True # Now plot the data if(gotData == True): plot( array(DOY),array(MPX),plotSymb[i],label=RXLabel[i]) havePlot = True # Now create a single station plot figSingle = figure() axSingle = figSingle.add_subplot(111) plot( array(DOY),array(MPX),plotSymb[i],label=RXLabel[i]) xlabel('Time [Days since 2018-01-01]') ylabel(file2type.get(fileSuffix) + ' [m]') title(file2Str.get(fileSuffix)) grid(True) tight_layout() # Prevent the axis numers having an offset axSingle.get_xaxis().get_major_formatter().set_useOffset(False) axSingle.get_yaxis().get_major_formatter().set_useOffset(False) show() # Save the data as a PNG file savefig(file2name.get(fileSuffix) + '-' + RXData[i][1] + '.png', dpi=150) close() # Now switch back to the combined plot figure(fig.number) if(havePlot == True): xlabel('Time [Days since 2018-01-01]') ylabel(file2type.get(fileSuffix) + ' [m]') title(file2Str.get(fileSuffix)) grid(True) legend() tight_layout() # Prevent the axis numers having an offset ax.get_xaxis().get_major_formatter().set_useOffset(False) ax.get_yaxis().get_major_formatter().set_useOffset(False) show() # Save the data as a PNG file savefig(file2name.get(fileSuffix) + '.png', dpi=150) close() ##### ##### ##### Code Start ##### ##### datetimeformat = "%Y-%m-%d %H:%M:%S" # Grab the date once in case we span into the next day start = datetime.datetime.now() #start = datetime.datetime.now() - datetime.timedelta(days=2) for RXNum in range(0,MaxRX): #start = datetime.datetime.strptime("2018-01-01 12:00:00",datetimeformat) #start = datetime.datetime.strptime("2018-02-13 12:00:00",datetimeformat) #stop = datetime.datetime.strptime("2018-02-14 13:00:00",datetimeformat) ###delta = ((stop-start).total_seconds()/86400) #for index in range((int)(ceil(delta))): if True: day = start.day month = start.month year = start.year # Bump the date for next time through ###start = start + datetime.timedelta(days=1) print(year,month,day,RXNum) if(True): # Set to False if you don't want to re-run Teqc os.system("rm %s/data.T04" % (year)) os.system("cp /net/daffy/mnt/data_drive/GNSSResults/%4d/%4d%02d%02d/%s_%4d%02d%02d.T04 %s/data.T04" % (year, year, month, day, RXData[RXNum][0], year, month, day, year)) os.system("cd %s ; ./doit.sh %d %s; cd .." % (year, year - 2000, RXData[RXNum][2])) os.system("mv %4d/dataDec-noGLN.%2dS %4d/%s-%4d%02d%02d-noGLN-L2E.MPX" % (year, year - 2000, year,RXData[RXNum][1], year, month, day)) os.system("mv %4d/dataDec-GLN_P.%2dS %4d/%s-%4d%02d%02d-GLN_P.MPX" % (year, year - 2000, year,RXData[RXNum][1], year, month, day)) os.system("mv %4d/dataDec-GLN_CA.%2dS %4d/%s-%4d%02d%02d-GLN_CA.MPX" % (year, year - 2000, year,RXData[RXNum][1], year, month, day)) os.system("mv %4d/dataDec-noGLN-L2C.%2dS %4d/%s-%4d%02d%02d-noGLN-L2C.MPX" % (year, year - 2000, year,RXData[RXNum][1], year, month, day)) # Cycle slip info os.system("mv %4d/dataDec-NoGLN.%2dslips %4d/%s-%4d%02d%02d-noGLN-L2E.Slips" % (year, year - 2000, year,RXData[RXNum][1], year, month, day)) os.system("mv %4d/dataDec-GLN_P.%2dslips %4d/%s-%4d%02d%02d-GLN_P.Slips" % (year, year - 2000, year,RXData[RXNum][1], year, month, day)) os.system("mv %4d/dataDec-GLN_CA.%2dslips %4d/%s-%4d%02d%02d-GLN_CA.Slips" % (year, year - 2000, year,RXData[RXNum][1], year, month, day)) os.system("mv %4d/dataDec-NoGLN-L2C.%2dslips %4d/%s-%4d%02d%02d-noGLN-L2C.Slips" % (year, year - 2000, year,RXData[RXNum][1], year, month, day)) ## YYYYMMDD-noGLN-L2E.Slips & .MPX ## - No GLONASS ## - No QZSS ## - No L2C ## ## YYYYMMDD-noGLN-L2C.Slips & .MPX ## - No GLONASS ## - No QZSS ## - No L2E ## ## YYYYMMDD-GLN-CA.Slips & .MPX ## - GLONASS Only ## - No GLONASS P (GLN CA present) ## ## YYYYMMDD-GLN-P.Slips & .MPX ## - GPS Only ## - No GLONASS CA (GLN P present) ## MPXX = [None] * 10 fileroot = "%4d/%s-%4d%02d%02d-noGLN-L2E" % (year,RXData[RXNum][1], year, month, day) resroot = "%4d/%s-noGLN-L2E" % (year,RXData[RXNum][1]) filename = "%s.MPX" % (fileroot) if(os.path.exists(filename)): parse17Sfile(fileroot,resroot) fileroot = "%4d/%s-%4d%02d%02d-noGLN-L2C" % (year,RXData[RXNum][1], year, month, day) resroot = "%4d/%s-noGLN-L2C" % (year,RXData[RXNum][1]) filename = "%s.MPX" % (fileroot) if(os.path.exists(filename)): parse17Sfile(fileroot,resroot) fileroot = "%4d/%s-%4d%02d%02d-GLN_P" % (year,RXData[RXNum][1], year, month, day) resroot = "%4d/%s-GLN_P" % (year,RXData[RXNum][1]) filename = "%s.MPX" % (fileroot) if(os.path.exists(filename)): parse17Sfile(fileroot,resroot) fileroot = "%4d/%s-%4d%02d%02d-GLN_CA" % (year,RXData[RXNum][1], year, month, day) resroot = "%4d/%s-GLN_CA" % (year,RXData[RXNum][1]) filename = "%s.MPX" % (fileroot) if(os.path.exists(filename)): parse17Sfile(fileroot,resroot) filename = "%4d/%s-noGLN-L2E-MP12.txt" % (year,RXData[RXNum][1]) plotMPxx(filename) filename = "%4d/%s-noGLN-L2C-MP12.txt" % (year,RXData[RXNum][1]) plotMPxx(filename) filename = "%4d/%s-GLN_P-MP12.txt" % (year,RXData[RXNum][1]) plotMPxx(filename) filename = "%4d/%s-GLN_CA-MP12.txt" % (year,RXData[RXNum][1]) plotMPxx(filename) # Plot the trend data for i in range(len(file2Str)): plotMPXTrend(list(file2Str.keys())[i])