""" Similar track plotting functionality to plt_tracking.py, but for multiple T04 files. The user specifies a space seperate list of files on the command line. Only the tracking (-t for plt_tracking) function is supported. The plots are all saved to PNG files in the current directory. If the file contains data for more than one antenna, data from both are plotted. """ import matplotlib from pylab import * import mutils as m import sys import numpy as np import sqlite3 import pandas as pd systemLUT = { 0:'GPS', 1:'SBAS', 2:'GLONASS', 3:'Galileo', 4:'QZSS', 9:'NavIC', 10:'BeiDou'} if __name__ == '__main__': # Allow running headless from the command line matplotlib.use("agg") # get the list of files from the command line fileList = sys.argv[1:] if(len(fileList) == 0): print("Usage: %s .T04 [.T04 .... .T04]" % sys.argv[0]) sys.exit colors = ['b','g','r','c','m','y','k'] # Create an empty pandas dataframe: df = pd.DataFrame(columns=['year', 'month', 'day', 'hour', 'sn', 'sat_type', 'freq', 'track', 'elMask', 'epochs', 'slips']) dataSetList = [] dataList = [] for thisFile in fileList: try: if(thisFile.endswith('.T04')): (data,key)=m.vd2arr(thisFile,rec='-d35:19') # The key should be consistent across all T04 files # Note - there could be a problem if the 35:19 schema # changes between files, i.e. due to different firmware if(len(data) > 0): # Extract the serial number and time from the filename tokens = thisFile.split('/') # Assumes SSSSSSSSSSYYYYMMDDHHMM.T04 serialNumber = tokens[-1][0:10] year = tokens[-1][10:14] month = tokens[-1][14:16] day = tokens[-1][16:18] hour = tokens[-1][18:20] timeStr = year + '-' + month + '-' + day + 'T' + hour + ':00:00' i = m.find( data[:,key.ANTENNA] == 0) j = m.find( data[:,key.ANTENNA] == 1) if( (len(i) > 0) and (len(j) > 0) ): # We have a dual antenna dataset dataList.append(data[i,:]) dataList.append(data[j,:]) dataSetList.append(serialNumber + ' Antenna 0') dataSetList.append(serialNumber + ' Antenna 1') elif( len(j) > 0): dataList.append(data[j,:]) # Unusual - only data from the second antenna dataSetList.append(serialNumber+ ' Antenna 1') else: dataList.append(data[i,:]) # Single antenna data dataSetList.append(serialNumber) except: print("Error processing %s" % thisFile) numDataSets = len(dataList) # We have the data # Empty numpy array to hold the unique satellite types sat_types = np.array([]) for thisData in dataList: sats_this_file = np.unique( thisData[:,key.SAT_TYPE] ) # Append the unique satellite types to the array sat_types = np.append(sat_types,sats_this_file) # Get a sorted unit array of the satellite types sat_types = np.sort(np.unique(sat_types)).astype(int) for thisSatType in sat_types: thisSatData = [] sig_list = [] sv_ids = np.array([]) for thisData in dataList: i = m.find(thisData[:,key.SAT_TYPE] == thisSatType) if(len(i) > 0): thisSatData.append(thisData[i,:]) sig_tuple = m.get_signals( thisData[i], key )[thisSatType] for thisSig in sig_tuple: sig_list.append(thisSig) svs_this_file = np.unique( thisData[i,key.SV] ) sv_ids = np.append(sv_ids,svs_this_file) else: # Need to maintain the data order, so add an empty array thisSatData.append([]) sv_ids = np.sort(np.unique(sv_ids)).astype(int) sig_list = list(set(sig_list)) sig_list.sort() for sig_id,(freq,track) in enumerate(sig_list): figure() for set_num,thisData in enumerate(thisSatData): if(len(thisData) == 0): continue i = m.find( (thisData[:,key.FREQ] == freq) & (thisData[:,key.TRACK] == track) ) if len(i) == 0: continue first = [False] * numDataSets for svNum, sv in enumerate(sv_ids): j = m.find( thisData[i,key.SV] == sv ) if len(j) == 0: continue iSlip = m.find( (thisData[i[j],key.MEAS].astype(int) & key.dMEAS_SLIP)!=0 ) xData = thisData[i[j],key.TIME].copy() # insert NaN for time gaps. That way we get gaps for tracking outages. dt = abs(np.diff(xData)) expected_dt = np.median(dt)+.001 igap = m.find( dt > expected_dt )+1 if len(igap) > 0: sv_xgap = np.insert(xData, igap, np.nan ) else: sv_xgap = xData # Create a y-axis data set for the plot. Generate an array # of the same length as the x-axis data, but with all values # set to the satellite ID # numDataSets + 1 so we have a space between each satellite sv_ygap = np.ones(len(sv_xgap))*svNum + set_num/(numDataSets+1) if(first[set_num] == False): first[set_num] = True label = dataSetList[set_num] else: label='_nolegend_' # plot tracking line plot( sv_xgap, sv_ygap, '.-', markersize=1, color=colors[set_num], label=label) # plot cycle slips plot( thisData[i[j[iSlip]],key.TIME], sv_ygap[0:len(iSlip)], 'o', mfc='none', label='_nolegend_', color=colors[set_num], markersize=3) allSlips = m.find( (thisData[i,key.MEAS].astype(int) & key.dMEAS_SLIP)!=0 ) # Add this data to the dataframe df = df.append(pd.DataFrame([[int(year), int(month), int(day), int(hour), dataSetList[set_num], # Serial number thisSatType, # current satellite type freq, track, # data type np.nan, # elevation mask - not available len(i), # number of epochs len(allSlips) # number of cycle slips ]], columns=df.columns), ignore_index=True) # Now create the metrics for an elevation mask of 10 degrees deg10Index = m.find( thisData[i,key.EL] >= 10 ) allSlips = m.find( (thisData[i[deg10Index],key.MEAS].astype(int) & key.dMEAS_SLIP)!=0 ) # Add this data to the dataframe df = df.append(pd.DataFrame([[int(year), int(month), int(day), int(hour), dataSetList[set_num], # Serial number thisSatType, # current satellite type freq, track, # data type 10, # elevation mask len(deg10Index), # number of epochs len(allSlips) # number of cycle slips ]], columns=df.columns), ignore_index=True) print('plotted %d %s %s' % (set_num, dataSetList[set_num], m.get_sub_type(thisSatType,freq,track).fullstr)) typeStr = m.get_sub_type(thisSatType,freq,track).fullstr # Set the legend font to 8pixels legend(fontsize=6) title(timeStr + ' - ' + typeStr) yticks( arange(len(sv_ids)), sv_ids.astype(int) ) grid(True) ylabel('SV ID') xlabel('Time [GPS Secs]') tight_layout() # Now save the plot # replace spaces with '-' and add .png to the data type string # use this as the filename filename = year + '-' + month + '-' + day + 'T' + hour + 'H-' filename += typeStr.replace(' ','-') + '.png' savefig(filename,dpi=600) close() print(df) # Write the DataFrame to a sqlite table called 'data' db = sqlite3.connect('trackData.db') df.to_sql('data', db, if_exists='append', index=False) # Remove duplicate entries from the SQL table db.execute(''' DELETE FROM data WHERE rowid NOT IN ( SELECT MIN(rowid) FROM data GROUP BY year, month, day, hour, sn, freq, track, elMask ) ''') db.commit() db.close() #db = sqlite3.connect('trackData.db') # Now read the table back into a DataFrame and sort it #df2 = pd.read_sql_query('SELECT * FROM data ORDER BY year, month, day, hour, freq, track, elMask, sn', db) #print('Data written to trackData.db')