- added mode to use L1C/B instead of L1C/A

- scan mode for RTCM decoder

Author: hirokawa

receiver/decode_jps.py

@@ -1011,6 +1011,9 @@ def main():
     parser.add_argument("-j", "--jobs", default=int(mp.cpu_count() / 2),
                         type=int, help='Max. number of parallel processes')
 
+    parser.add_argument("--useL1CB", action='store_true',
+                        help="use L1C/B as like L1C/A for QZS")
+
     # Retrieve all command line arguments
     #
     args = parser.parse_args()
@@ -1033,6 +1036,8 @@ def main():
 
     opt.flg_gpslnav = True
 
+    opt.useL1CB = args.useL1CB
+
     # Start processing pool
     #
     with mp.Pool(processes=args.jobs) as pool:

receiver/decode_lgr.py

@@ -8,6 +8,8 @@
 @author Rui Hirokawa
 """
 
+from skyfield.framelib import itrs
+from skyfield.api import load
 import argparse
 from glob import glob
 import multiprocessing as mp
@@ -16,7 +18,8 @@
 from pathlib import Path
 import struct as st
 from crccheck.crc import Crc24LteA
-from cssrlib.gnss import uGNSS, uTYP, prn2sat, Obs, rSigRnx, gpst2time, uSIG, rCST, gtime_t
+from cssrlib.gnss import uGNSS, uTYP, prn2sat, Obs, rSigRnx, gpst2time, \
+    uSIG, rCST, gtime_t, time2gpst
 from cssrlib.rawnav import rcvDec, rcvOpt
 
 
@@ -83,6 +86,11 @@ def __init__(self, opt=None, prefix='', gnss_t='GE'):
 
         self.nav = navRec()
 
+        self.fn = open('lgr-nav.txt', 'w')
+
+        self.fn.write("# week, tow, nsat, x, y, z, vx, vy, vz, cb, cd, " +
+                      "ggto, sigp, sigv, sigt, pdop, hdop, vdop\n")
+
     def sync(self, buff, k):
         return buff[k] == 0x71  # 'q'
 
@@ -170,6 +178,18 @@ def decode_nav(self, buff, k=10):
         self.nav.std = np.array([stdp, stdv, stdt])
         self.nav.dops = np.array([gdop, pdop, hdop, vdop, tdop])
 
+        return self.nav
+
+    def output_nav(self, nav):
+        f = self.fn
+        week, tow = time2gpst(nav.t)
+        f.write(f"{week:4d} {tow:6.1f} {nav.nsat:2d} ")
+        f.write(f"{nav.pos[0]:11.1f} {nav.pos[1]:11.1f} {nav.pos[2]:11.1f} ")
+        f.write(f"{nav.vel[0]:8.1f} {nav.vel[1]:8.1f} {nav.vel[2]:8.1f} ")
+        f.write(f"{nav.clkb:10.1f} {nav.clkd:8.1f} {nav.ggto:4.1f} ")
+        f.write(f"{nav.std[0]:8.1f} {nav.std[1]:8.1f} {nav.std[2]:8.1f} ")
+        f.write(f"{nav.dops[1]:8.1f} {nav.dops[2]:8.1f} {nav.dops[3]:8.1f}\n")
+
     def decode_obs(self, buff, k=10):
         """ decode RAW message """
         obs = Obs()
@@ -281,7 +301,8 @@ def decode(self, buff, len_, sys=[], prn=[]):
                 self.re.rnx_obs_header(obs.time, self.fh_rnxobs)
                 self.re.rnx_obs_body(obs, self.fh_rnxobs)
         elif mt == b'NAV':
-            self.decode_nav(buff)
+            nav = self.decode_nav(buff)
+            self.output_nav(nav)
         elif mt == b'IQS':
             self.decode_iqs(buff)
         return 0
@@ -333,6 +354,61 @@ def decode(f, opt, args):
     lgrdec.file_close()
 
 
+def ecef_to_lunar_fixed(x_m, y_m, z_m, target_time=None):
+    """
+    Converts ECEF (ITRS) coordinates to the Lunar-fixed coordinate system (Selenocentric).
+    """
+    # 1. Load ephemeris and timescale data
+    # DE421 is a standard JPL ephemeris; DE440 is a more recent alternative.
+    eph = load('de440.bsp')
+    earth = eph['earth']
+    moon = eph['moon']
+    ts = load.timescale()
+
+    if target_time is None:
+        t = ts.now()
+    else:
+        t = target_time
+
+    # 2. Define ECEF coordinates as Earth-fixed (ITRS)
+    # This creates a position object relative to the Earth's center in the ITRS frame.
+    ecef_pos = itrs.at(t, x_m=x_m, y_m=y_m, z_m=z_m)
+
+    # 3. Get the point's position in the Inertial Coordinate System (ICRF)
+    # Calling .at(t) on an ITRS position transforms it into the inertial ICRF frame.
+    pos_icrf = ecef_pos
+
+    # 4. Get the Moon's center position in the ICRF frame
+    moon_icrf = moon.at(t)
+
+    # 5. Calculate the relative vector from the Moon's center to the point
+    # Relative Vector = Point Position (ICRF) - Moon Position (ICRF)
+    relative_vec_icrf = pos_icrf - moon_icrf
+
+    # 6. Transform to the Lunar-fixed frame (considering rotation and libration)
+    # Using the "Mean Earth/Polar Axis" (ME) frame via the moon_pa_de421 kernel.
+    # Note: .frame_xyz() rotates the inertial vector into the specified body-fixed frame.
+    lunar_fixed_pos = relative_vec_icrf.frame_xyz(load('moon_pa_de421.tf'))
+
+    # Retrieve the coordinates in meters
+    x_l, y_l, z_l = lunar_fixed_pos.m
+
+    return x_l, y_l, z_l
+
+
+# --- Example Usage ---
+# Example: ECEF coordinates for a point on Earth (e.g., Tokyo area)
+x_ecef = -3957224.0
+y_ecef = 3310210.0
+z_ecef = 3737512.0
+
+lx, ly, lz = ecef_to_lunar_fixed(x_ecef, y_ecef, z_ecef)
+
+print(f"Time (UTC): {load.timescale().now().utc_strftime()}")
+print(f"Input ECEF (m): X={x_ecef}, Y={y_ecef}, Z={z_ecef}")
+print(f"Output Lunar Fixed (m): X={lx:.2f}, Y={ly:.2f}, Z={lz:.2f}")
+
+
 def main():
 
     # Parse command line arguments

receiver/decode_rtcm.py

@@ -14,7 +14,7 @@
 import os
 from pathlib import Path
 
-from cssrlib.gnss import uGNSS, uTYP, rSigRnx, Obs, gtime_t, timediff
+from cssrlib.gnss import uGNSS, uTYP, rSigRnx, Obs, gtime_t, timediff, sat2prn
 from cssrlib.rawnav import rcvDec, rcvOpt
 
 from cssrlib.rtcm import rtcm
@@ -41,11 +41,11 @@ def __init__(self, opt=None, prefix='', gnss_t='GECJ'):
 
         if 'G' in gnss_t:
             self.sig_tab[uGNSS.GPS] = {
-                uTYP.C: [rSigRnx('GC1C'), rSigRnx('GC1W'), rSigRnx('GC2W'),
+                uTYP.C: [rSigRnx('GC1C'), rSigRnx('GC2W'),
                          rSigRnx('GC2L'), rSigRnx('GC5Q'), rSigRnx('GC1L')],
-                uTYP.L: [rSigRnx('GL1C'), rSigRnx('GL1W'), rSigRnx('GL2W'),
+                uTYP.L: [rSigRnx('GL1C'), rSigRnx('GL2W'),
                          rSigRnx('GL2L'), rSigRnx('GL5Q'), rSigRnx('GL1L')],
-                uTYP.S: [rSigRnx('GS1C'), rSigRnx('GS1W'), rSigRnx('GS2W'),
+                uTYP.S: [rSigRnx('GS1C'), rSigRnx('GS2W'),
                          rSigRnx('GS2L'), rSigRnx('GS5Q'), rSigRnx('GS1L')],
             }
 
@@ -108,16 +108,18 @@ def __init__(self, opt=None, prefix='', gnss_t='GECJ'):
 
         self.rtcm = rtcm(foutname=foutname)
         self.time_p = gtime_t()
-        self.obs = None
 
         if opt is not None:
             self.init_param(opt=opt, prefix=prefix)
 
     def add_obs(self, obs):
+        sys = list(obs.sig)[0]
+        if sys not in self.sig_tab:
+            return
+
         self.obs.sat = np.hstack((self.obs.sat, obs.sat))
         nsat = len(obs.sat)
 
-        sys = list(obs.sig)[0]
         if sys not in self.obs.sig:
             self.obs.sig[sys] = obs.sig[sys]
 
@@ -132,14 +134,30 @@ def add_obs(self, obs):
         obs.L[np.isnan(obs.L)] = 0.0
 
         for k, sig in enumerate(obs.sig[sys][uTYP.L]):
+
             if sig in self.sig_tab[sys][uTYP.L]:
                 idx = self.sig_tab[sys][uTYP.L].index(sig)
+
                 obs_.P[:, idx] = obs.P[:, k]
                 obs_.L[:, idx] = obs.L[:, k]
                 obs_.S[:, idx] = obs.S[:, k]
                 obs_.D[:, idx] = obs.D[:, k]
                 obs_.lli[:, idx] = obs.lli[:, k]
 
+        if self.useL1CB and sys == uGNSS.QZS:  # QZS L1E -> L1C
+            k = obs_.P[:, 1] != 0
+            obs_.P[k, 0] = obs_.P[k, 1]
+            obs_.L[k, 0] = obs_.L[k, 1]
+            obs_.S[k, 0] = obs_.S[k, 1]
+            obs_.D[k, 0] = obs_.D[k, 1]
+            obs_.lli[k, 0] = obs_.lli[k, 1]
+
+            obs_.P[k, 1] = 0.0
+            obs_.L[k, 1] = 0.0
+            obs_.S[k, 1] = 0.0
+            obs_.D[k, 1] = 0.0
+            obs_.lli[k, 1] = 0
+
         self.obs.P = np.vstack((self.obs.P, obs_.P))
         self.obs.L = np.vstack((self.obs.L, obs_.L))
         self.obs.S = np.vstack((self.obs.S, obs_.S))
@@ -167,10 +185,19 @@ def init_obs(self, time):
         self.obs.sat = np.empty(0, dtype=int)
         self.obs.sig = {}
 
-    def decode(self, buff, len_, sys=[], prn=[]):
+    def decode(self, buff, len_, sys=[], prn=[], scanmode=False):
         """ decode RTCM binary messages """
 
-        _, obs, eph, geph, seph = self.rtcm.decode(buff, len_)
+        _, obs, eph, geph, seph = self.rtcm.decode(buff, scanmode=scanmode)
+
+        if scanmode:
+            self.re.anttype = self.rtcm.ant_desc
+            self.re.rectype = self.rtcm.rcv_type
+            self.re.rec = self.rtcm.rcv_serial
+            self.re.recver = self.rtcm.firm_ver
+            if self.rtcm.pos_arp is not None:
+                self.re.pos = self.rtcm.pos_arp
+            self.re.glo_bias = self.rtcm.glo_bias
 
         if self.flg_rnxobs and obs is not None:
             self.time = obs.time
@@ -188,42 +215,36 @@ def decode(self, buff, len_, sys=[], prn=[]):
 
             if not self.re.rnx_obs_header_sent:
                 for sys in obs.sig:
-                    self.re.sig_tab[sys] = obs.sig[sys]
+                    if sys in self.sig_tab:
+                        self.re.sig_tab[sys] = obs.sig[sys]
 
             self.add_obs(obs)
 
             self.time_p = self.time
 
-        if eph is not None:
-            self.re.rnx_nav_body(eph, self.fh_rnxnav)
-
-        if geph is not None:
-            self.re.rnx_gnav_body(geph, self.fh_rnxnav)
-
-        if seph is not None:
-            self.re.rnx_snav_body(seph, self.fh_rnxnav)
+        if self.flg_rnxnav:
 
+            if eph is not None:
+                sys, prn = sat2prn(eph.sat)
+                if sys in self.sig_tab:
+                    self.re.rnx_nav_body(eph, self.fh_rnxnav)
 
-def decode(f, opt, args):
-
-    print("Decoding {}".format(f))
+            if geph is not None:
+                sys, prn = sat2prn(geph.sat)
+                if sys in self.sig_tab:
+                    self.re.rnx_gnav_body(geph, self.fh_rnxnav)
 
-    bdir, fname = os.path.split(f)
+            if seph is not None:
+                sys, prn = sat2prn(seph.sat)
+                if sys in self.sig_tab:
+                    self.re.rnx_snav_body(seph, self.fh_rnxnav)
 
-    prefix = fname[4:].removesuffix('.rtcm3')+'_'
-    prefix = str(Path(bdir) / prefix) if bdir else prefix
-    rtcmdec = rtcmDec(opt=opt, prefix=prefix, gnss_t=args.gnss)
-    rtcmdec.monlevel = 2
 
-    rtcmdec.rtcm.week = args.weekref
-
-    path = str(Path(bdir) / fname) if bdir else fname
-    blen = os.path.getsize(path)
+def rtcm_decode(rtcmdec, path, blen, scanmode=False):
 
     with open(path, 'rb') as f:
         msg = f.read(blen)
         maxlen = len(msg)-5
-        # maxlen = 400000
         k = 0
         for _ in range(maxlen):
             if k > maxlen:
@@ -237,9 +258,29 @@ def decode(f, opt, args):
                 continue
 
             len_ = rtcmdec.rtcm.len+3
-            rtcmdec.decode(msg[k:k+len_], len_)
+            rtcmdec.decode(msg[k:k+len_], len_, scanmode=scanmode)
             k += rtcmdec.rtcm.dlen
 
+
+def decode(f, opt, args):
+
+    bdir, fname = os.path.split(f)
+
+    prefix = fname[4:].removesuffix('.rtcm3')+'_'
+    prefix = str(Path(bdir) / prefix) if bdir else prefix
+    rtcmdec = rtcmDec(opt=opt, prefix=prefix, gnss_t=args.gnss)
+    rtcmdec.monlevel = 2
+
+    rtcmdec.rtcm.week = args.weekref
+
+    path = str(Path(bdir) / fname) if bdir else fname
+    blen = os.path.getsize(path)
+
+    print(f"Pre-scanning {f}")
+    rtcm_decode(rtcmdec, path, blen, True)
+    print(f"Decoding     {f}")
+    rtcm_decode(rtcmdec, path, blen)
+
     rtcmdec.file_close()
 
 # python decode_rtcm.py ..\data\doy2025-298\CL07298j.rtc --weekref=2389
@@ -270,6 +311,9 @@ def main():
     parser.add_argument("-j", "--jobs", default=int(mp.cpu_count() / 2),
                         type=int, help='Max. number of parallel processes')
 
+    parser.add_argument("--useL1CB", action='store_true',
+                        help="use L1C/B as like L1C/A for QZS")
+
     # Retrieve all command line arguments
     #
     args = parser.parse_args()
@@ -278,10 +322,18 @@ def main():
 
     opt.flg_rnxobs = True
     opt.flg_rnxnav = True
+    opt.useL1CB = args.useL1CB
+
+    # args.weekref = 2397  # 2025/352
+    # args.inpFileName = '..\data\doy2025-352\sept352a.rtc'
+    # args.inpFileName = '../data/doy2025-352/tr92352a.rtc'
+    # args.gnss = 'GJ'
+    # opt.useL1CB = True
 
     s = args.inpFileName
     opt.foutname = s[:s.rfind('.')]+'.log'
 
+    # decode(args.inpFileName, opt, args)
     # Start processing pool
     #
     with mp.Pool(processes=args.jobs) as pool: