#!/usr/bin/env python # # analyzeTracked.py # # Desc: The user provides this script with a T04 containing observables, the # script will load the obs record from the T04, decimating on import. # The start/stop time will be extracted and RINEX for that day (it does # not handle a file that straddles a day) will be downloaded and imported. # We have a two level cache (1) the RINEX data is stored locally, if the # file exists it is used (2) we use GeoRINEX to import the RINEX data, it # has a faster file format for re-import of the data, we save this faster # format and use it as the perferred data source if it exists. # # The RINEX is an aggregation of all broadcast data. # # Once we have the data and RINEX loaded we scan a system at a time # and compute the expected vs actual for each satellite / code type. We have # a dictionary that provides a LUT of which satellites support which frequeny/code. # The script is aware of different level of support of the different GNSS satellites. # We use a hard table, we could convert to e JSON file, and this information may # be something we could extract from the RINEX??? # # # Limitations: # Only GPS has been implemented. # # ToDo: # - Extend the GNSSSignals table # - Check the ephemeris routine - during an early test this was within a meter # of the C/C++ firmware version # - Add BDS GEO support in the ephemeris calculation # - GLONASS uses a different ephemeris format - GeoRINEX can read it - implement the ephemeris calc. # - Add per system physical parameters to the ephemeris routines - may have negligible impact # - We look at three sources of RINEX, #2 & #3 both have issues (#2 is GPS, #3 may be EU only). only # #1 is global and all systems (including IRNSS). We should: # - Look for other sources (data is available in RINEX 4.0 - it may be from the same DLR source and # GeoRINEX does not support 4.0 yet) # - Instead of selecting a source, read them all in and merge. # nav = xarray.merge((nav1, nav2)) # - the code spends a lot of time in the ephemeris calc, if we have a high elevation, we know that # we won't set quickly and can skip the calculation (only needed to get the elevation angle to check # against the mask). A crude low elevation cut off has been added, we can optimize this and call the # ephemeris calculation much less frequently. # - Data is printed to stdout - append to a summary file or export in a structured way # # Copyright Trimble Inc 2022 # import math import mutils as m import datetime import numpy as np import argparse import rinSupport as rin from colorama import Fore, Style ########################## # Satellite / Signal LUT ########################## # This will change over time, we could add date fields when the satellites change, we # could also make this a loadable JSON file. GNSSSignals = {'E':{'E1':{'svs':[ 1, 2, 3, 4, 5, 7, 8, 9,10, 11,12,13,14,15, 18,19, 21, 24,25,26,27, 30, 31, 33,34, 36], 'T04Freq':0, 'T04Code':23}, 'E5A':{'svs':[ 1, 2, 3, 4, 5, 7, 8, 9,10, 11,12,13,14,15, 18,19, 21, 24,25,26,27, 30, 31, 33,34, 36], 'T04Freq':2, 'T04Code':11}, 'E5B':{'svs':[ 1, 2, 3, 4, 5, 7, 8, 9,10, 11,12,13,14,15, 18,19, 21, 24,25,26,27, 30, 31, 33,34, 36], 'T04Freq':3, 'T04Code':11}, 'E5Alt':{'svs':[ 1, 2, 3, 4, 5, 7, 8, 9,10, 11,12,13,14,15, 18,19, 21, 24,25,26,27, 30, 31, 33,34, 36], 'T04Freq':4, 'T04Code':14}, 'E6':{'svs':[ 1, 2, 3, 4, 5, 7, 8, 9,10, 11,12,13,14,15, 18,19, 21, 24,25,26,27, 30, 31, 33,34, 36], 'T04Freq':5, 'T04Code':36} }, 'G':{'L1CA':{'svs':[ 1, 2, 3, 4, 5, 6, 7, 8, 9,10, 11,12,13,14,15,16,17,18,19,20, 21,22,23,24,25,26,27,28,29,30, 31,32], 'T04Freq':0, 'T04Code':0}, 'L1C':{'svs':[ 4, 11, 14, 18, 23], 'T04Freq':0, 'T04Code':20}, 'L2E':{'svs':[ 1, 2, 3, 4, 5, 6, 7, 8, 9,10, 11,12,13,14,15,16,17,18,19,20, 21,22,23,24,25,26,27,28,29,30, 31,32], 'T04Freq':1, 'T04Code':2}, 'L2C':{'svs':[1, 3, 4, 5, 6, 7, 8, 10, 11,12, 14,15, 17,18, 23,24,25,26,27,28,29,30, 31,32], 'T04Freq':1, 'T04Code':5}, # Make this a list to cover CM, CL, and CM+CL? 'L5':{'svs':[1, 3, 4, 6, 8, 10, 11, 14, 18, 23,24,25,26,27,28, 30, 32], 'T04Freq':2, 'T04Code':8} }, # Make this a list to cover D / P / D+P 'J':{'L1CA':{'svs':[ 1, 2, 3, 4, 7 ], 'T04Freq':0, 'T04Code':0}, 'L1C':{'svs':[ 1, 2, 3, 4, 7 ], 'T04Freq':0, 'T04Code':20}, 'L2C':{'svs':[ 1, 2, 3, 4, 7 ], 'T04Freq':1, 'T04Code':5}, 'L5':{'svs':[ 1, 2, 3, 4, 7 ], 'T04Freq':2, 'T04Code':8}}, 'I':{'L5CA':{'svs':[ 2, 3, 4, 5, 6, 7, 9 ], 'T04Freq':2, 'T04Code':0}, 'S1CA':{'svs':[ 2, 3, 4, 5, 6, 7, 9 ], 'T04Freq':11, 'T04Code':0}}, 'C':{'B1I':{'svs':[ 1, 2, 3, 4, 5, 6, 7, 8, 9,10, 11,12,13,14, 16, 19,20, 21,22,23,24,25,26,27,28,29,30, 31,32,33,34,35,36,37,38,39,40, 41,42,43,44,45,46, 56,57,58,59,60, 61 ], 'T04Freq':6, 'T04Code':26}, 'B1C':{'svs':[ 19,20, 21,22,23,24,25,26,27,28,29,30, 31,32,33,34,35,36,37,38,39,40, 41,42,43,44,45,46, 56,57,58], 'T04Freq':0, 'T04Code':20}, 'B2I':{'svs':[ 1, 2, 3, 4, 5, 6, 7, 8, 9,10, 11,12,13,14, 16], # BDS-II 'T04Freq':3, 'T04Code':28}, 'B2A':{'svs':[ 19,20, 21,22,23,24,25,26,27,28,29,30, 31,32,33,34,35,36,37,38,39,40, 41,42,43,44,45,46, 56,57,58], 'T04Freq':2, 'T04Code':8}, 'B2B':{'svs':[ 19,20, 21,22,23,24,25,26,27,28,29,30, 31,32,33,34,35,36,37,38,39,40, 41,42,43,44,45,46, 56,57,58,59,60, 61 ], 'T04Freq':3, 'T04Code':6}, 'B3I':{'svs':[ 1, 2, 3, 4, 5, 6, 7, 8, 9,10, 11,12,13,14, 16, 19,20, 21,22,23,24,25,26,27,28,29,30, 31,32,33,34,35,36,37,38,39,40, 41,42,43,44,45,46, 56,57,58,59,60, 61 ], 'T04Freq':7, 'T04Code':29} } } # The signals we will process and the RINEX / T04 mapping sattype = [ {'RINEX':'G', 'T04':0}, # GPS {'RINEX':'E', 'T04':3}, # Galileo {'RINEX':'C', 'T04':10}, # BeiDou {'RINEX':'J', 'T04':4}, # QZSS {'RINEX':'I', 'T04':9}] # IRNSS # Caller must provide the year, day of year, start and stop seconds over which we'll # analyze data. decRate is the decimation rate of the analysis and input data in # seconds. data and key are record 35:19 data, elevMask is the analysis elevation mask. # Finally, provide the user position we'll use to compute the elevation angle. def getTracked(year,dayofYear,start,stop,decRate, data,key, xUser,yUser,zUser, elevMask=[10], only_sattypes=None): results = [] if(isinstance(elevMask,list) == False): # Old script did not use a list, force to a list elevMask = [elevMask] # Script assumes elevMask goes from low to high elevMask.sort() print('Output yield for elevMask:',elevMask[0]) # loop around for each system, then loop on the satellite ID (found in the RINEX file), and # the most inner loop iterates from the start to the stop time incrementing by the decimation # rate. for thisSystem in sattype: thisType = thisSystem['RINEX'] if only_sattypes is not None and only_sattypes.find(thisType) < 0: continue print('Loading RINEX: ' + thisType) nav = rin.getRINEXNav(year,dayofYear,thisType) print('Got RINEX: ' + thisType) sigKeys = list(GNSSSignals[thisType].keys()) if(len(nav)>0): for thisSV in nav['sv'].data: try: # For some systems, e.g. Galileo, we may get ephemeris from # multiple sources (E1, E5A, ...). Georinex adds a "_" to the # extra ephemeris sources, e.g. E1_1, E1_2. Remove the secondary # ephemeris and use only the primary (e.g. E1). if('_' in thisSV): continue # RINEX QZSS numbers from 1, PRNs start at 193 - adjust so that # the T04 data and RINEX agree sv = int(thisSV[1:]) if(thisType == 'J'): sv+=192 # Get data for this satellite from the primary antenna i = m.find( (data[:,key.ANTENNA] == 0) & (data[:,key.SAT_TYPE]==thisSystem['T04']) & (data[:,key.SV] == sv) ) if(len(i) == 0): svData = data[i,:] else: cols = np.array([key.TIME,key.FREQ,key.TRACK,key.MEAS]) svData = data[i,:] svData = svData[:,cols] # Remove this data so the primary array is smaller next time around #data = np.delete(data, i, axis = 0) # We've selected just a subset of the data to reduce the memory footprint timeKey = 0 freqKey = 1 trackKey = 2 measKey = 3 expectCount = [] actualCount = [] actualCountPLL = [] for index in range(len(elevMask)): expectCount.append([0] * len(sigKeys)) actualCount.append([0] * len(sigKeys)) actualCountPLL.append([0] * len(sigKeys)) # ToDo - we could speed up this loop. If we know that the satellite has a large # magnitude negative elevation angle, we can jump further forward than the # decRate! For now a quick hack, if the elev is < (elevMask-20) deg add # 1,000 seconds. Need to analyze the maximum elevation angle rate and adjust thisTime = start deltaT = 1000000 # Force collection of a new ephemeris elev = -1 computeT = 1000000 removeDataCount = 0 #for thisTime in range(start,stop,decRate): while(thisTime <= stop): # We only need to search for a new ephemeris if the current one is # getting stale if(deltaT >= (3*3600)): # This is an expensive function call [eph,deltaT] = rin.getEph(nav, thisSV, thisTime) # If deltaT is still large, we don't have data for the satellite. # Jump forward in time to avoid trying to load an ephemeris # every decRate seconds if(deltaT > (3*3600)): thisTime += deltaT - 3*3600 - decRate if(deltaT < (3*3600)): # If we are at a higher elevation, we can reduce the rate # at which we compute the elevation angle as the satellite # will take a while to dip below the mask #if(True): if( (elev < (elevMask[-1] + 20)) or (computeT >= 600)): elev = rin.getElev(eph,thisSV,thisTime,xUser,yUser,zUser) elev *= 180/math.pi computeT = 0 if(elev < (elevMask[0]-20)): # Low elevation (relative to the mask) satellite - bump by 1000s # ToDo - should be able to jump by a larger time thisTime += 1000 deltaT += 1000 # Keep this in sync so we grab a new ephemeris when needed else: for elevIndex in range(len(elevMask)): thisElev = elevMask[elevIndex] if(elev > thisElev): # Check whether we tracked this satellite # ToDo check key.FREQ and key.TRACK and find which signals # are tracked. Compare to what we expect in GNSSSignals. for thisKeyIndex in range(len(sigKeys)): expectCount[elevIndex][thisKeyIndex] += 1 # Get track data for this timestamp i = m.find(svData[:,timeKey] == thisTime) if(len(i) > 0): for thisKeyIndex in range(len(sigKeys)): thisKey = sigKeys[thisKeyIndex] gotKey = False for thisIndex in i: if( (svData[thisIndex,freqKey] == GNSSSignals[thisType][thisKey]['T04Freq']) and (svData[thisIndex,trackKey] == GNSSSignals[thisType][thisKey]['T04Code']) ): actualCount[elevIndex][thisKeyIndex] += 1 # Now check if this data is phase locked if(int(svData[thisIndex,measKey]) & 1): # Test phase valid bit actualCountPLL[elevIndex][thisKeyIndex] += 1 gotKey = True break if(gotKey == True): break thisTime += decRate deltaT += decRate computeT += decRate # We've looped around for the current satellite, provide a summary # Collate all data, but only print to stdout data for the lowest # elevation satellite if(sum(expectCount[0]) > 0): # we should have some data - need to check whether this satellite supports # each signal. print(thisSV + ': ',end='') for thisKeyIndex in range(len(sigKeys)): thisKey = sigKeys[thisKeyIndex] thisSatDataType = GNSSSignals[thisType][thisKey] for sv in thisSatDataType['svs']: if(sv == int(thisSV[1:])): percent = 100.0*actualCount[0][thisKeyIndex]/expectCount[0][thisKeyIndex] if(percent == 100): color = Fore.GREEN else: color = Fore.RED print("%5s %s%5.1f%s " % (thisKey,color,percent,Style.RESET_ALL),end='') # actualCount - all epochs (FLL or PLL) # actualCountPLL - must have phase valid (PLL) set # expected - based on the ephemeris (from RINEX) number of epochs we should get # with perfect tracking. for elevIndex in range(len(elevMask)): thisElevStr = str(elevMask[elevIndex]) epochs = {'signal':thisKey, 'actual_' + thisElevStr:actualCount[elevIndex][thisKeyIndex], 'actualPLL_' + thisElevStr:actualCountPLL[elevIndex][thisKeyIndex], 'expected_' + thisElevStr:expectCount[elevIndex][thisKeyIndex]} results.append({thisSV:epochs}) print("") except: print("Problem with:",thisSV) print('RINEX complete') return(results) if __name__ == "__main__": ###################################################################### # Parse arguments parser = argparse.ArgumentParser(description='Analyzes observable yield against public RINEX NAV (ephemeris) data\n' +' The first record 35:2 is used by default as the station position') parser.add_argument('filename', help='T04 filename') parser.add_argument('-d','--decimation', help='Observable decimation rate in seconds (default = 10s)') parser.add_argument('-e','--elevMask', help='Elevation Mask (default = 10 degrees or multiple "-e 5,10,15")') parser.add_argument('-s','--station_xyz', help='Station X,Y,Z - default (first T04 35:2 record)') parser.add_argument('-l','--station_llh', help='Station Lat,Lon,Hgt - default (first T04 35:2 record)') parser.add_argument('-b','--station_default', help='When set, use a Sunnyvale default position', action="store_true") parser.add_argument('-o','--only_sattypes', help='list of RINEX satellite types to process (e.g.: G,I,E)') args = parser.parse_args() ###################################################################### T04File = args.filename if(args.decimation): decRate = int(args.decimation) else: decRate = 10 if(args.elevMask): elevStr = args.elevMask.split(',') elevMask = [int(thisStr) for thisStr in elevStr] # Code assumes the elevation mask is in assending order elevMask.sort() else: elevMask = [10] if(args.station_xyz): xyz_str = args.station_xyz [xUser,yUser,zUser] = np.array([float(x) for x in xyz_str.split(',')]) elif(args.station_llh): llh_str = args.station_llh [Lat,Lon,Hgt] = np.array([float(x) for x in llh_str.split(',')]) [xUser,yUser,zUser] = m.conv_lla_ecef( Lat, Lon, Hgt) elif(args.station_default): # Sunnyvale station position xUser = -2689303.1536 yUser = -4302871.9778 zUser = 3851431.0448 else: # By default use the position from the T04 file - this is a change to an earlier # version of this script [Lat,Lon,Hgt] = m.getRefPosition(T04File) [xUser,yUser,zUser] = m.conv_lla_ecef( Lat, Lon, Hgt) print('File:',T04File) print('Decimation Rate:',decRate,'s') print('Elev Mask:',elevMask,'deg') print('User Position (XYZ):',xUser,yUser,zUser) [Lat,Lon,Hgt] = m.conv_ecef_lla(xUser,yUser,zUser) print('User Position (LLH): %.10f %.10f %.3f' % (Lat,Lon,Hgt)) # Load the T04 file decimating during import (data,key)=m.vd2arr(T04File,rec='-d35:19',opt=['--dec=' + str(decRate * 1000)]) # Find the start/stop window start = data[0,key.TIME] stop = data[-1,key.TIME] if(start > stop): print("ToDo - week wrap") exit() # Make sure we are on the decimation rate bound start = (int(start/decRate) + 1)*decRate stop = int(stop/decRate)*decRate # Find the date of the test - use this to get the correct RINEX file week = data[0,key.WN] # ToDo - leap seconds? start_time = datetime.datetime(1980, 1, 6) + datetime.timedelta(seconds = (start + 7*86400*week)) year = start_time.year month = start_time.month day = start_time.day # date to day within the year - RINEX filenames need this dayofYear = m.date_to_doy(year,month,day) print('Year: %d\nMonth:%d\nDay: %d\nDOY: %d' % (year,month,day,dayofYear)) print('WN: %d\nStartTime:%d\nStopTime: %d' % (week,start,stop)) res = getTracked(year,dayofYear,start,stop,decRate,data,key,xUser,yUser,zUser,elevMask,args.only_sattypes)