From e8689ee93b52be6bbf22fdf32ecb397b2d76d983 Mon Sep 17 00:00:00 2001 From: Rui Hirokawa Date: Sat, 18 Jan 2025 19:44:31 +0900 Subject: [PATCH] - fixed sign of BCNAV1 - added STO/EOP/ION support for RINEX 4.02 --- src/cssrlib/gnss.py | 34 +++++++++- src/cssrlib/rawnav.py | 32 +++++++--- src/cssrlib/rinex.py | 145 +++++++++++++++++++++++++++++------------- 3 files changed, 154 insertions(+), 57 deletions(-) diff --git a/src/cssrlib/gnss.py b/src/cssrlib/gnss.py index 041f155..f968941 100644 --- a/src/cssrlib/gnss.py +++ b/src/cssrlib/gnss.py @@ -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 """ @@ -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 diff --git a/src/cssrlib/rawnav.py b/src/cssrlib/rawnav.py index 07a4488..69d8a5f 100644 --- a/src/cssrlib/rawnav.py +++ b/src/cssrlib/rawnav.py @@ -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 """ @@ -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): @@ -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: @@ -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 @@ -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 @@ -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 diff --git a/src/cssrlib/rinex.py b/src/cssrlib/rinex.py index ca2fb18..ec08a0e 100644 --- a/src/cssrlib/rinex.py +++ b/src/cssrlib/rinex.py @@ -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: @@ -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 """ @@ -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': @@ -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))