From a51f0c555a6559a563fe87411b1f2c959fc886b2 Mon Sep 17 00:00:00 2001 From: Rui Hirokawa Date: Sun, 5 Jan 2025 16:08:12 +0900 Subject: [PATCH] - support RINEX 4.02 - added NavIC L1NAV support - added Glonass L1OC, L3OC support - fixed BDS D1/D2 navigation message decoder --- src/cssrlib/gnss.py | 29 ++- src/cssrlib/rawnav.py | 404 ++++++++++++++++++++++++++++++++++++++++-- src/cssrlib/rinex.py | 403 +++++++++++++++++++++++++++++++---------- 3 files changed, 726 insertions(+), 110 deletions(-) diff --git a/src/cssrlib/gnss.py b/src/cssrlib/gnss.py index 3ee1c28..77b3f83 100644 --- a/src/cssrlib/gnss.py +++ b/src/cssrlib/gnss.py @@ -100,6 +100,7 @@ class rCST(): P2_20 = 9.536743164062500E-07 P2_21 = 4.768371582031250E-07 P2_24 = 5.960464477539063e-08 + P2_26 = 1.490116119384766e-08 P2_27 = 7.450580596923828e-09 P2_28 = 3.725290298461914E-09 P2_29 = 1.862645149230957E-09 @@ -691,15 +692,37 @@ class Geph(): sva = 0 age = 0.0 toe = gtime_t() + toes = 0.0 tof = gtime_t() pos = np.zeros(3) vel = np.zeros(3) acc = np.zeros(3) taun = 0.0 # SV clock bias [s] - gamn = 0.0 # relative frq bias + gamn = 0.0 # SV clock drift [s/s] + beta = 0.0 # SV clock drift rate [s/s^2] dtaun = 0.0 # delta between L1 and L2 [s] mode = 0 - status = 0 + status = 0 # data validity + flag = 0 + + # for CDMA + urai = np.zeros(2, dtype=int) + dpos = np.zeros(3) + + psi = 0.0 # yaw angle [rad] + sn = 0 # sign flag + win = 0.0 # angular rate [rad/s] + dw = 0.0 # angular accel [rad/s^2] + wmax = 0.0 # max angular rate[rad/s] + aode = 0 + aodc = 0 # age of data orbit/clock [days] + tin = 0.0 + tau1 = 0.0 + tau2 = 0.0 + + src = 0 # source flags (b0-1: Rt, b2-3: Re) + sattype = 0 # 0 - M(L3), 1 - K1(L3), 3 - K1(L2/L3), 2 - K2 (L1/L2/L3) + isc = np.zeros(3) # 0: ISC_L1OC, 1: ISC_L2OC, 2: ISC_L3OC def __init__(self, sat=0): self.sat = sat @@ -782,6 +805,8 @@ def __init__(self, nf=2): self.thresar = 3.0 # AR acceptance threshold self.elmaskar = np.deg2rad(20.0) # elevation mask for AR + self.leaps = 18 # leap seconds [s] + # Select tropospheric model # self.trpModel = uTropoModel.SAAST diff --git a/src/cssrlib/rawnav.py b/src/cssrlib/rawnav.py index 2ca0005..c226587 100644 --- a/src/cssrlib/rawnav.py +++ b/src/cssrlib/rawnav.py @@ -8,8 +8,8 @@ [1-3] NAVSTAR GPS Space Segment/Navigation User Segment L1C Interfaces (IS-GPS-800) Rev.J, 2022 [2] Quasi-Zenith Satellite System Interface Specification - Satellite Positioning, Navigation and Timing Service (IS-QZSS-PNT-005), - 2022 + Satellite Positioning, Navigation and Timing Service (IS-QZSS-PNT-006), + 2024 [3] Galileo Open Service Signal-in-Space Interface Control Document (OS SIS ICD) Issue 2.1, 2023 [4] BeiDou Navigation Satellite System Signal In Space Interface Control @@ -17,17 +17,44 @@ [5] GLONASS Interface Control Document (Edition 5.1), 2008 [6] Indian Regional NAvigation Satellite System Signal in Space ICD for Standard Positioning Service (Version 1.1), 2017 +[7] NAVIC Signal in Space ICD for Standard Positioning Service + in L1 Frequency (Version 1.0), 2023 """ import numpy as np import bitstruct.c as bs from cssrlib.gnss import gpst2time, gst2time, bdt2time, rCST -from cssrlib.gnss import prn2sat, uGNSS, sat2prn, bdt2gpst, utc2gpst -from cssrlib.gnss import Eph, Geph, uTYP +from cssrlib.gnss import prn2sat, uGNSS, sat2prn, bdt2gpst, utc2gpst, time2gpst +from cssrlib.gnss import Eph, Geph, uTYP, copy_buff, gtime_t, timeadd from cssrlib.rinex import rnxenc, rSigRnx +def gep2time(N4, Nt, sod): + """ calculate time from Glonass epoch """ + time = gtime_t() + + if Nt <= 366: + j = 1 + doy = Nt + elif Nt <= 731: + j = 2 + doy = Nt-366 + elif Nt <= 1096: + j = 3 + doy = Nt-731 + else: + j = 4 + doy = Nt-1096 + + year = 1996 + 4*(N4-1) + j-1 + days = (year-1970)*365+(year-1969)//4+doy-1 + sec = int(sod) + time.time = days*86400+sec + time.sec = sod-sec + return utc2gpst(timeadd(time, -10800.0)) + + class RawNav(): def __init__(self, opt=None, prefix=''): self.gps_lnav = {} @@ -70,6 +97,12 @@ def __init__(self, opt=None, prefix=''): for k in range(uGNSS.GLOMAX): self.glo_ca[k] = bytearray(200) + self.glo_l1oc = {} + self.glo_l2oc = {} + self.glo_l3oc = {} + + self.monlevel = 1 + self.irn_nav = {} for k in range(uGNSS.IRNMAX): self.irn_nav[k] = bytearray(200) @@ -574,8 +607,7 @@ def decode_gps_cnav(self, week, time, sat, msg): eph.M0 = M0*rCST.P2_32*rCST.SC2RAD eph.e = e*rCST.P2_34 eph.omg = omg*rCST.P2_32*rCST.SC2RAD - eph.integ = isf - # eph.esf = esf + eph.integ = (alert << 2) | (esf << 1) | isf eph.wn_op = wn//256*256 + wn_op @@ -595,7 +627,7 @@ def decode_gps_cnav(self, week, time, sat, msg): eph.mode = 1 # CNAV eph.toc = gpst2time(eph.week, toc) eph.toe = gpst2time(eph.week, eph.toes) - eph.tot = bdt2time(eph.week, tow) + eph.tot = gpst2time(eph.week, tow) eph.top = gpst2time(eph.wn_op, eph.tops) if wn_op < 0: @@ -681,8 +713,10 @@ def decode_gps_cnav2(self, week, time, sat, msg): eph.M0 = M0*rCST.P2_32*rCST.SC2RAD eph.e = e*rCST.P2_34 eph.omg = omg*rCST.P2_32*rCST.SC2RAD - eph.integ = isf - # eph.esf = esf + if sys == uGNSS.GPS: + eph.integ = isf + else: + eph.integ = (esf << 1) | isf eph.wn_op = wn//256*256 + wn_op @@ -702,7 +736,7 @@ def decode_gps_cnav2(self, week, time, sat, msg): eph.mode = 2 # CNAV/2 eph.toe = gpst2time(eph.week, eph.toes) eph.toc = eph.toe - eph.tot = bdt2time(eph.week, tow) + eph.tot = gpst2time(eph.week, tow) eph.top = gpst2time(eph.wn_op, eph.tops) return eph @@ -823,6 +857,10 @@ def decode_bds_b1c(self, week, time_, sat, msg): eph.week, how, eph.iodc, eph.iode = bs.unpack_from( 'u13u8u10u8', msg, i) i += 39 + + if (eph.iodc & 0xff) != eph.iode: + return None + tow = how*3600+soh eph.tot = bdt2time(eph.week, tow) @@ -921,6 +959,10 @@ def decode_bds_b2a(self, week, time_, sat, msg): eph.iodc, tgd_b2ap, isc_b2ad = bs.unpack_from('u10s12s12', buff, i) i += 34 + + if (eph.iodc & 0xff) != eph.iode: + return None + i = self.decode_bds_cnav_iono(buff, i) tgd_b1cp = bs.unpack_from('u12', buff, i)[0] @@ -1047,7 +1089,7 @@ def decode_bds_d1(self, week, time, sat, msg): eph.crs = self.getbits2(buff, i+224, 8, i+240, 10)*rCST.P2_6 sqrtA = self.getbitu2(buff, i+250, 12, i+270, 20)*rCST.P2_19 eph.A = sqrtA**2 - toe1 = bs.unpack_from('s2', buff, i+290)[0] + toe1 = bs.unpack_from('u2', buff, i+290)[0] # subframe 3 i = 320*2 @@ -1252,7 +1294,8 @@ def decode_irn_lnav(self, week, time, sat, msg): eph.af1 = af1*rCST.P2_43 eph.af2 = af2*rCST.P2_55 - eph.urai = self.urai2sva(urai) + eph.urai = urai + eph.sva = self.urai2sva(urai) toc *= 16.0 eph.tgd = tgd*rCST.P2_31 eph.deln = deln*rCST.P2_41*rCST.SC2RAD @@ -1287,6 +1330,177 @@ def decode_irn_lnav(self, week, time, sat, msg): eph.mode = 0 return eph + def decode_irn_l1nav_iono_grid(self, msg, i): + mask, regid = bs.unpack_from('u10u4', msg, i) + i += 14 + for k in range(15): + givei, givd = bs.unpack_from('u10u4', msg, i) + i += 13 + iodi = bs.unpack_from('u3', msg, i)[0] + i += 3+31 + return i + + def decode_irn_l1nav_alm(self, msg, i): + wna, e, toa, i0, OMGd, sqrtA, OMG0, omg, M0, af0, af1, prn_a = \ + bs.unpack_from('u13u20u16s24s19u24s24s24s24s14s11u6', msg, i) + i += 243 + return i + + def decode_irn_l1nav_iono_nequick(self, msg, i): + for k in range(3): + modip_mac, modip_min, mlon_max, mlon_min, a0, a1, a2, idf = \ + bs.unpack_from('s6s6s7s7u11s11s14u1', msg, i) + i += 63 + + iodn = bs.unpack_from('u2', msg, i)[0] + i += 2+52 + return i + + def decode_irn_l1nav_iono_klob(self, msg, i): + alp0, alp1, alp2, alp3 = \ + bs.unpack_from('s8s8s10s12', msg, i) + i += 38 + bet0, bet1, bet2, bet3 = \ + bs.unpack_from('s8s8s11s14', msg, i) + i += 41 + + lon_max, lon_min, lat_max, lat_min, iodk = \ + bs.unpack_from('u6u6s5s5u2', msg, i) + i += 24+2 + return i + + def decode_irn_l1nav_eop(self, msg, i): + teop, pmx, pmxd, pmy, pmyd, dut1, dut1d = \ + bs.unpack_from('u16s21s15s21s15s31s19', msg, i) + i += 138 + return i + + def decode_irn_l1nav_utc(self, msg, i): + iodt, tug, wnug, dtls, wnlsf, dn, dtlsf = \ + bs.unpack_from('u3u8u13s8u13u4s8', msg, i) + i += 57 + + a0utc, a1utc, a2utc = bs.unpack_from('s16s13s7', msg, i) + i += 36 + + flg_utc = bs.unpack_from('u1', msg, i)[0] + i += 1 + if flg_utc: + a0n, a1n, a2n = bs.unpack_from('s16s13s7', msg, i) + i += 36 + + for k in range(3): + gnss, flg_gnss = bs.unpack_from('u3u1', msg, i) + i += 4 + if flg_gnss: + a0g, a1g = bs.unpack_from('s16s13', msg, i) + i += 29 + i += 14 + return i + + def decode_irn_l1nav(self, week, time, sat, msg): + """ NavIC L1NAV Message decoder """ + + sys, prn_ = sat2prn(sat) + eph = Eph(sat) + + # subframe 1 (52 syms) + toi = bs.unpack_from('u9', msg, 0)[0] + + # subframe 2 (600) + i = 52 + wn, itow, alrt = bs.unpack_from('u13u8u1', msg, i) + i += 22 + # data 548bit + svh, iodec, urai, toec, dA, Adot, deln, delnd, M0 = \ + bs.unpack_from('u1u4s5u11s26s26s19s23s33', msg, i) + i += 148 + e, omg, OMG0, OMGd, i0, idot, cis, cic, crs, crc, cus, cuc = \ + bs.unpack_from('u33s33s33s25s33s15s16s16s24s24s21s21', msg, i) + i += 294 + af0, af1, af2, tgd, isc1, isc2, rsf = \ + bs.unpack_from('s29s22s15s12s12s12u1', msg, i) + i += 103+3 + prn = bs.unpack_from('u6', msg, i)[0] + + tow = itow*7200+toi*18 + + # subframe 3 (274) + i = 52+1200 + msgid, ivld = bs.unpack_from('u6u1', msg, i) + i += 7 + + # message data + if msgid == 5: # iono grid parameters + i = self.decode_irn_l1nav_iono_grid(msg, i) + elif msgid == 6: # almanac + i = self.decode_irn_l1nav_alm(msg, i) + elif msgid == 8: # NeQuick-N Iono parameters + i = self.decode_irn_l1nav_iono_nequick(msg, i) + elif msgid == 10: # Klobuchar like Iono coefficient & EOP + i = self.decode_irn_l1nav_eop(msg, i) + i = self.decode_irn_l1nav_iono_klob(msg, i) + elif msgid == 17: # NavIC Time offsets + i = self.decode_irn_l1nav_utc(msg, i) + elif msgid == 0: # null message + None + + i += 243 + + eph.isc = np.zeros(6) + + eph.week = wn + eph.code = 2 + eph.svh = svh + + eph.iodc = iodec + eph.iode = iodec + + # clock + eph.af2 = af2*rCST.P2_66 + eph.af1 = af1*rCST.P2_50 + eph.af0 = af0*rCST.P2_35 + eph.urai = urai + eph.sva = self.urai2sva(urai) + + # group-delay + eph.tgd = tgd*rCST.P2_35 # L1CA + eph.isc[4] = isc2*rCST.P2_35 # L1D + eph.isc[5] = isc1*rCST.P2_35 # L1P if rsf=0, S if rsf=1 + + # ephemeris + # eph.tops = top*300.0 + A0 = 26559710.0 if sys == uGNSS.GPS else 42164200.0 + eph.A = A0 + dA*rCST.P2_9 + eph.Adot = Adot*rCST.P2_21 + eph.deln = deln*rCST.P2_44*rCST.SC2RAD + eph.delnd = delnd*rCST.P2_57*rCST.SC2RAD + eph.M0 = M0*rCST.P2_32*rCST.SC2RAD + eph.e = e*rCST.P2_34 + eph.omg = omg*rCST.P2_32*rCST.SC2RAD + + eph.integ = rsf + + eph.toes = toec*300.0 + eph.OMG0 = OMG0*rCST.P2_32*rCST.SC2RAD + eph.i0 = i0*rCST.P2_32*rCST.SC2RAD + eph.OMGd = OMGd*rCST.P2_44*rCST.SC2RAD + eph.idot = idot*rCST.P2_44*rCST.SC2RAD + + eph.cis = cis*rCST.P2_30 + eph.cic = cic*rCST.P2_30 + eph.crs = crs*rCST.P2_8 + eph.crc = crc*rCST.P2_8 + eph.cus = cus*rCST.P2_30 + eph.cuc = cuc*rCST.P2_30 + + eph.mode = 2 # L1NAV + eph.toe = gst2time(eph.week, eph.toes) + eph.toc = eph.toe + eph.tot = gst2time(eph.week, tow) + + return eph + def decode_glo_fdma(self, week, time, sat, msg, freq=0): """ Glonass FDMA navigation message decoder """ @@ -1396,6 +1610,172 @@ def decode_glo_fdma(self, week, time, sat, msg, freq=0): geph.mode = 0 return geph + def decode_glo_l1oc(self, week, time, sat, msg): + """ Glonass CDMA L1OC navigation message decoder """ + sid_t = {10: 0, 11: 1, 12: 2, 16: 3, 25: 4} + + sys, prn = sat2prn(sat) + pre, sid = bs.unpack_from('u12u6', msg, 0) + if pre != 0x5f1: + return None + if sid not in (10, 11, 12): + return None + + if sat not in self.glo_l1oc.keys(): + self.glo_l1oc[sat] = bytearray(32*5) + + buff = self.glo_l1oc[sat] + i0 = sid_t[sid]*32*8 + copy_buff(msg, buff, 12, i0, 222) + geph = self.decode_glo_cdma(week, time, sat, buff, 1) + return geph + + def decode_glo_l3oc(self, week, time, sat, msg): + """ Glonass CDMA L3OC navigation message decoder """ + sid_t = {10: 0, 11: 1, 12: 2, 16: 3, 25: 4} + + pre, sid = bs.unpack_from('u20u6', msg, 0) + if pre != 0x494e: + return None + + if sid not in sid_t: + return None + + if sat not in self.glo_l3oc.keys(): + self.glo_l3oc[sat] = bytearray(32*5) + + buff = self.glo_l3oc[sat] + i0 = sid_t[sid]*32*8 + copy_buff(msg, buff, 20, i0, 256) + geph = self.decode_glo_cdma(week, time, sat, buff, 3) + return geph + + def decode_glo_cdma(self, week, time, sat, buff, stype): + """ Glonass CDMA L3OC navigation message decoder """ + + sid10 = bs.unpack_from('u6', buff, 0)[0] + sid11 = bs.unpack_from('u6', buff, 32*8)[0] + sid12 = bs.unpack_from('u6', buff, 32*2*8)[0] + # sid16 = bs.unpack_from('u6', buff, 32*3*8)[0] + + # note: + # currently, only MT10,11,12 are broadcast in L1OC, L3OC + # it need to be fixed to include MT16 once it is available. + if sid10 != 10 or sid11 != 11 or sid12 != 12: + return None + + geph = Geph(sat) + + # MT10 + i = 6 + # service field + if stype == 1: # L1OS + svid, svh, valid, P1, P2, KP, A, ts = \ + bs.unpack_from('u6u1u1u4u1u2u1u16', buff, i) + i += 32 + else: # L3OS + ts, svid, svh, valid, P1, P2, KP, A = \ + bs.unpack_from('u15u6u1u1u4u1u2u1', buff, i) + i += 31 + + N4, Nt, M, PS, tb, Ee, Et, Re, Rt, Fe, Ft = \ + bs.unpack_from('u5u11u3u6u10u8u8u2u2s5s5', buff, i) + i += 65 + + # N4: four-year interval + # Nt: number of day + + # PS: pseudoframe size + + # Ee, Et: age of ephemeris, clock + # Re, Rt: regime for generation of ephemeris/clock + # Fe, Ft: ephemeris, clock accuracy index urai orb/clk + + tau, gam, bet, tau_c, tc = \ + bs.unpack_from('u32u19u15u40u13', buff, i) + + # MT11 + i = 32*8+6 + + ts, svid, svh, valid, P1, P2, KP, A = \ + bs.unpack_from('u15u6u1u1u4u1u2u1', buff, i) + + i += 31 + x, y, z, vx, vy = \ + bs.unpack_from('s40s40s40s35s35', buff, i) + + # MT12 + i = 32*2*8+6 + + ts, svid, svh, valid, P1, P2, KP, A = \ + bs.unpack_from('u15u6u1u1u4u1u2u1', buff, i) + i += 31 + vz, ax, ay, az, dx, dy, dz, dtau_L3, tau_gps = \ + bs.unpack_from('s35s15s15s15s13s13s13s18s30', buff, i) + + # MT16 + i = 32*3*8+6 + + # ts, svid, svh, valid, P1, P2, KP, A = \ + # bs.unpack_from('u15u6u1u1u4u1u2u1', buff, i) + i += 31 + tin, psi, sn, wmax, win, dw, tau1, tau2 = \ + bs.unpack_from('u22u15u1u17u17u15u13u17', buff, i) + + geph.psi = psi*rCST.P2_14*rCST.SC2RAD + geph.win = win*rCST.P2_26*rCST.SC2RAD + geph.wmax = wmax*rCST.P2_26*rCST.SC2RAD + geph.dw = dw*rCST.P2_30*rCST.SC2RAD + geph.sn = sn # sign flag of maneuver + + geph.tin = tin*rCST.P2_5 + geph.tau1 = tau1*rCST.P2_5 + geph.tau2 = tau2*rCST.P2_5 + + geph.pos[0] = x*rCST.P2_20*1e3 + geph.pos[1] = y*rCST.P2_20*1e3 + geph.pos[2] = z*rCST.P2_20*1e3 + geph.vel[0] = vx*rCST.P2_30*1e3 + geph.vel[1] = vy*rCST.P2_30*1e3 + geph.vel[2] = vz*rCST.P2_30*1e3 + geph.acc[0] = ax*rCST.P2_39*1e3 + geph.acc[1] = ay*rCST.P2_39*1e3 + geph.acc[2] = az*rCST.P2_39*1e3 + + geph.dpos[0] = dx*rCST.P2_10 + geph.dpos[1] = dy*rCST.P2_10 + geph.dpos[2] = dz*rCST.P2_10 + + # KP: expected UTC(SU) correection + # A: L1OCd time correction in next string 0: no plan + + geph.taun = tau*rCST.P2_38 + geph.gamn = gam*rCST.P2_48 + geph.beta = bet*rCST.P2_57 + + geph.sattype = M # 0:M(L3),1:K1(L3),3:K1(L2/L3),2:K2(L1/L2/L3) + geph.svh = svh # 0:health, 1:unhealthy + geph.status = valid # data validity 0:valid, 1:invalid + geph.flag = P2 # 0: yaw steering, 1: rate-limited yaw manoeuvre + geph.src = (Re << 2) | Rt + + geph.aode = Ee*0.25 # age of ephemeris [days] + geph.aodc = Et*0.25 # age of clock [days] + + geph.urai[0] = Fe + geph.urai[1] = Ft + + geph.isc[2] = dtau_L3*rCST.P2_38 + + geph.toe = gep2time(N4, Nt, tb*90.0) + _, geph.toes = time2gpst(geph.toe) + geph.tof = gep2time(N4, Nt, ts*3.0) + geph.iode = tb + + geph.mode = stype # FDMA:0,L1OC:1,L2OC:2,L3OC:3 + + return geph + class rcvOpt(): flg_qzslnav = False diff --git a/src/cssrlib/rinex.py b/src/cssrlib/rinex.py index 0e274c1..a6cdbfe 100644 --- a/src/cssrlib/rinex.py +++ b/src/cssrlib/rinex.py @@ -1,5 +1,8 @@ """ module for RINEX 3.0x processing + +[1] RINEX: The Receiver Independent Exchange Format Version 4.02, 2024 + """ import numpy as np @@ -8,7 +11,7 @@ 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 +from cssrlib.gnss import timeadd, id2sat, gpst2utc, Seph class pclk_t: @@ -161,6 +164,8 @@ def decode_nav(self, navfile, nav, append=False): if line[0:4] == 'GPSB' or line[0:4] == 'QZSB': for k in range(4): nav.ion[1, k] = self.flt(line[5+k*12:5+(k+1)*12]) + elif line[60:72] == 'LEAP SECONDS': + nav.leaps = int(line[:6]) for line in fnav: @@ -205,12 +210,14 @@ def decode_nav(self, navfile, nav, append=False): elif line[0:5] == '> ION': # iono (TBD) sys = char2sys(line[6]) itype = line[10:14] + stype = '' if len(line) < 20 else line[15:19] line = fnav.readline() ttm = self.decode_time(line, 4) if sys == uGNSS.GAL and itype == 'IFNV': # Nequick-G - for k in range(3): + for k in range(3): # ai0, ai1, ai2 nav.ion[0, k] = self.flt(line, k+1) line = fnav.readline() + # disturbance flags nav.ion[0, 3] = int(self.flt(line, 0)) elif sys == uGNSS.BDS and itype == 'CNVX': # BDGIM ttm = self.decode_time(line, 4) @@ -223,6 +230,41 @@ def decode_nav(self, navfile, nav, append=False): line = fnav.readline() for k in range(2): nav.ion_gim[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) + line = fnav.readline() + for k in range(4): + nav.ion[1, k] = self.flt(line, k) + line = fnav.readline() + nav.ion_region = np.zeros(4) + for k in range(4): + nav.ion_region[k] = self.flt(line, k) + + elif stype == 'NEQN': + nav.ion_region = np.zeros((3, 4)) + iodn = self.flt(line, 1) + + 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) + 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) + + 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] + else: # Klobucher (LNAV, D1D2, CNVX) self.ion_gim = np.zeros(9) for k in range(3): @@ -243,12 +285,16 @@ def decode_nav(self, navfile, nav, append=False): m = line[10:14] if m == 'CNAV' or m == 'CNV1' or m == 'FNAV': self.mode_nav = 1 - elif m == 'CNV2': + elif m == 'CNV2' or m == 'L1NV': self.mode_nav = 2 elif m == 'CNV3': self.mode_nav = 3 elif m == 'FDMA': self.mode_nav = 0 + elif m == 'L1OC': + self.mode_nav = 1 + elif m == 'L3OC': + self.mode_nav = 3 elif m == 'SBAS': self.mode_nav = 0 line = fnav.readline() @@ -265,9 +311,6 @@ def decode_nav(self, navfile, nav, append=False): if sys == uGNSS.GLO: prn = int(line[1:3]) sat = prn2sat(sys, prn) - pos = np.zeros(3) - vel = np.zeros(3) - acc = np.zeros(3) geph = Geph(sat) geph.mode = self.mode_nav @@ -279,42 +322,41 @@ def decode_nav(self, navfile, nav, append=False): geph.taun = -self.flt(line, 1) geph.gamn = self.flt(line, 2) - t0 = self.flt(line, 3) - - tod = t0 % 86400.0 - tof = gpst2time(week, tod + dow*86400.0) - tof = self.adjday(tof, toc) - - geph.toe = utc2gpst(toc) - geph.tof = utc2gpst(tof) - - # iode = Tb(7bit) - geph.iode = int(((tocs+10800.0) % 86400)/900.0+0.5) + if self.mode_nav == 0: # FDMA + t0 = self.flt(line, 3) + else: # L1OC, L3OC + bet_ = self.flt(line, 3) # clock drift rate line = fnav.readline() # line #1 - pos[0] = self.flt(line, 0)*1e3 - vel[0] = self.flt(line, 1)*1e3 - acc[0] = self.flt(line, 2)*1e3 + geph.pos[0] = self.flt(line, 0)*1e3 + geph.vel[0] = self.flt(line, 1)*1e3 + geph.acc[0] = self.flt(line, 2)*1e3 geph.svh = int(self.flt(line, 3)) line = fnav.readline() # line #2 - pos[1] = self.flt(line, 0)*1e3 - vel[1] = self.flt(line, 1)*1e3 - acc[1] = self.flt(line, 2)*1e3 - geph.frq = int(self.flt(line, 3)) + geph.pos[1] = self.flt(line, 0)*1e3 + geph.vel[1] = self.flt(line, 1)*1e3 + geph.acc[1] = self.flt(line, 2)*1e3 - if geph.frq > 128: - geph.frq -= 256 + if self.mode_nav == 0: # FDMA + geph.frq = int(self.flt(line, 3)) + + if geph.frq > 128: + geph.frq -= 256 + else: # L1OC + dvalid = int(self.flt(line, 3)) line = fnav.readline() # line #3 - pos[2] = self.flt(line, 0)*1e3 - vel[2] = self.flt(line, 1)*1e3 - acc[2] = self.flt(line, 2)*1e3 - geph.age = int(self.flt(line, 3)) + geph.pos[2] = self.flt(line, 0)*1e3 + geph.vel[2] = self.flt(line, 1)*1e3 + geph.acc[2] = self.flt(line, 2)*1e3 - geph.pos = pos - geph.vel = vel - geph.acc = acc + if self.mode_nav == 0: # FDMA + geph.age = int(self.flt(line, 3)) + elif self.mode_nav == 1: # L1OC + tgd_L2OCp = self.flt(line, 3) # tgd_L2OCp + elif self.mode_nav == 3: # L3OC + isc_L3OCp = self.flt(line, 3) # isc_L3OCp # Use GLONASS line #4 only from RINEX v3.05 onwards # @@ -322,16 +364,88 @@ def decode_nav(self, navfile, nav, append=False): line = fnav.readline() # line #4 - # b7-8: M, b6: P4, b5: P3, b4: P2, b2-3: P1, b0-1: P - geph.status = int(self.flt(line, 0)) - geph.dtaun = -self.flt(line, 1) - geph.urai = int(self.flt(line, 2)) - # svh = int(self.flt(line, 3)) + if self.mode_nav == 0: # FDMA + # b7-8: M, b6: P4, b5: P3, b4: P2, b2-3: P1, b0-1: P + geph.status = int(self.flt(line, 0)) + geph.dtaun = -self.flt(line, 1) + geph.urai = int(self.flt(line, 2)) + if len(line) >= 80: + geph.svh = int(self.flt(line, 3)) + else: # L1OC,L3OC + sattype = int(self.flt(line, 0)) + src = int(self.flt(line, 1)) + aode_ee = int(self.flt(line, 2)) + aode_et = int(self.flt(line, 3)) + + line = fnav.readline() # line #5 + P2 = int(self.flt(line, 0)) # attitude flag + t0 = self.flt(line, 1) # sec of day, UTC(SU) + tau1 = self.flt(line, 2) + tau2 = self.flt(line, 3) + + line = fnav.readline() # line #6 + yaw = self.flt(line, 0) + sgn = int(self.flt(line, 1)) + win = self.flt(line, 2) + dw = self.flt(line, 3) + + line = fnav.readline() # line #7 + wmax = self.flt(line, 0) + dxpc = self.flt(line, 1) + dypc = self.flt(line, 2) + dzpc = self.flt(line, 3) + + line = fnav.readline() # line #8 + urai_orb = int(self.flt(line, 0)) + urai_clk = int(self.flt(line, 1)) + tot = self.flt(line, 2) + + tod = t0 % 86400.0 + tof = gpst2time(week, tod + dow*86400.0) + tof = self.adjday(tof, toc) + + geph.toe = utc2gpst(toc) + geph.tof = utc2gpst(tof) + + # iode = Tb(7bit) + geph.iode = int(((tocs+10800.0) % 86400)/900.0+0.5) nav.geph.append(geph) continue - elif sys not in (uGNSS.GPS, uGNSS.GAL, uGNSS.QZS, uGNSS.BDS): + elif sys == uGNSS.SBS: + prn = int(line[1:3])+100 + sat = prn2sat(sys, prn) + seph = Seph(sat) + + seph.toc = self.decode_time(line, 4) + seph.af0 = self.flt(line, 1) + seph.af1 = self.flt(line, 2) + seph.tot = self.flt(line, 3) + + line = fnav.readline() # line #1 + seph.pos[0] = self.flt(line, 0)*1e3 + seph.vel[0] = self.flt(line, 1)*1e3 + seph.pos[0] = self.flt(line, 2)*1e3 + seph.svh = int(self.flt(line, 3)) + + line = fnav.readline() # line #2 + seph.pos[1] = self.flt(line, 0)*1e3 + seph.vel[1] = self.flt(line, 1)*1e3 + seph.pos[1] = self.flt(line, 2)*1e3 + seph.sva = self.flt(line, 3) + + line = fnav.readline() # line #3 + seph.pos[2] = self.flt(line, 0)*1e3 + seph.vel[2] = self.flt(line, 1)*1e3 + seph.pos[2] = self.flt(line, 2)*1e3 + seph.iodn = int(self.flt(line, 3)) + + nav.seph.append(seph) + continue + + elif sys not in (uGNSS.GPS, uGNSS.GAL, uGNSS.QZS, uGNSS.BDS, + uGNSS.IRN): continue prn = int(line[1:3]) @@ -353,7 +467,8 @@ def decode_nav(self, navfile, nav, append=False): line = fnav.readline() # line #1 - if sys == uGNSS.GAL: + if sys == uGNSS.GAL or \ + (sys == uGNSS.IRN and self.mode_nav == 0): eph.iode = int(self.flt(line, 0)) eph.iodc = eph.iode else: @@ -374,7 +489,15 @@ def decode_nav(self, navfile, nav, append=False): eph.A = sqrtA**2 line = fnav.readline() # line #3 - eph.toes = int(self.flt(line, 0)) + if sys == uGNSS.IRN and self.mode_nav == 2: + eph.iode = int(self.flt(line, 0)) + eph.iode = eph.iodc + else: + if (sys == uGNSS.GPS or sys == uGNSS.QZS) and \ + self.mode_nav > 0: # CNAV, CNAV/2 + eph.tops = self.flt(line, 0) + else: + eph.toes = self.flt(line, 0) eph.cic = self.flt(line, 1) eph.OMG0 = self.flt(line, 2) eph.cis = self.flt(line, 3) @@ -395,12 +518,18 @@ def decode_nav(self, navfile, nav, append=False): if sys == uGNSS.GAL and self.ver < 4.0: eph.mode = 1 if eph.code & 0x2 else 0 + elif sys == uGNSS.IRN and self.mode_nav == 0: + eph.week = int(self.flt(line, 2)) + else: eph.delnd = self.flt(line, 1) if sys == uGNSS.BDS: eph.sattype = int(self.flt(line, 2)) eph.tops = int(self.flt(line, 3)) - else: + elif sys == uGNSS.IRN and self.mode_nav == 2: + eph.integ = int(self.flt(line, 3)) # rsf + else: # CNAV, CNAV/2 + eph.urai = [0, 0, 0, 0] eph.urai[0] = int(self.flt(line, 2)) eph.urai[1] = int(self.flt(line, 3)) @@ -410,6 +539,13 @@ def decode_nav(self, navfile, nav, append=False): eph.sisai[1] = int(self.flt(line, 1)) # ocb 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)) + eph.svh = int(self.flt(line, 1)) + if self.mode_nav == 2 and eph.integ == 1: + eph.tgd = int(self.flt(line, 3)) + else: + eph.tgd = int(self.flt(line, 2)) else: eph.sva = int(self.flt(line, 0)) eph.svh = int(self.flt(line, 1)) @@ -418,7 +554,7 @@ def decode_nav(self, navfile, nav, append=False): if self.mode_nav == 0: eph.iodc = int(self.flt(line, 3)) else: - eph.urai[2] = int(self.flt(line, 3)) + eph.urai[2] = int(self.flt(line, 3)) # URAI_NED2 eph.urai[3] = eph.sva # URAI_ED elif sys == uGNSS.GAL: tgd_b = float(self.flt(line, 3)) @@ -437,53 +573,72 @@ def decode_nav(self, navfile, nav, append=False): if self.mode_nav < 3: line = fnav.readline() # line #7 if sys == uGNSS.BDS: - if self.mode_nav == 0: + if self.mode_nav == 0: # D1/D2 tot = self.flt(line, 0) eph.iodc = int(self.flt(line, 1)) - else: - if self.mode_nav == 1: + else: # CNAV-1,2,3 + if self.mode_nav == 1: # CNAV-1 eph.isc[0] = float(self.flt(line, 0)) # B1Cd - elif self.mode_nav == 2: + elif self.mode_nav == 2: # CNAV-2 eph.isc[1] = float(self.flt(line, 1)) # B2ad eph.tgd = float(self.flt(line, 2)) # tgd_B1Cp eph.tgd_b = float(self.flt(line, 3)) # tgd_B2ap + elif sys == uGNSS.IRN: + if self.mode_nav > 0: + if eph.integ == 0: # rsf + eph.isc[5] = float(self.flt(line, 0)) # S + eph.isc[4] = float(self.flt(line, 1)) # L1D + else: + eph.isc[5] = float(self.flt(line, 2)) # L1P + eph.isc[4] = float(self.flt(line, 3)) # L1D + + line = fnav.readline() # line #8 + + tot = self.flt(line, 0) + elif sys == uGNSS.GAL: - tot = int(self.flt(line, 0)) + tot = self.flt(line, 0) + + elif sys in (uGNSS.GPS, uGNSS.QZS): + if self.mode_nav > 0: # CNAV, CNAV/2 + eph.isc[0] = self.flt(line, 0) # ISC_L1CA + eph.isc[1] = self.flt(line, 1) # ISC_L2C + eph.isc[2] = self.flt(line, 2) # ISC_L5I5 + eph.isc[3] = self.flt(line, 3) # ISC_L5Q5 + else: # LNAV + tot = self.flt(line, 0) + if line(line) >= 42: + eph.fit = int(self.flt(line, 1)) - else: - if self.mode_nav > 0 and sys != uGNSS.GAL: - eph.isc[0] = self.flt(line, 0) - eph.isc[1] = self.flt(line, 1) - eph.isc[2] = self.flt(line, 2) - eph.isc[3] = self.flt(line, 3) - line = fnav.readline() + if sys in (uGNSS.GPS, uGNSS.QZS): + if self.mode_nav > 0: # CNAV, CNAV/2 + line = fnav.readline() # line #8 + if self.mode_nav == 2: # CNAV/2 + eph.isc[4] = self.flt(line, 0) # ISC_L1Cd + eph.isc[5] = self.flt(line, 1) # ISC_L1Cp - if self.mode_nav == 2: - eph.isc[4] = self.flt(line, 0) - eph.isc[5] = self.flt(line, 1) - line = fnav.readline() + line = fnav.readline() # line #9 tot = int(self.flt(line, 0)) - if self.mode_nav > 0: - eph.week = int(self.flt(line, 1)) - elif len(line) >= 42: - eph.fit = int(self.flt(line, 1)) + eph.wn_op = int(self.flt(line, 1)) + if len(line) >= 61: # optional + eph.integ = int(self.flt(line, 2)) - if sys == uGNSS.BDS and self.mode_nav > 0: + elif sys == uGNSS.BDS and self.mode_nav > 0: # CNAV-1,2,3 line = fnav.readline() # line #8 eph.sismai = int(self.flt(line, 0)) eph.svh = int(self.flt(line, 1)) eph.integ = int(self.flt(line, 2)) - if self.mode_nav < 3: + if self.mode_nav < 3: # CNAV-1,2 eph.iodc = int(self.flt(line, 3)) - else: - eph.tgd_b = float(self.flt(line, 3)) + else: # CNAV-3 + eph.tgd_b = float(self.flt(line, 3)) # tgd_B2bI line = fnav.readline() # line #9 - tot = int(self.flt(line, 0)) - if self.mode_nav < 3: + tot = self.flt(line, 0) + if self.mode_nav < 3: # CNAV-1,2 eph.iode = int(self.flt(line, 3)) if sys == uGNSS.BDS: @@ -759,7 +914,7 @@ def __init__(self, sig_tab=None): self.rec_eph = {} - def rnx_nav_header(self, fh=None, ver=4.00): + def rnx_nav_header(self, fh=None, ver=4.02): tutc = timeget() tgps = utc2gpst(tutc) leaps = timediff(tgps, tutc) @@ -777,7 +932,7 @@ def rnx_nav_header(self, fh=None, ver=4.00): fh.write("{:60s}{:20s}\n". format("", "END OF HEADER")) - def rnx_obs_header(self, ts: gtime_t, fh=None, ver=4.00): + def rnx_obs_header(self, ts: gtime_t, fh=None, ver=4.02): if self.rnx_obs_header_sent: return @@ -905,19 +1060,16 @@ def rnx_obs_body(self, obs=None, fh=None): def rnx_nav_body(self, eph=None, fh=None): if eph.sat in self.rec_eph.keys(): - if eph.iode in self.rec_eph[eph.sat].keys() and \ - self.rec_eph[eph.sat][eph.iode][0] == eph.mode: + if eph.mode in self.rec_eph[eph.sat].keys() and \ + self.rec_eph[eph.sat][eph.mode][0] == eph.iode: return else: self.rec_eph[eph.sat] = {} - self.rec_eph[eph.sat][eph.iode] = [eph.mode, eph.toes] + self.rec_eph[eph.sat][eph.mode] = [eph.iode, eph.toes] id_ = sat2id(eph.sat) sys, prn = sat2prn(eph.sat) - if sys == uGNSS.BDS and (eph.iodc & 0xff) != eph.iode: - return - if sys == uGNSS.BDS: ep = time2epoch(gpst2bdt(eph.toc)) week, tot_ = time2bdt(eph.tot) @@ -956,9 +1108,14 @@ def rnx_nav_body(self, eph=None, fh=None): v1 = float(eph.iode) v2 = eph.toes elif sys == uGNSS.IRN: - lbl = "LNAV" - v1 = float(eph.iode) - v2 = eph.toes + if eph.mode == 0: + lbl = "LNAV" + v1 = float(eph.iode) + v2 = eph.toes + else: + lbl = "L1NV" + v1 = eph.Adot + v2 = eph.iode else: return @@ -1037,8 +1194,8 @@ def rnx_nav_body(self, eph=None, fh=None): fh.write(" {:19.12E}{:19.12E}{:19s}{:19s}\n". format(float(eph.isc[4]), float(eph.isc[5]), "", "")) - fh.write(" {:19.12E}{:19.12E}{:19s}{:19s}\n". - format(tot_, float(eph.wn_op), "", "")) + fh.write(" {:19.12E}{:19.12E}{:19.12E}{:19s}\n". + format(tot_, float(eph.wn_op), float(eph.integ), "")) if sys == uGNSS.GAL: # I/NAV, F/NAV fh.write(" {:19.12E}{:19.12E}{:19.12E}{:19s}\n". @@ -1059,15 +1216,36 @@ def rnx_nav_body(self, eph=None, fh=None): "")) fh.write(" {:19.12E}{:19s}{:19s}{:19s}\n". format(tot_, "", "", "")) + elif eph.mode == 2: # L1NV + rsf = eph.integ + fh.write(" {:19.12E}{:19.12E}{:19s}{:19.12E}\n". + format(eph.idot, eph.delnd, "", rsf)) + if rsf == 0: + fh.write(" {:19.12E}{:19.12E}{:19.12E}{:19s}\n". + format(float(eph.urai), float(eph.svh), eph.tgd, + "")) + fh.write(" {:19.12E}{:19.12E}{:19s}{:19s}\n". + format(float(eph.isc[5]), float(eph.isc[4]), "", + "")) + else: # rsf = 1 + fh.write(" {:19.12E}{:19.12E}{:19s}{:19.12E}\n". + format(float(eph.urai), float(eph.svh), "", + eph.tgd)) + fh.write(" {:19s}{:19s}{:19.12E}{:19.12E}\n". + format("", "", float(eph.isc[5]), + float(eph.isc[4]))) + + fh.write(" {:19.12E}{:19s}{:19s}{:19s}\n". + format(tot_, "", "", "")) def rnx_gnav_body(self, geph=None, fh=None): if geph.sat in self.rec_eph.keys(): - if geph.iode in self.rec_eph[geph.sat].keys() and \ - self.rec_eph[geph.sat][geph.iode][0] == geph.mode: + if geph.mode in self.rec_eph[geph.sat].keys() and \ + self.rec_eph[geph.sat][geph.mode][0] == geph.iode: return else: self.rec_eph[geph.sat] = {} - self.rec_eph[geph.sat][geph.iode] = [geph.mode, geph.toes] + self.rec_eph[geph.sat][geph.mode] = [geph.iode, geph.toes] id_ = sat2id(geph.sat) sys, prn = sat2prn(geph.sat) @@ -1075,25 +1253,58 @@ def rnx_gnav_body(self, geph=None, fh=None): if sys != uGNSS.GLO: return - lbl = "FDMA" ep = time2epoch(gpst2utc(geph.toe)) week, tot_ = time2gpst(geph.tof) + if geph.mode == 0: + lbl = "FDMA" + v1 = tot_ + v2 = float(geph.frq) + v3 = float(geph.age) + elif geph.mode == 1: + lbl = "L1OC" + v1 = geph.beta + v2 = float(geph.status) + v3 = geph.isc[1] # tgd_L2OCp + else: + lbl = "L3OC" + v1 = geph.beta + v2 = float(geph.status) + v3 = geph.isc[2] # ISC_L3OC + fh.write("> {:2s} {:3s} {:2s}\n".format("EPH", id_, lbl)) fh.write("{:3s} {:4d} {:02d} {:02d} {:02d} {:02d} {:02d}". format(id_, int(ep[0]), int(ep[1]), int(ep[2]), int(ep[3]), int(ep[4]), int(ep[5]))) fh.write("{:19.12E}{:19.12E}{:19.12E}\n". - format(-geph.taun, geph.gamn, tot_)) + format(-geph.taun, geph.gamn, v1)) fh.write(" {:19.12E}{:19.12E}{:19.12E}{:19.12E}\n". - format(geph.pos[0]*1e-3, geph.vel[0]*1e-3, geph.acc[0]*1e-3, - float(geph.svh))) + format(geph.pos[0]*1e-3, geph.vel[0]*1e-3, + geph.acc[0]*1e-3, float(geph.svh))) fh.write(" {:19.12E}{:19.12E}{:19.12E}{:19.12E}\n". - format(geph.pos[1]*1e-3, geph.vel[1]*1e-3, geph.acc[1]*1e-3, - float(geph.frq))) + format(geph.pos[1]*1e-3, geph.vel[1]*1e-3, + geph.acc[1]*1e-3, v2)) fh.write(" {:19.12E}{:19.12E}{:19.12E}{:19.12E}\n". - format(geph.pos[2]*1e-3, geph.vel[2]*1e-3, geph.acc[2]*1e-3, - float(geph.age))) - fh.write(" {:19.12E}{:19.12E}{:19.12E}{:19s}\n". - format(float(geph.flag), float(geph.dtaun), float(geph.sva), - "")) + format(geph.pos[2]*1e-3, geph.vel[2]*1e-3, + geph.acc[2]*1e-3, v3)) + + if geph.mode == 0: # FDMA + fh.write(" {:19.12E}{:19.12E}{:19.12E}{:19s}\n". + format(float(geph.flag), float(geph.dtaun), + float(geph.sva), "")) + else: # L1OC, L3OC + fh.write(" {:19.12E}{:19.12E}{:19.12E}{:19.12E}\n". + format(float(geph.sattype), float(geph.src), + geph.aode, geph.aodc)) + fh.write(" {:19.12E}{:19.12E}{:19.12E}{:19.12E}\n". + format(float(geph.flag), geph.tin, + geph.tau1, geph.tau2)) + fh.write(" {:19.12E}{:19.12E}{:19.12E}{:19.12E}\n". + format(geph.psi, float(geph.sn), + geph.win, geph.dw)) + fh.write(" {:19.12E}{:19.12E}{:19.12E}{:19.12E}\n". + format(geph.wmax, float(geph.dpos[0]), + geph.dpos[1], geph.dpos[2])) + fh.write(" {:19.12E}{:19.12E}{:19s}{:19.12E}\n". + format(float(geph.urai[0]), float(geph.urai[1]), "", + tot_))