Skip to content

Commit

Permalink
- fixed sign of BCNAV1
Browse files Browse the repository at this point in the history 
- added STO/EOP/ION support for RINEX 4.02
  • Loading branch information
@hirokawa
hirokawa committed Jan 18, 2025
1 parent 423014b commit e8689ee
Show file tree
Hide file tree
Showing 3 changed files with 154 additions and 57 deletions.
34 changes: 32 additions & 2 deletions src/cssrlib/gnss.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -616,6 +616,31 @@ def __gt__(self, other):
(self.time == other.time and self.sec > other.sec)


class STOParam():
""" System Time and UTC Office """
sbas = 0 # SBAS ID
prm = [0, 0] # System time offset parameter
t_ot = None # reference epoch
t_t = 0.0 # transmission time of message (Time of week [sec])
a = np.zeros(3) # a0, a1, a2


class EOPParam():
""" Earth Orientation Parameter """
prm = np.zeros(9)
# EOP parameters (xp,dxp,ddxp,yp,dyp,ddyp,dut1,ddut1,dddut1)
t_ot = None # reference epoch
t_t = 0.0 # transmission time of message (Time of week [sec])


class IONParam():
""" Ionospheric delay model Parameter """
iod = 0
prm = np.zeros(9) # ION parameters
t_tm = None # transmission time
region = None


class Obs():
""" class to define the observation """

Expand Down Expand Up @@ -787,8 +812,13 @@ def __init__(self, nf=2):
[0.1167E+06, -0.2294E+06, -0.1311E+06, 0.1049E+07]])
self.ion_gim = np.zeros(9)
self.ion_region = 0 # 0: wide-area, 1: Japan-aera (QZSS only)
self.sto = np.zeros(3)
self.sto_prm = np.zeros(4, dtype=int)
# self.sto = np.zeros(3)
# self.sto_prm = np.zeros(4, dtype=int)

self.sto_prm = {}
self.eop_prm = {}
self.ion_prm = {}

self.eop = np.zeros(9)
self.elmin = np.deg2rad(15.0)
self.tidecorr = False
Expand Down
32 changes: 22 additions & 10 deletions src/cssrlib/rawnav.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -142,10 +142,14 @@ def getbitg(self, buff, pos, len_):

def urai2sva(self, urai, sys=uGNSS.GPS):
""" GPS/QZSS SV accuracy [1] 20.3.3.3.1.3, [2] Tab. 5.4.3-2 """
urai_t = [2.40, 3.40, 4.85, 6.85, 9.65, 13.65, 24.00, 48.00, 96.00,
192.00, 384.00, 768.00, 1536.00, 3072.00, 6144.00, -1.0]
# sva_t = [2.40, 3.40, 4.85, 6.85, 9.65, 13.65, 24.00, 48.00, 96.00,
# 192.00, 384.00, 768.00, 1536.00, 3072.00, 6144.00, -1.0]

return urai_t[urai]
# from RINEX 4.02
sva_t = [2.0, 2.8, 4.0, 5.7, 8.0, 11.3, 16.0, 32.0,
64.0, 128.0, 256.0, 512.0, 1024.0, 2048.0, 4096.0, 8192.0]

return sva_t[urai]

def sisa2sva(self, sisa):
""" Galileo SISA Index Values [3] Tab. 89 """
Expand Down Expand Up @@ -873,26 +877,32 @@ def decode_bds_b1c(self, week, time_, sat, msg):

# decode Subframe 3
i = 608
page, eph.svh, eph.integ, eph.sismai = bs.unpack_from('u6u2u3u4',
page, eph.svh, eph.integ, eph.sismai = bs.unpack_from('u6u2u3s4',
msg, i)
i += 15
if page == 1: # Iono, BDT-UTC
eph.sisai[0] = bs.unpack_from('u5', msg, i)[0]
eph.sisai[0] = bs.unpack_from('u5', msg, i)[0] # SISAI_oe
i += 5
i = self.decode_bds_cnav_sisai(msg, eph, i)
elif page == 2: # Reduced almanac
i = self.decode_bds_cnav_sisai(msg, eph, i)
i += 22
elif page == 3: # EOP, BGTO
eph.sisai[0] = bs.unpack_from('u5', msg, i)[0]
eph.sisai[0] = bs.unpack_from('u5', msg, i)[0] # SISAI_oe
i += 5
# EOP
# BGTO
elif page == 4: # Midi almanac
i = self.decode_bds_cnav_sisai(msg, eph, i)
i += 22
# Midi Almanac

# i += 264-15
eph.mode = 1 # B-CNAV1

if page != 1: # page 1 include SISAI*
return None

return eph

def decode_bds_b2a(self, week, time_, sat, msg):
Expand All @@ -918,7 +928,7 @@ def decode_bds_b2a(self, week, time_, sat, msg):
id4, sow4 = bs.unpack_from('u6u18', buff, 320*3+6)
id5, sow5 = bs.unpack_from('u6u18', buff, 320*4+6)

if id1 != 10 or id2 != 11 or id3 != 30 or id4 != 34 or id5 != 40:
if id1 != 10 or id2 != 11 or id3 != 30 or id5 != 40:
return None

if sow2 != sow1+1 or sow3 != sow2+1:
Expand All @@ -936,11 +946,13 @@ def decode_bds_b2a(self, week, time_, sat, msg):
if id5 == 40:
i = 320*4+42
eph.sisai[0] = bs.unpack_from('u5', buff, i)[0]
i += 5
i = self.decode_bds_cnav_sisai(buff, eph, i)

# decode MT10
i = 30
eph.week, ib2a, eph.sismai, ib1c, eph.iode = bs.unpack_from(
'u13u3u4u3u8', buff, i)
'u13u3s4u3u8', buff, i)
i += 31
i = self.decode_bds_cnav_eph1(buff, eph, i)
eph.integ = (ib2a << 3)+ib1c
Expand All @@ -962,7 +974,7 @@ def decode_bds_b2a(self, week, time_, sat, msg):
return None

i = self.decode_bds_cnav_iono(buff, i)
tgd_b1cp = bs.unpack_from('u12', buff, i)[0]
tgd_b1cp = bs.unpack_from('s12', buff, i)[0]

eph.tgd = tgd_b1cp*rCST.P2_34
eph.isc[1] = isc_b2ad*rCST.P2_34
Expand Down Expand Up @@ -1009,7 +1021,7 @@ def decode_bds_b2b(self, week, time_, sat, msg, ofst=12):
i = self.decode_bds_cnav_eph1(buff, eph, i)
i = self.decode_bds_cnav_eph2(buff, eph, i)

eph.integ, eph.sismai = bs.unpack_from('u3u4', buff, i)
eph.integ, eph.sismai = bs.unpack_from('u3s4', buff, i)
i += 7

# decode MT30
Expand Down
145 changes: 100 additions & 45 deletions src/cssrlib/rinex.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,8 @@
from cssrlib.gnss import gpst2time, bdt2time, epoch2time, timediff, gtime_t
from cssrlib.gnss import prn2sat, char2sys, timeget, utc2gpst, time2epoch
from cssrlib.gnss import Eph, Geph, Obs, sat2id, sat2prn, gpst2bdt, time2gpst
from cssrlib.gnss import timeadd, id2sat, gpst2utc, Seph
from cssrlib.gnss import timeadd, id2sat, gpst2utc, Seph, STOParam, EOPParam
from cssrlib.gnss import IONParam


class pclk_t:
Expand Down Expand Up @@ -49,6 +50,13 @@ def __init__(self):
self.mode_nav = 0
self.glo_ch = {}

self.ofst_src = {'GP': uGNSS.GPS, 'GL': uGNSS.GLO,
'GA': uGNSS.GAL, 'BD': uGNSS.BDS,
'QZ': uGNSS.QZS, 'IR': uGNSS.IRN,
'SB': uGNSS.SBS, 'UT': uGNSS.NONE}
self.itype_t = {'LNAV': 0, 'FDMA': 1, 'IFNV': 2, 'D1D2': 3,
'SBAS': 4, 'CNVX': 5, 'L1NV': 6, 'LXOC': 7}

def setSignals(self, sigList):
""" define the signal list for each constellation """

Expand Down Expand Up @@ -172,111 +180,155 @@ def decode_nav(self, navfile, nav, append=False):
if self.ver >= 4.0:

if line[0:5] == '> STO': # system time offset (TBD)
ofst_src = {'GP': uGNSS.GPS, 'GL': uGNSS.GLO,
'GA': uGNSS.GAL, 'BD': uGNSS.BDS,
'QZ': uGNSS.QZS, 'IR': uGNSS.IRN,
'SB': uGNSS.SBS, 'UT': uGNSS.NONE}

sys = char2sys(line[6])
itype = line[10:14]

if sys not in nav.sto_prm:
nav.sto_prm[sys] = {}

if itype not in self.itype_t:
fnav.readline()
fnav.readline()
continue

im = self.itype_t[itype]

if im not in nav.sto_prm[sys]:
nav.sto_prm[sys][im] = STOParam()

line = fnav.readline()
ttm = self.decode_time(line, 4)
nav.sto_prm[sys][im].t_ot = self.decode_time(line, 4)
mode = line[24:28]
if mode[0:2] in ofst_src and mode[2:4] in ofst_src:
nav.sto_prm[0] = ofst_src[mode[0:2]]
nav.sto_prm[1] = ofst_src[mode[2:4]]
if mode[0:2] in self.ofst_src and \
mode[2:4] in self.ofst_src:
nav.sto_prm[sys][im].prm[0] = \
self.ofst_src[mode[0:2]]
nav.sto_prm[sys][im].prm[1] = \
self.ofst_src[mode[2:4]]

line = fnav.readline()
ttm = self.flt(line, 0)
nav.sto_prm[sys][im].t_t = self.flt(line, 0)
for k in range(3):
nav.sto[k] = self.flt(line, k+1)
nav.sto_prm[sys][im].a[k] = self.flt(line, k+1)
continue

elif line[0:5] == '> EOP': # earth orientation param
sys = char2sys(line[6])
itype = line[10:14]

if sys not in nav.eop_prm:
nav.eop_prm[sys] = {}

if itype not in self.itype_t:
fnav.readline()
fnav.readline()
fnav.readline()
continue

im = self.itype_t[itype]

if im not in nav.eop_prm[sys]:
nav.eop_prm[sys][im] = EOPParam()

line = fnav.readline()
ttm = self.decode_time(line, 4)
nav.eop_prm[sys][im].t_eop = self.decode_time(line, 4)
for k in range(3):
nav.eop[k] = self.flt(line, k+1)
nav.eop_prm[sys][im].prm[k] = self.flt(line, k+1)
line = fnav.readline()
for k in range(3):
nav.eop[k+3] = self.flt(line, k+1)
nav.eop_prm[sys][im].prm[k+3] = self.flt(line, k+1)
line = fnav.readline()
ttm = self.flt(line, 0)
nav.eop_prm[sys][im].t_t = self.flt(line, 0)
for k in range(3):
nav.eop[k+6] = self.flt(line, k+1)
nav.eop_prm[sys][im].prm[k+6] = self.flt(line, k+1)
continue

elif line[0:5] == '> ION': # iono (TBD)
sys = char2sys(line[6])
itype = line[10:14]
stype = '' if len(line) < 20 else line[15:19]

if sys not in nav.ion_prm:
nav.ion_prm[sys] = {}

im = self.itype_t[itype]
nav.ion_prm[sys][im] = IONParam()

line = fnav.readline()
ttm = self.decode_time(line, 4)
nav.ion_prm[sys][im].t_tm = self.decode_time(line, 4)
if sys == uGNSS.GAL and itype == 'IFNV': # Nequick-G
for k in range(3): # ai0, ai1, ai2
nav.ion[0, k] = self.flt(line, k+1)
nav.ion_prm[sys][im].prm[k] = \
self.flt(line, k+1)
line = fnav.readline()
# disturbance flags
nav.ion[0, 3] = int(self.flt(line, 0))
nav.ion_prm[sys][im].prm[3] = \
int(self.flt(line, 0))
elif sys == uGNSS.BDS and itype == 'CNVX': # BDGIM
ttm = self.decode_time(line, 4)
self.ion_gim = np.zeros(9)
for k in range(3):
nav.ion_gim[k] = self.flt(line, k+1)
nav.ion_prm[sys][im].prm[k] = \
self.flt(line, k+1)
line = fnav.readline()
for k in range(4):
nav.ion_gim[k+3] = self.flt(line, k)
nav.ion_prm[sys][im].prm[k+3] = \
self.flt(line, k)
line = fnav.readline()
for k in range(2):
nav.ion_gim[k+7] = self.flt(line, k)
nav.ion_prm[sys][im].prm[k+7] = \
self.flt(line, k)
elif sys == uGNSS.IRN and itype == 'L1NV': # L1NAV
if stype == 'KLOB': #
iodk = self.flt(line, 1)
line = fnav.readline()
for k in range(4):
nav.ion[0, k] = self.flt(line, k)
nav.ion_prm[sys][im].prm[k] = \
self.flt(line, k)
line = fnav.readline()
for k in range(4):
nav.ion[1, k] = self.flt(line, k)
nav.ion_prm[sys][im].prm[k+4] = \
self.flt(line, k)
nav.ion_prm[sys][im].iod = iodk
line = fnav.readline()
nav.ion_region = np.zeros(4)
nav.ion_prm[sys][im].region = np.zeros(4)
for k in range(4):
nav.ion_region[k] = self.flt(line, k)

nav.ion_prm[sys][im].region[k] = \
self.flt(line, k)
elif stype == 'NEQN':
nav.ion_region = np.zeros((3, 4))
iodn = self.flt(line, 1)

nav.ion_prm[sys][im].iod = self.flt(line, 1)
prm = np.zeros((3, 8))
for j in range(3):
line = fnav.readline()
nav.ion = np.zeros((3, 4))
for k in range(4): # a0, a1, a2, idf
nav.ion[j, k] = self.flt(line, k)
prm[j, k] = self.flt(line, k)
line = fnav.readline()
# lon_min, lon_max, mopid_min, mopid_max
for k in range(4):
nav.ion_region[j, k] = \
self.flt(line, k)
prm[j, k+4] = self.flt(line, k)
nav.ion_prm[sys][im].prm = prm

elif sys == uGNSS.GLO and itype == 'LXOC':
c_A = self.flt(line, 1)
c_F10_7 = self.flt(line, 2)
c_Ap = self.flt(line, 3)
nav.ion[0, 0:3] = [c_A, c_F10_7, c_Ap]
nav.ion_prm[sys][im].prm[0:3] = \
[c_A, c_F10_7, c_Ap]

else: # Klobucher (LNAV, D1D2, CNVX)
self.ion_gim = np.zeros(9)
nav.ion_prm[sys][im].prm = np.zeros(9)

for k in range(3):
nav.ion[0, k] = self.flt(line, k+1)
nav.ion_prm[sys][im].prm[k] = \
self.flt(line, k+1)
line = fnav.readline()
nav.ion[0, 3] = self.flt(line, 0)
for k in range(3):
nav.ion[1, k] = self.flt(line, k+1)
for k in range(4):
nav.ion_prm[sys][im].prm[k+3] = \
self.flt(line, k)
line = fnav.readline()
nav.ion[1, 3] = self.flt(line, 0)
nav.ion_prm[sys][im].prm[7] = self.flt(line, 0)
if len(line) >= 42:
nav.ion_region = int(self.flt(line, 1))
nav.ion_prm[sys][im].prm[8] = \
int(self.flt(line, 1))
continue

elif line[0:5] == '> EPH':
Expand Down Expand Up @@ -540,7 +592,10 @@ def decode_nav(self, navfile, nav, append=False):
eph.sisai[2] = int(self.flt(line, 2)) # oc1
eph.sisai[3] = int(self.flt(line, 3)) # oc2
elif sys == uGNSS.IRN:
eph.urai = int(self.flt(line, 0))
if self.mode_nav == 0:
eph.sva = self.flt(line, 0)
else: # L1NV
eph.urai = self.flt(line, 0)
eph.svh = int(self.flt(line, 1))
if self.mode_nav == 2 and eph.integ == 1:
eph.tgd = int(self.flt(line, 3))
Expand Down

0 comments on commit e8689ee

Please sign in to comment.