Merge pull request #41 from AndreHauschild/devel

Use CODE MGEX products

Author: hirokawa

data/igs_files.txt

@@ -13,17 +13,17 @@ brdc  ftp://gdc.cddis.eosdis.nasa.gov/pub/gnss/data/daily/2023/brdc/BRD400DLR_S_
 brdc  ftp://gdc.cddis.eosdis.nasa.gov/pub/gnss/data/daily/2023/brdc/BRD400DLR_S_20232230000_01D_MN.rnx.gz
 brdc  ftp://gdc.cddis.eosdis.nasa.gov/pub/gnss/data/daily/2023/brdc/BRD400DLR_S_20231890000_01D_MN.rnx.gz
 
-igs   ftp://gdc.cddis.eosdis.nasa.gov/pub/gps/products/2353/COD0OPSFIN_20250460000_01D_05M_ORB.SP3.gz
-igs   ftp://gdc.cddis.eosdis.nasa.gov/pub/gps/products/2353/COD0OPSFIN_20250460000_01D_30S_CLK.CLK.gz
-igs   ftp://gdc.cddis.eosdis.nasa.gov/pub/gps/products/2353/COD0OPSFIN_20250460000_01D_01D_OSB.BIA.gz
+igs   ftp://gdc.cddis.eosdis.nasa.gov/pub/gps/products/2353/COD0MGXFIN_20250460000_01D_05M_ORB.SP3.gz
+igs   ftp://gdc.cddis.eosdis.nasa.gov/pub/gps/products/2353/COD0MGXFIN_20250460000_01D_30S_CLK.CLK.gz
+igs   ftp://gdc.cddis.eosdis.nasa.gov/pub/gps/products/2353/COD0MGXFIN_20250460000_01D_01D_OSB.BIA.gz
 
-igs   ftp://gdc.cddis.eosdis.nasa.gov/pub/gps/products/2274/COD0OPSFIN_20232230000_01D_01D_OSB.BIA.gz
-igs   ftp://gdc.cddis.eosdis.nasa.gov/pub/gps/products/2274/COD0OPSFIN_20232230000_01D_05M_ORB.SP3.gz
-igs   ftp://gdc.cddis.eosdis.nasa.gov/pub/gps/products/2274/COD0OPSFIN_20232230000_01D_30S_CLK.CLK.gz
+igs   ftp://gdc.cddis.eosdis.nasa.gov/pub/gps/products/2274/COD0MGXFIN_20232230000_01D_01D_OSB.BIA.gz
+igs   ftp://gdc.cddis.eosdis.nasa.gov/pub/gps/products/2274/COD0MGXFIN_20232230000_01D_05M_ORB.SP3.gz
+igs   ftp://gdc.cddis.eosdis.nasa.gov/pub/gps/products/2274/COD0MGXFIN_20232230000_01D_30S_CLK.CLK.gz
 
-igs   ftp://gdc.cddis.eosdis.nasa.gov/pub/gps/products/2269/COD0OPSFIN_20231890000_01D_01D_OSB.BIA.gz
-igs   ftp://gdc.cddis.eosdis.nasa.gov/pub/gps/products/2269/COD0OPSFIN_20231890000_01D_05M_ORB.SP3.gz
-igs   ftp://gdc.cddis.eosdis.nasa.gov/pub/gps/products/2269/COD0OPSFIN_20231890000_01D_30S_CLK.CLK.gz
+igs   ftp://gdc.cddis.eosdis.nasa.gov/pub/gps/products/2269/COD0MGXFIN_20231890000_01D_01D_OSB.BIA.gz
+igs   ftp://gdc.cddis.eosdis.nasa.gov/pub/gps/products/2269/COD0MGXFIN_20231890000_01D_05M_ORB.SP3.gz
+igs   ftp://gdc.cddis.eosdis.nasa.gov/pub/gps/products/2269/COD0MGXFIN_20231890000_01D_30S_CLK.CLK.gz
 
 igs   ftp://gdc.cddis.eosdis.nasa.gov/pub/gps/products/2176/COD0MGXFIN_20212650000_01D_01D_OSB.BIA.gz
 igs   ftp://gdc.cddis.eosdis.nasa.gov/pub/gps/products/2176/COD0MGXFIN_20212650000_01D_05M_ORB.SP3.gz

samples/ssr2sp3.py

@@ -20,9 +20,10 @@
 from cssrlib.peph import atxdec
 from cssrlib.peph import peph, peph_t, apc2com
 from cssrlib.cssrlib import sCSSRTYPE as sc
-from cssrlib.cssrlib import cssr
+#from cssrlib.cssrlib import cssr
 from cssrlib.cssr_bds import cssr_bds
 from cssrlib.cssr_has import cssr_has
+from cssrlib.cssr_mdc import cssr_mdc
 from cssrlib.rinex import rnxdec
 
 
@@ -94,14 +95,14 @@ def write_bsx(bsxfile, ac, data):
             f.write(line+'\n')
 
 
-def file2time(year, fileName):
+def file2time(fileName):
     """
     Convert hourly SBF filename to epoch
     """
 
-    doy = os.path.basename(fileName).split('_')[0]
-    hour = ord(doy[3])-ord('a')
-    doy = int(doy[0:3])
+    year = int(os.path.dirname(fileName).split('/')[-1][3:7])
+    doy = int(os.path.dirname(fileName).split('/')[-1][8:11])
+    hour = ord(os.path.basename(fileName)[3])-ord('a')
     time = epoch2time([year, 1, 1, hour, 0, 0])
 
     return timeadd(time, (doy-1)*86400)
@@ -117,27 +118,21 @@ def file2time(year, fileName):
 if len(sys_argv) > 1:
     ssrfiles = sys_argv[1:]
 else:
-    ssrfiles = ['../data/doy2023-189/gale6_189e.txt', ]
+    ssrfiles = ['../data/doy2025-046/046r_gale6.txt', ]
 
 # Start time
 #
-if "_189e" in ssrfiles[0]:
-    time = epoch2time([2023, 7, 8, 4, 0, 0])
-else:
-    time = file2time(2023, ssrfiles[0])
-    """
-    print("ERROR: unknown epoch!")
-    sys_exit(1)
-    """
+time = file2time(ssrfiles[0])
 
 ep = time2epoch(time)
+
 year = ep[0]
 hour = ep[3]
 doy = int(time2doy(time))
 
 if "qzsl6" in ssrfiles[0]:
 
-    name = 'QZS0CLSOPS'
+    name = 'QZS0MDCOPS'
     step = "10S"
 
     dtype = [('wn', 'int'), ('tow', 'int'), ('prn', 'int'),
@@ -219,8 +214,7 @@ def file2time(year, fileName):
     file_gm = baseDirName+'Galileo-HAS-SIS-ICD_1.0_Annex_B_Reed_Solomon_Generator_Matrix.txt'
     gMat = np.genfromtxt(file_gm, dtype="u1", delimiter=",")
 elif 'qzsl6' in ssrfiles[0]:
-    cs = cssr()
-    cs.cssrmode = sc.QZS_MADOCA
+    cs = cssr_mdc()
 elif "bdsb2b" in ssrfiles[0]:
     cs = cssr_bds()
 else:
@@ -412,7 +406,11 @@ def file2time(year, fileName):
 
             # Convert to CoM using ANTEX PCO corrections
             #
-            rs[0, :] += apc2com(nav, sat, time, rs[0, :], sig0, k=0)
+            pco = apc2com(nav, sat, time, rs[0, :], sig0, k=0)
+            if pco is None:
+                continue
+
+            rs[0, :] += pco
 
             for i in range(3):
                 peph.pos[sat-1, i] = rs[0, i]

samples/test_eph.py

@@ -2,40 +2,48 @@
 import cartopy.crs as ccrs
 import numpy as np
 from cssrlib.rinex import rnxdec
-from cssrlib.gnss import Nav, epoch2time, prn2sat, uGNSS, sat2prn,\
+from cssrlib.gnss import Nav, epoch2time, prn2sat, uGNSS, sat2prn, \
     timeadd, ecef2pos
-from cssrlib.ephemeris import findeph, eph2pos
+from cssrlib.ephemeris import findeph, eph2pos, geph2pos
 
 rnx_ver = 4
 
 if rnx_ver == 3:  # RINEX 3
-    navfile = '../data/doy2021-078/30340780.21q'
-    t0 = epoch2time([2021, 3, 19, 0, 0, 0])
+    navfile = '../data/brdc/BRDC00IGS_R_20231890000_01D_MN.rnx'
+    t0 = epoch2time([2023, 7, 8, 0, 0, 0])
     sys_ref = uGNSS.QZS
     prn_ref = 194
     mode_ref = 0  # 0: LNAV, 1: CNAV, 2: CNAV2
 elif rnx_ver == 4:  # RINEX 4
     navfile = '../data/brdc/BRD400DLR_S_20231890000_01D_MN.rnx'
-    t0 = epoch2time([2023, 7, 8, 4, 0, 0])
-    sys_ref = uGNSS.BDS
-    prn_ref = 35
-    mode_ref = 1  # 0: LNAV, 1: CNAV, 2: CNAV2
+    t0 = epoch2time([2023, 7, 8, 0, 0, 0])
+    sys_ref = uGNSS.GLO
+    prn_ref = 7
+    mode_ref = 0  # 0: LNAV, 1: CNAV, 2: CNAV2
 
 dec = rnxdec()
 nav = dec.decode_nav(navfile, Nav())
 
-n = 24*3600//300
+step = 300
+n = 24*3600//step
 flg_plot = True
 
 
 if True:
     t = t0
     sat = prn2sat(sys_ref, prn_ref)
-    eph = findeph(nav.eph, t, sat, mode=mode_ref)
+    if sys_ref == uGNSS.GLO:
+        eph = findeph(nav.geph, t, sat, mode=mode_ref)
+    else:
+        eph = findeph(nav.eph, t, sat, mode=mode_ref)
     if eph is not None:
-        rs, vs, dts = eph2pos(t, eph, True)
+        if sys_ref == uGNSS.GLO:
+            rs, vs, dts = geph2pos(t, eph, True)
+        else:
+            rs, vs, dts = eph2pos(t, eph, True)
 
 if flg_plot:
+
     lon0 = 135
     plt.figure(figsize=(6, 6))
     ax = plt.axes(projection=ccrs.Orthographic(central_longitude=lon0,
@@ -51,11 +59,22 @@
         sys, prn = sat2prn(sat)
         if sys != sys_ref:
             continue
+        """
+        if prn != prn_ref:
+            continue
+        """
         for i in range(n):
-            t = timeadd(t0, i*300)
+            t = timeadd(t0, i*step)
+            if sys_ref == uGNSS.GLO:
+                eph = findeph(nav.geph, t, sat, mode=mode_ref)
+            else:
+                eph = findeph(nav.eph, t, sat, mode=mode_ref)
             if eph is None:
                 continue
-            rs, dts = eph2pos(t, eph)
+            if sys_ref == uGNSS.GLO:
+                rs, dts = geph2pos(t, eph)
+            else:
+                rs, dts = eph2pos(t, eph)
             pos[i, :] = ecef2pos(rs)
             pos[i, 0] = np.rad2deg(pos[i, 0])
             pos[i, 1] = np.rad2deg(pos[i, 1])

samples/test_pntpos.py

@@ -7,36 +7,71 @@
 from sys import stdout
 
 from cssrlib.rinex import rnxdec
-from cssrlib.gnss import ecef2pos, timediff, ecef2enu
-from cssrlib.gnss import rSigRnx, epoch2time, Nav, time2str
+from cssrlib.gnss import Nav, rSigRnx, ecef2pos, ecef2enu
+from cssrlib.gnss import epoch2time, time2doy, time2str, timediff
 from cssrlib.pntpos import stdpos
 
-if False:
-    xyz_ref = [-3962108.673,   3381309.574,   3668678.638]
-    navfile = '../data/doy2021-078/SEPT078M.21P'
-    obsfile = '../data/doy2021-078/SEPT078M.21O'
+# Start epoch and number of epochs
+#
+dataset = 3
+
+if dataset == 0:  # SETP078M.21O
     ep = [2021, 3, 19, 12, 0, 0]
-else:
-    xyz_ref = [-3962108.6726, 3381309.4719, 3668678.6264]
+    xyz_ref = [-3962108.6617, 3381309.5232, 3668678.6410]
+elif dataset == 1:  # SETP1890.23O
+    ep = [2023, 7, 8, 4, 0, 0]
+    xyz_ref = [-3962108.7063, 3381309.5703, 3668678.6690]
+elif dataset == 2:  # SETP223Z.23O
     ep = [2023, 8, 11, 21, 0, 0]
-    navfile = '../data/brdc/BRD400DLR_S_20232230000_01D_MN.rnx'
-    # navfile = '../data/brdc/BRDC00IGS_R_20232230000_01D_MN.rnx'
-    # navfile = '../data/doy2023-223/NAV223.23p'
-    # obsfile = '../data/doy2023-223/SEPT223Z.23O'  # MOSAIC-CLAS
-    obsfile = '../data/doy2023-223/SEPT223Y.23O'  # PolaRX5
-
-pos_ref = ecef2pos(xyz_ref)
-nep = 360
+    xyz_ref = [-3962108.7063, 3381309.5703, 3668678.6690]
+elif dataset == 3:  # 046r_rnx.obs
+    ep = [2025, 2, 15, 17, 0, 0]
+    xyz_ref = [-3962108.6819, 3381309.5707, 3668678.6750]
+else:
+    print("ERROR: no RINEX data set selected!")
+    exit(1)
 
 time = epoch2time(ep)
+year = ep[0]
+doy = int(time2doy(time))
+
+nep = 900*4
+
+bdir = '../data/doy{:04d}-{:03d}/'.format(year, doy)
+
+if dataset == 0:
+    let = chr(ord('A')+ep[3])
+    navfile = bdir+'SEPT{:03d}{}.{:02d}P'.format(doy, let, year % 2000)
+    obsfile = bdir+'SEPT{:03d}{}.{:02d}O'.format(doy, let, year % 2000)
+if dataset == 1:
+    let = '0'
+    navfile = bdir+'SEPT{:03d}{}.{:02d}P'.format(doy, let, year % 2000)
+    obsfile = bdir+'SEPT{:03d}{}.{:02d}O'.format(doy, let, year % 2000)
+elif dataset == 2:
+    # let = 'Z' # MOSAIC-CLAS
+    let = 'Y'  # PolaRX5
+    navfile = '../data/brdc/' + \
+        'BRD400DLR_S_{:04d}{:03d}0000_01D_MN.rnx'.format(year, doy)
+    obsfile = bdir+'SEPT{:03d}{}.{:02d}O'.format(doy, let, year % 2000)
+else:
+    let = chr(ord('a')+ep[3])
+    navfile = bdir+'{:03d}{}_rnx.nav'.format(doy, let)
+    obsfile = bdir+'{:03d}{}_rnx.obs'.format(doy, let)
+
+pos_ref = ecef2pos(xyz_ref)
 
 # Define signals to be processed
 #
-sigs = [rSigRnx("GC1C"), rSigRnx("EC1C"), rSigRnx("JC1C"),
-        rSigRnx("GL1C"), rSigRnx("EL1C"), rSigRnx("JL1C"),
-        rSigRnx("GS1C"), rSigRnx("ES1C"), rSigRnx("JS1C")]
-
-# sigs = [rSigRnx("GC1C"), rSigRnx("GL1C"), rSigRnx("GS1C")]
+gnss = "GEJ"
+sigs = []
+if 'G' in gnss:
+    sigs.extend([rSigRnx("GC1C"), rSigRnx("GL1C"), rSigRnx("GS1C")])
+if 'E' in gnss:
+    sigs.extend([rSigRnx("EC1C"), rSigRnx("EL1C"), rSigRnx("ES1C")])
+if 'C' in gnss:
+    sigs.extend([rSigRnx("CC2I"), rSigRnx("CL2I"), rSigRnx("CS2I")])
+if 'J' in gnss:
+    sigs.extend([rSigRnx("JC1C"), rSigRnx("JL1C"), rSigRnx("JS1C")])
 
 rnx = rnxdec()
 rnx.setSignals(sigs)
@@ -108,8 +143,19 @@
         if obs.t.time == 0:
             break
 
+    # Send line-break to stdout
+    #
+    stdout.write('\n')
+
+    # Close RINEX observation file
+    #
     rnx.fobs.close()
 
+    # Close output file
+    #
+    if nav.fout is not None:
+        nav.fout.close()
+
 
 fig_type = 1
 ylim_h = 2.0

samples/test_pntpos2.py

@@ -20,9 +20,16 @@
 
 # Define signals to be processed
 #
-sigs = [rSigRnx("GC1C"), rSigRnx("EC1C"),
-        rSigRnx("GL1C"), rSigRnx("EL1C"),
-        rSigRnx("GS1C"), rSigRnx("ES1C")]
+gnss = "GE"
+sigs = []
+if 'G' in gnss:
+    sigs.extend([rSigRnx("GC1C"), rSigRnx("GL1C"), rSigRnx("GS1C")])
+if 'E' in gnss:
+    sigs.extend([rSigRnx("EC1C"), rSigRnx("EL1C"), rSigRnx("ES1C")])
+if 'C' in gnss:
+    sigs.extend([rSigRnx("CC2I"), rSigRnx("CL2I"), rSigRnx("CS2I")])
+if 'J' in gnss:
+    sigs.extend([rSigRnx("JC1C"), rSigRnx("JL1C"), rSigRnx("JS1C")])
 
 rnx = rnxdec()
 rnx.setSignals(sigs)

samples/test_pppigs.py

@@ -5,7 +5,6 @@
 import matplotlib.pyplot as plt
 import matplotlib.dates as md
 import numpy as np
-from os.path import exists
 from sys import exit as sys_exit
 from sys import stdout
 
@@ -67,7 +66,7 @@
     navfile = bdir+'{:03d}{}_rnx.nav'.format(doy, let)
     obsfile = bdir+'{:03d}{}_rnx.obs'.format(doy, let)
 
-ac = 'COD0OPSFIN'
+ac = 'COD0MGXFIN'
 
 orbfile = '../data/igs/{}_{:4d}{:03d}0000_01D_05M_ORB.SP3'\
     .format(ac, year, doy)
@@ -78,17 +77,6 @@
 bsxfile = '../data/igs/{}_{:4d}{:03d}0000_01D_01D_OSB.BIA'\
     .format(ac, year, doy)
 
-if not exists(clkfile):
-    orbfile = orbfile.replace('COD0OPSFIN', 'COD0OPSRAP')
-    clkfile = clkfile.replace('COD0OPSFIN', 'COD0OPSRAP')
-    bsxfile = bsxfile.replace('COD0OPSFIN', 'COD0OPSRAP')
-if not exists(clkfile):
-    orbfile = orbfile.replace('COD0OPSRAP', 'COD0MGXFIN')
-    clkfile = clkfile.replace('COD0OPSRAP', 'COD0MGXFIN')
-    bsxfile = bsxfile.replace('COD0OPSRAP', 'COD0MGXFIN')
-if not exists(orbfile):
-    orbfile = orbfile.replace('_05M_', '_15M_')
-
 # Define signals to be processed
 #
 gnss = "GE"
@@ -168,10 +156,10 @@
     #
     ppp = pppos(nav, rnx.pos, 'test_pppigs.log')
     nav.ephopt = 4  # IGS
-    nav.armode = 3
+    nav.armode = 3  # 1: continuous, 3: fix-and-hold
+    nav.thresar = 2.0
 
     nav.elmin = np.deg2rad(10.0)
-    nav.thresar = 2.0
 
     # Get equipment information
     #