diff --git a/data/ewss/jma/tab4_1_2_45.txt b/data/ewss/jma/tab4_1_2_45.txt index ef3266c..6c448d7 100644 --- a/data/ewss/jma/tab4_1_2_45.txt +++ b/data/ewss/jma/tab4_1_2_45.txt @@ -203,6 +203,7 @@ 830304004400 神田川(東京都) 830304004700 妙正寺川(東京都) 830304006400 入間川流域(埼玉県) +830304006403 入間川中流部（埼玉県） 830305000100 多摩川(東京都・神奈川県) 830305000500 野川・仙川(東京都) 830305002000 浅川(東京都) diff --git a/data/ewss/jma/tab4_1_2_6.txt b/data/ewss/jma/tab4_1_2_6.txt index b2d3c5a..901f1a6 100644 --- a/data/ewss/jma/tab4_1_2_6.txt +++ b/data/ewss/jma/tab4_1_2_6.txt @@ -27,6 +27,7 @@ 147 沖合で高い津波を観測したため津波警報を切り替えました。 148 沖合で高い津波を観測したため予想される津波の高さを切り替えました。 149 ただちに避難してください。 +150 南海トラフ地震臨時情報を発表しています。 201 強い揺れに警戒してください。 211 津波警報等(大津波警報・津波警報あるいは津波注意報 )を発表中です。 212 この地震により、日本の沿岸では若干の海面変動があるかもしれませんが、被害の心配はありません。 diff --git a/data/igs_files.txt b/data/igs_files.txt index 66fce1c..dbe6794 100644 --- a/data/igs_files.txt +++ b/data/igs_files.txt @@ -4,9 +4,9 @@ antex https://files.igs.org/pub/station/general/igs14.atx antex https://files.igs.org/pub/station/general/pcv_archive/igs20_2353.atx antex https://www.gsc-europa.eu/sites/default/files/sites/all/files/has14_2345.atx -antex http://ftp.aiub.unibe.ch/CODE_MGEX/CODE/I20.ATX -antex http://ftp.aiub.unibe.ch/CODE_MGEX/CODE/M20.ATX -antex http://ftp.aiub.unibe.ch/CODE_MGEX/CODE/M14.ATX +antex https://www.aiub.unibe.ch/download/CODE_MGEX/CODE/I20.ATX +antex https://www.aiub.unibe.ch/download/CODE_MGEX/CODE/M20.ATX +antex https://www.aiub.unibe.ch/download/CODE_MGEX/CODE/M14.ATX brdc ftp://gdc.cddis.eosdis.nasa.gov/pub/gnss/data/daily/2023/brdc/BRDC00IGS_R_20231890000_01D_MN.rnx.gz brdc ftp://gdc.cddis.eosdis.nasa.gov/pub/gnss/data/daily/2023/brdc/BRDC00IGS_R_20231800000_01D_MN.rnx.gz diff --git a/data/ldpc/H_bds_b1c_sf2.txt b/data/ldpc/H_bds_b1c_sf2.txt new file mode 100644 index 0000000..ea8ab12 --- /dev/null +++ b/data/ldpc/H_bds_b1c_sf2.txt @@ -0,0 +1,50 @@ +11 62 102 150 9 60 100 148 0 51 142 197 22 80 116 154 +4 90 131 177 47 95 138 191 51 79 146 195 44 75 142 190 +13 57 135 198 24 65 120 173 6 88 129 179 7 89 130 176 +6 58 106 158 8 60 108 160 44 92 139 188 4 56 104 156 +10 61 101 149 39 87 123 168 15 67 105 167 50 78 145 194 +17 98 151 187 46 94 137 190 14 66 104 166 7 59 107 159 +21 83 119 153 31 87 114 167 2 49 140 199 12 64 106 164 +40 53 132 159 19 96 149 185 16 68 112 168 14 58 132 199 +34 69 125 162 23 75 119 175 42 96 144 192 8 63 103 151 +23 81 117 155 24 93 111 182 20 72 116 172 17 69 113 169 +34 82 130 182 1 53 101 153 46 73 140 188 13 65 107 165 +2 54 102 154 18 70 114 170 26 67 122 175 29 77 125 177 +36 84 120 169 25 94 108 183 39 89 137 185 21 73 117 173 +28 76 124 176 36 90 138 186 33 68 124 161 12 56 134 197 +29 85 112 165 45 93 136 189 27 64 123 172 28 84 115 164 +25 66 121 174 37 85 121 170 3 50 141 196 48 76 147 192 +35 70 126 163 32 80 128 180 0 52 100 152 43 52 135 158 +35 83 131 183 10 62 110 162 19 71 115 171 15 59 133 196 +33 81 129 181 41 54 133 156 20 82 118 152 38 86 122 171 +30 78 126 178 9 61 109 161 26 95 109 180 45 72 143 191 +1 48 143 198 40 98 146 194 18 99 148 184 5 57 105 157 +41 99 147 195 31 79 127 179 3 55 103 155 22 74 118 174 +37 91 139 187 5 91 128 178 30 86 113 166 43 97 145 193 +16 97 150 186 11 63 111 163 32 71 127 160 42 55 134 157 +38 88 136 184 47 74 141 189 49 77 144 193 27 92 110 181 +35 13 51 60 1 44 53 24 1 45 15 6 45 15 6 1 +1 44 53 24 1 45 15 6 35 46 56 15 6 1 45 15 +15 6 1 45 44 53 24 1 24 1 44 30 1 45 15 6 +30 24 1 44 24 1 44 30 45 15 6 1 17 38 49 11 +24 1 44 30 24 1 44 53 24 1 44 53 30 24 1 44 +33 42 14 24 33 42 14 24 45 15 6 1 1 45 15 6 +30 24 1 44 24 1 44 53 1 44 30 24 57 25 9 41 +1 45 15 6 1 45 15 6 42 36 12 57 6 1 45 15 +24 1 44 53 24 1 44 30 1 45 15 6 1 45 15 6 +44 53 24 1 30 24 1 44 1 44 30 24 53 24 1 44 +1 44 53 24 27 28 30 31 53 24 1 44 24 1 44 30 +45 15 6 1 30 24 1 44 1 45 15 6 26 22 14 2 +35 13 18 60 45 15 6 1 30 1 44 7 6 1 45 15 +6 1 45 15 53 24 1 44 24 1 44 53 30 24 1 44 +1 44 30 24 44 53 24 1 53 24 1 44 44 30 24 1 +30 24 1 44 1 44 30 24 1 44 30 24 41 16 29 51 +1 44 30 24 38 23 22 7 44 53 24 1 1 45 15 6 +30 24 1 44 53 24 1 44 6 1 45 15 24 1 44 53 +35 46 56 15 5 33 42 14 54 7 38 23 1 45 15 6 +44 30 24 1 6 1 45 15 53 24 1 44 44 53 24 1 +1 44 53 24 1 44 30 24 44 30 24 1 1 44 53 24 +45 15 6 1 6 1 45 15 1 44 53 24 42 47 37 32 +51 60 35 13 29 28 30 31 6 1 45 15 24 1 44 53 +44 53 24 1 44 30 24 1 38 49 11 17 44 30 24 1 +24 1 44 30 24 1 44 30 1 44 53 24 53 24 1 44 \ No newline at end of file diff --git a/data/ldpc/H_bds_b1c_sf3.txt b/data/ldpc/H_bds_b1c_sf3.txt new file mode 100644 index 0000000..f8120d6 --- /dev/null +++ b/data/ldpc/H_bds_b1c_sf3.txt @@ -0,0 +1,22 @@ +14 35 56 70 11 29 55 73 13 39 53 69 15 34 57 71 +1 27 45 54 23 41 63 87 2 20 46 68 6 24 50 61 +2 26 61 79 9 33 59 77 4 30 48 74 22 42 59 76 +12 38 52 68 23 43 58 77 19 21 63 64 11 25 65 82 +17 39 44 75 9 35 49 72 19 29 66 84 13 36 56 82 +17 43 67 81 22 40 62 86 3 21 47 69 10 24 64 83 +0 37 70 86 5 31 49 75 4 40 53 84 5 41 52 85 +18 28 67 85 0 26 44 55 10 28 54 72 7 30 50 81 +1 36 71 87 16 38 45 74 8 34 48 73 8 32 58 76 +12 37 57 83 6 31 51 80 15 33 47 79 16 42 66 80 +7 25 51 60 3 27 60 78 14 32 46 78 18 20 62 65 +30 24 1 44 24 1 44 30 40 32 61 18 53 24 1 44 +51 60 35 13 18 15 32 61 15 6 1 45 30 24 1 44 +6 1 45 15 45 15 6 1 1 45 15 6 1 44 53 24 +24 1 44 53 44 30 24 1 34 33 45 36 55 9 34 3 +1 44 53 24 61 47 20 8 53 24 1 44 15 6 1 45 +13 18 60 35 45 15 6 1 24 1 44 53 37 32 52 47 +44 53 24 1 39 36 34 33 44 35 31 50 12 25 36 14 +15 35 46 56 53 24 1 44 1 44 53 24 24 1 44 30 +44 30 24 1 15 6 1 45 30 24 1 44 2 50 22 14 +33 42 14 5 34 3 55 9 44 35 61 50 15 6 1 45 +45 15 6 1 1 44 30 24 6 1 45 15 1 44 53 24 diff --git a/data/ldpc/H_bds_b2a.txt b/data/ldpc/H_bds_b2a.txt new file mode 100644 index 0000000..01dfdd0 --- /dev/null +++ b/data/ldpc/H_bds_b2a.txt @@ -0,0 +1,24 @@ +19 46 49 76 5 29 53 71 17 30 64 72 22 36 59 82 +22 41 68 94 20 44 54 75 9 41 61 86 6 47 60 89 +8 40 60 87 15 26 66 81 19 24 67 95 2 26 50 72 +5 38 70 89 16 34 64 92 21 45 55 74 0 24 48 78 +23 37 58 83 15 43 56 91 18 47 48 77 14 42 57 90 +6 30 54 76 14 27 67 80 17 35 65 93 7 46 61 88 +1 25 49 79 12 45 69 79 18 25 66 94 23 40 69 95 +8 36 51 84 3 38 56 86 0 29 62 85 2 39 57 87 +11 33 59 81 20 43 74 93 13 32 63 91 11 35 52 83 +16 31 65 73 4 28 52 70 1 28 63 84 12 33 62 90 +21 42 75 92 7 31 55 77 9 37 50 85 10 34 53 82 +4 39 71 88 13 44 68 78 3 27 51 73 10 32 58 80 +1 45 15 6 1 44 53 24 45 15 6 1 30 24 1 44 +18 15 32 61 3 55 9 34 35 31 50 44 45 15 6 1 +24 1 44 53 30 24 1 44 32 42 47 37 6 1 45 15 +44 53 24 1 39 36 34 33 44 53 24 1 44 53 24 1 +45 15 6 1 6 1 45 15 24 1 44 53 9 41 57 58 +32 61 18 40 1 45 15 6 22 14 2 50 24 1 44 30 +30 24 1 44 15 46 45 44 45 15 6 1 1 44 30 24 +24 1 44 53 15 6 1 45 53 24 1 44 7 38 23 54 +1 45 15 6 44 53 24 1 57 25 9 41 35 13 51 60 +33 45 36 34 6 1 45 15 6 1 45 15 6 1 45 15 +44 35 31 50 26 27 37 5 24 1 44 30 33 42 14 5 +24 1 44 30 24 1 44 30 1 44 53 24 1 44 30 24 diff --git a/data/ldpc/H_bds_b2b.txt b/data/ldpc/H_bds_b2b.txt new file mode 100644 index 0000000..805684e --- /dev/null +++ b/data/ldpc/H_bds_b2b.txt @@ -0,0 +1,42 @@ +19 67 109 130 27 71 85 161 31 78 96 122 2 44 83 125 +26 71 104 132 30 39 93 154 4 46 85 127 21 62 111 127 +13 42 101 146 18 66 108 129 27 72 100 153 29 70 84 160 +23 61 113 126 8 50 89 131 34 74 111 157 12 44 100 145 +22 60 112 128 0 49 115 151 6 47 106 144 33 53 82 140 +3 45 84 126 38 80 109 147 9 60 96 141 1 43 82 124 +20 77 88 158 37 54 122 159 3 65 104 149 5 47 86 128 +0 42 81 123 32 79 97 120 35 72 112 158 15 57 93 138 +22 75 107 143 24 69 102 133 1 50 116 152 24 57 119 135 +17 59 95 140 7 45 107 145 34 51 83 138 14 43 99 144 +21 77 106 142 16 58 94 139 20 68 110 131 2 48 114 150 +10 52 91 133 25 70 103 134 32 41 95 153 14 56 91 137 +33 73 113 156 28 73 101 154 4 63 102 147 6 48 87 129 +8 46 105 146 30 80 98 121 41 68 119 150 35 52 81 139 +16 63 114 124 13 55 90 136 31 40 94 155 10 61 97 142 +36 56 121 161 29 74 99 155 5 64 103 148 18 75 89 156 +36 78 110 148 19 76 87 157 15 65 116 123 11 53 92 134 +25 58 117 136 39 66 117 151 11 62 98 143 9 51 90 132 +38 55 120 160 7 49 88 130 17 64 115 125 0 0 0 0 +28 69 86 159 23 76 105 141 12 54 92 135 0 0 0 0 +40 67 118 152 37 79 108 149 26 59 118 137 0 0 0 0 +46 45 44 15 15 24 50 37 24 50 37 15 15 32 18 61 +58 56 60 62 37 53 61 29 46 58 18 6 36 19 3 57 +54 7 38 23 51 59 63 47 9 3 43 29 56 8 46 13 +26 22 14 2 63 26 41 12 17 32 58 37 38 23 55 22 +35 1 31 44 44 51 35 13 30 1 44 7 27 5 2 62 +16 63 20 9 27 56 8 43 1 44 30 24 5 26 27 37 +42 47 37 32 38 12 25 51 43 34 48 57 39 9 30 48 +63 13 54 10 2 46 56 35 47 20 33 26 62 54 56 60 +1 21 25 7 43 58 19 49 28 4 52 44 46 44 14 15 +41 48 2 27 49 21 7 35 40 21 44 17 24 23 45 11 +46 25 22 48 13 29 53 61 52 17 24 61 29 41 10 16 +60 24 4 50 32 49 58 19 43 34 48 57 29 7 10 16 +25 11 7 1 32 49 58 19 42 14 24 33 39 56 30 48 +13 27 56 8 53 40 61 18 8 43 27 56 18 40 32 61 +60 48 2 27 50 54 60 62 58 19 32 49 9 3 63 43 +53 35 16 13 23 25 30 16 18 6 61 21 15 1 42 45 +20 16 63 9 27 37 5 26 29 7 10 16 11 60 6 49 +43 47 18 20 42 14 24 33 43 22 41 20 22 15 12 33 +9 41 57 58 5 31 51 30 9 3 63 43 0 0 0 0 +37 53 61 29 6 45 56 19 33 45 36 34 0 0 0 0 +19 24 42 14 1 45 15 6 8 43 27 56 0 0 0 0 diff --git a/data/sc134/msg/MT54_10_DFi209=0.bin b/data/sc134/msg/MT54_10_DFi209=0.bin new file mode 100644 index 0000000..6a1b93b Binary files /dev/null and b/data/sc134/msg/MT54_10_DFi209=0.bin differ diff --git a/data/sc134/msg/MT54_10_DFi209=1.bin b/data/sc134/msg/MT54_10_DFi209=1.bin new file mode 100644 index 0000000..3c025df Binary files /dev/null and b/data/sc134/msg/MT54_10_DFi209=1.bin differ diff --git a/data/sc134/msg/MT54_10_DFi209=2.bin b/data/sc134/msg/MT54_10_DFi209=2.bin new file mode 100644 index 0000000..15e6b88 Binary files /dev/null and b/data/sc134/msg/MT54_10_DFi209=2.bin differ diff --git a/data/sc134/msg/MT54_9.bin b/data/sc134/msg/MT54_9.bin new file mode 100644 index 0000000..e72a6d3 Binary files /dev/null and b/data/sc134/msg/MT54_9.bin differ diff --git a/data/sc134/msg/ssr/ROMAMSG.bin b/data/sc134/msg/ssr/ROMAMSG.bin new file mode 100644 index 0000000..a80cac4 Binary files /dev/null and b/data/sc134/msg/ssr/ROMAMSG.bin differ diff --git a/data/sc134/msg/ssr/ROVRMSG.bin b/data/sc134/msg/ssr/ROVRMSG.bin new file mode 100644 index 0000000..360f4a4 Binary files /dev/null and b/data/sc134/msg/ssr/ROVRMSG.bin differ diff --git a/data/sc134/msg/ssr/SSRTEST_20260206_CORR_v123.bin b/data/sc134/msg/ssr/SSRTEST_20260206_CORR_v123.bin new file mode 100644 index 0000000..ce54f7e Binary files /dev/null and b/data/sc134/msg/ssr/SSRTEST_20260206_CORR_v123.bin differ diff --git a/receiver/cssr2rtcm.py b/receiver/cssr2rtcm.py new file mode 100644 index 0000000..d62b7ec --- /dev/null +++ b/receiver/cssr2rtcm.py @@ -0,0 +1,179 @@ +""" +Compact SSR messages to RTCM 3 messages converter + +[1] RTCM Standard 10403.4 Differential GNSS Services - + Version 3 with Amendment 1, November, 2024 + +@author: Rui Hirokawa +""" + +import argparse +from glob import glob +import multiprocessing as mp +import numpy as np +from binascii import unhexlify + +from cssrlib.rtcm import rtcme +from cssrlib.cssrlib import cssr, sCType +from cssrlib.gnss import load_config, char2sys + +config = load_config('config_ppprtk.yml') + + +def decode_msg(v, tow, prn_ref, l6_ch=0, prn_ref_ext=0, l6_ch_ext=0): + """ find valid correction message """ + + msg, msg_e = None, None + + vi_ = v[v['tow'] == tow] + vi = vi_[(vi_['type'] == l6_ch) & (vi_['prn'] == prn_ref)] + if len(vi) > 0: + msg = unhexlify(vi['nav'][0]) + + # load regional STEC info (experimental) + if prn_ref_ext > 0: + vi = vi_[(vi_['type'] == l6_ch_ext) & + (vi_['prn'] == prn_ref_ext)] + if len(vi) > 0: + msg_e = unhexlify(vi['nav'][0]) + + return msg, msg_e + + +def encode_msg(cs, re, msgtype, maxlen=1024): + """ encode RTCM message """ + + k = 0 + buff = bytearray(maxlen) + msg = bytearray(maxlen) + + re.msgtype = msgtype + re.udi = cs.udi + re.datum = cs.datum + re.nsat_n = cs.nsat_n + re.sat_n = cs.sat_n + re.lc = cs.lc + re.grid = cs.grid + i = re.encode(buff) + + len_ = (i+7)//8 + if k+len_+6 >= maxlen: + return -1 + + re.set_body(msg, buff, k, len_) + re.set_sync(msg, k) + re.set_len(msg, k, len_) + re.set_checksum(msg, k) + k += re.dlen + + return msg[:k] + + +def out_ssr(cs, re, fc, sct, sys_t=None, inet=0): + """ output SSR messages """ + + if cs.lc[inet].cstat & (1 << sct) == 0: + if sct not in [sCType.META, sCType.GRID]: + return + + if sys_t is None: + mt = re.sct2mt(sct) + msg = encode_msg(cs, re, mt) + fc.write(msg) + else: + for sys in sys_t: + mt = re.sct2mt(sct, sys) + if mt < 0: + return + msg = encode_msg(cs, re, mt) + fc.write(msg) + cs.lc[inet].cstat ^= (1 << sct) + + +def process(infile, outfile, args): + dtype = [('wn', 'int'), ('tow', 'int'), ('prn', 'int'), + ('type', 'int'), ('len', 'int'), ('nav', 'S500')] + v = np.genfromtxt(args.inpFileName, dtype=dtype) + + griddef = config['griddef'] + + cs = cssr() + cs.monlevel = config['cs']['monlevel'] + cs.read_griddef(griddef) + + fc = open(outfile, 'wb') + if not fc: + print("RTCM message file cannot open.") + + re = rtcme() + re.gtype = 1 + + re.gid = args.gid + prn_ref = args.prnref + l6_ch = args.l6ch # L6D + tow = v[0]['tow']-1 + nep = 3600 + + re.inet = re.gid + # maxlen = len(cs.buff) + + sys_t = char2sys(args.gnss) + + for ne in range(nep): + tow += 1 + msg_, _ = decode_msg(v, tow, prn_ref, l6_ch) + re.tow = tow + if msg_ is not None: + cs.decode_l6msg(msg_, 0) + if cs.fcnt == 5: # end of sub-frame + cs.decode_cssr(bytes(cs.buff), 0) + + if ne == 0: + out_ssr(cs, re, fc, sCType.META) + out_ssr(cs, re, fc, sCType.GRID) + + if cs.lc[0].cstat & (1 << sCType.MASK) == 0: + continue + + out_ssr(cs, re, fc, sCType.CLOCK, sys_t) + out_ssr(cs, re, fc, sCType.ORBIT, sys_t) + out_ssr(cs, re, fc, sCType.URA, sys_t) + out_ssr(cs, re, fc, sCType.CBIAS, sys_t) + out_ssr(cs, re, fc, sCType.PBIAS, sys_t, inet=re.inet) + out_ssr(cs, re, fc, sCType.TROP, inet=re.inet) + out_ssr(cs, re, fc, sCType.STEC, sys_t, inet=re.inet) + + fc.close() + + +if __name__ == "__main__": + + # Parse command line arguments + # + parser = argparse.ArgumentParser( + description="QZS L6 (CSSR) to RTCM SSR converter") + + parser.add_argument("inpFileName", + help="Input QZS L6 file(s) (wildcards allowed)") + parser.add_argument("-g", "--gnss", default='GEJ', + help="GNSS [GEJ]") + parser.add_argument("--prnref", type=int, default=199, + help="QZS satellite PRN [193-210] (default:199)") + parser.add_argument("--l6ch", type=int, default=0, + help="QZS satellite L6 channel [0|1] (default:0)") + parser.add_argument("--gid", type=int, default=7, + help="Network ID [1-12] (default:7)") + parser.add_argument("-j", "--jobs", default=int(mp.cpu_count() / 2), + type=int, help='Max. number of parallel processes') + + args = parser.parse_args() + + # args.inpFileName = '../data/doy2025-233/233h_qzsl6.txt' + foutname = args.inpFileName.removesuffix('.txt')+'.rtcm3' + + # process(args.inpFileName, foutname, args) + # Start processing pool + # + with mp.Pool(processes=args.jobs) as pool: + pool.starmap(process, [(f, foutname, args) + for f in glob(args.inpFileName)]) diff --git a/receiver/decode_jps.py b/receiver/decode_jps.py index 23c2b26..f5eb34f 100755 --- a/receiver/decode_jps.py +++ b/receiver/decode_jps.py @@ -127,7 +127,7 @@ class jps(rcvDec): [1, 6, 2, 4, 3, 0], [1, 1, 4, 2, 3, 1], [1, 6, 2, 3, 3, 1], [0, 0, 0, 0, 3, 1]] - sys_t = { + gnss2sys_t = { GNSS.GPS: uGNSS.GPS, GNSS.GLO: uGNSS.GLO, GNSS.GAL: uGNSS.GAL, GNSS.BDS: uGNSS.BDS, GNSS.QZS: uGNSS.QZS, GNSS.SBS: uGNSS.SBS, GNSS.IRN: uGNSS.IRN, GNSS.GLO_C: uGNSS.GLO, @@ -352,13 +352,15 @@ def decode_obs(self): else: prn = self.prn[k] - sys = self.sys_t[self.sys[k]] + sys = self.gnss2sys_t[self.sys[k]] if sys not in self.sig_tab.keys(): continue if sys == uGNSS.GLO and prn == 255: continue if sys == uGNSS.SBS and self.prn[k] > 156: continue + if sys == uGNSS.IRN: + pass sat = prn2sat(sys, prn) if sat in obs.sat: jn = obs.sat.index(sat) @@ -398,9 +400,15 @@ def decode_obs(self): obs.D[np.isnan(obs.D)] = 0 obs.S[np.isnan(obs.S)] = 0 + obs.sat = np.array(obs.sat, dtype=int) + obs.sort() + return obs def decode_nd(self, buff, sys=uGNSS.GPS): + if sys not in self.sys_t: + return 0 + prn, time_, type_, len_ = st.unpack_from('= 2: @@ -542,15 +550,14 @@ def decode(self, buff, len_, sys=[], prn=[]): if self.monlevel > 1: print("[xd] prn={:d} tow={:d} type={:d} len={:d}". format(prn, time_, type_, len_)) - if self.week >= 0: - if self.flg_qzsl6: - if self.prn_ref > 0 and prn != self.prn_ref: - return - msg_l6 = buff[12:12+len_] - self.fh_qzsl6.write( - "{:4d}\t{:6d}\t{:3d}\t{:1d}\t{:3d}\t{:s}\n". - format(self.week, time_, prn, type_, len_, - hexlify(msg_l6).decode())) + if self.week >= 0 and self.flg_qzsl6: + if self.prn_ref > 0 and prn != self.prn_ref: + return + msg_l6 = buff[12:12+len_] + self.fh_qzsl6.write( + "{:4d}\t{:6d}\t{:3d}\t{:1d}\t{:3d}\t{:s}\n". + format(self.week, time_, prn, type_, len_, + hexlify(msg_l6).decode())) # errCorr = st.unpack_from('= 59: # B2b: BDS PPP self.fh_bdsb2b.write( - "{:4d}\t{:6d}\t{:3d}\t{:1d}\t{:3d}\t{:s}\n". - format(self.week, time_, prn, type_, len_*4, + "{:4d}\t{:6d}\t{:3d}\t{:3d}\t{:s}\n". + format(self.week, time_, prn, len_*4, hexlify(b).decode())) elif head == 'gd': # GPS Navigation data @@ -606,6 +615,8 @@ def decode(self, buff, len_, sys=[], prn=[]): b = bytes(np.array(msg, dtype='uint32')) if self.flg_rnxnav: + if uGNSS.IRN not in self.sys_t: + return 0 eph = None if type_ == 0: eph = self.rn.decode_irn_lnav(self.week, time_, sat, b) @@ -648,6 +659,8 @@ def decode(self, buff, len_, sys=[], prn=[]): bs.pack_into('u2', buff, 25*k, (d >> 23) & 0x3) if self.flg_rnxnav and type_ == 0: + if uGNSS.GLO not in self.sys_t: + return 0 geph = self.rn.decode_glo_fdma( self.week, self.tow, sat, buff, fcn) @@ -665,6 +678,8 @@ def decode(self, buff, len_, sys=[], prn=[]): b = bytes(np.array(msg, dtype='uint32')) if self.flg_rnxnav: + if uGNSS.GLO not in self.sys_t: + return 0 geph = None if type_ == 0: # L1OC geph = self.rn.decode_glo_l1oc(self.week, self.tow, sat, b) @@ -699,6 +714,8 @@ def decode(self, buff, len_, sys=[], prn=[]): if type_ == 0 or type_ == 2: # INAV b = buff[12:] if self.flg_rnxnav: + if uGNSS.GAL not in self.sys_t: + return 0 eph = self.rn.decode_gal_inav( self.week, time_, sat, type_, b) if eph is not None: @@ -711,20 +728,22 @@ def decode(self, buff, len_, sys=[], prn=[]): elif type_ == 1: # FNAV b = buff[12:] if self.flg_rnxnav: + if uGNSS.GAL not in self.sys_t: + return 0 eph = self.rn.decode_gal_fnav( self.week, time_, sat, type_, b) if eph is not None: self.re.rnx_nav_body(eph, self.fh_rnxnav) if self.flg_galfnav and self.week >= 0: self.fh_galfnav.write( - "{:4d}\t{:6d}\t{:3d}\t{:1d}\t{:3d}\t{:s}\n" - .format(self.week, time_, prn, type_, len_, + "{:4d}\t{:6d}\t{:3d}\t{:3d}\t{:s}\n" + .format(self.week, time_, prn, len_, hexlify(b).decode())) elif type_ == 6: # CNAV if self.flg_gale6 and self.week >= 0: self.fh_gale6.write( - "{:4d}\t{:6d}\t{:3d}\t{:1d}\t{:3d}\t{:s}\n" - .format(self.week, time_, prn, type_, len_, + "{:4d}\t{:6d}\t{:3d}\t{:3d}\t{:s}\n" + .format(self.week, time_, prn, len_, hexlify(buff[12:]).decode())) elif head == 'WD': # SBAS Navigation data @@ -857,7 +876,7 @@ def decode(self, buff, len_, sys=[], prn=[]): nsat = (len_-6)//8 pr_ = np.array(st.unpack_from('d'*nsat, buff, 5)) if self.navic_work_around: - pr_[pr_ < 0 & (np.array(self.sys) == GNSS.IRN)] = np.nan + pr_[((pr_ < 0) | (pr_ > 0.2)) & (np.array(self.sys) == GNSS.IRN)] = np.nan if head[1] == 'X': # [RX] self.PR_REF[:nsat] = pr_ else: @@ -869,7 +888,7 @@ def decode(self, buff, len_, sys=[], prn=[]): nsat = (len_-6)//8 cp_ = np.array(st.unpack_from('d'*nsat, buff, 5)) if self.navic_work_around: - cp_[cp_ < 0 & (np.array(self.sys) == GNSS.IRN)] = np.nan + cp_[((cp_ < 0) | (cp_ > 1.0)) & (np.array(self.sys) == GNSS.IRN)] = np.nan self.cp[:nsat, ch] = cp_ elif head[0] == 'c' and head[1] in self.ch_t.keys(): @@ -953,7 +972,7 @@ def decode(f, opt, args): bdir, fname = os.path.split(f) - prefix = fname[4:].removesuffix('.jps')+'_' + prefix = fname.removesuffix('.jps') prefix = str(Path(bdir) / prefix) if bdir else prefix jpsdec = jps(opt=opt, prefix=prefix, gnss_t=args.gnss) jpsdec.monlevel = 1 @@ -1003,12 +1022,15 @@ def main(): parser.add_argument("--antenna", default='unknown', help="Antenna type [unknown]") - parser.add_argument("-g", "--gnss", default='GRECIJ', - help="GNSS [GRECIJ]") + parser.add_argument("-g", "--gnss", default='GRECJ', + help="GNSS [GRECJ]") 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() @@ -1020,7 +1042,7 @@ def main(): opt.flg_gale6 = True opt.flg_galinav = True - opt.flg_galfnav = True + opt.flg_galfnav = False opt.flg_qzsl6 = True @@ -1031,6 +1053,11 @@ def main(): opt.flg_gpslnav = True + opt.useL1CB = args.useL1CB + + #args.inpFileName = "../data/doy2026-038/jav3038a.jps" + #decode(args.inpFileName, opt, args) + # Start processing pool # with mp.Pool(processes=args.jobs) as pool: diff --git a/receiver/decode_lgr.py b/receiver/decode_lgr.py new file mode 100644 index 0000000..5ed6773 --- /dev/null +++ b/receiver/decode_lgr.py @@ -0,0 +1,455 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Qascom LuGRE messages decoder + + [1] Quascom, LuGRE Receiver Interface Control Document Issue 2.0, 2025 + +@author Rui Hirokawa +""" + +from skyfield.framelib import itrs +from skyfield.api import load +import argparse +from glob import glob +import multiprocessing as mp +import numpy as np +import os +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, time2gpst +from cssrlib.rawnav import rcvDec, rcvOpt + + +class navRec(): + """ class for navigation record """ + + def __init__(self): + self.t = gtime_t() + self.nsat = 0 + self.pos = np.zeros(3) + self.vel = np.zeros(3) + self.std = np.zeros(3) + self.clkb = 0.0 + self.clkd = 0.0 + self.ggto = 0.0 + self.dops = np.zeros(5) + + +class lgr(rcvDec): + """ class for LuGRE Binary Format decoder """ + tow = -1 + week = -1 + + def __init__(self, opt=None, prefix='', gnss_t='GE'): + super().__init__(opt, prefix, gnss_t) + + sig_tbl = { + uGNSS.GPS: {0: uSIG.L1C, 1: uSIG.L5Q}, + uGNSS.GAL: {2: uSIG.L1C, 3: uSIG.L5Q, 4: uSIG.L7Q} + } + + obs_tbl = { + uGNSS.GPS: [uSIG.L1C, uSIG.L5Q], + uGNSS.GAL: [uSIG.L1C, uSIG.L5Q, uSIG.L7Q] + } + + self.sig_t = {} + for sys in sig_tbl.keys(): + if sys not in self.sig_tab.keys(): + continue + + self.sig_t[sys] = {} + for key in sig_tbl[sys].keys(): + sig = sig_tbl[sys][key] + self.sig_t[sys][key] = { + uTYP.C: rSigRnx(sys, uTYP.C, sig), + uTYP.L: rSigRnx(sys, uTYP.L, sig), + uTYP.D: rSigRnx(sys, uTYP.D, sig), + uTYP.S: rSigRnx(sys, uTYP.S, sig), + } + + for sys in sig_tbl.keys(): + if sys not in self.sig_tab.keys(): + continue + self.sig_tab[sys][uTYP.C] = [] + self.sig_tab[sys][uTYP.L] = [] + self.sig_tab[sys][uTYP.D] = [] + self.sig_tab[sys][uTYP.S] = [] + for sig in obs_tbl[sys]: + self.sig_tab[sys][uTYP.C].append(rSigRnx(sys, uTYP.C, sig)) + self.sig_tab[sys][uTYP.L].append(rSigRnx(sys, uTYP.L, sig)) + self.sig_tab[sys][uTYP.D].append(rSigRnx(sys, uTYP.D, sig)) + self.sig_tab[sys][uTYP.S].append(rSigRnx(sys, uTYP.S, sig)) + + 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' + + def msg_len(self, msg, k): + self.len = st.unpack_from(' 0 and k+len_ >= maxlen: + return False + cs = Crc24LteA.calc(msg[k:k+len_-3]) + cs_ = msg[len_-1] << 16 | msg[len_-2] << 8 | msg[len_-3] + + # if self.monlevel > 0 and cs != cs_: + # print(f"checksum error: len={len_}") + self.len = len_ + self.dlen = len_ + # return cs == cs_ + return True + + def svid2prn(self, svid, sigid): + """ convert from svid/sigid to sys/prn """ + sys = uGNSS.GPS if sigid <= 1 else uGNSS.GAL + prn = svid + return sys, prn + + def decode_iqs(self, buff, k=10): + """ decode I/Q Sample message """ + trx = st.unpack_from('> 4) & 0x3fff + bid = (blk >> 18) & 0x3fff + + # stype 0:1ch,real,1:1ch,complex,2:2ch,real,3:2ch,complex + # ns*qb/4 is stype = complex, ceil(ns*qb/8) otherwise + if stype == 1 or stype == 3: + sz = ns*qb/4 + else: + sz = int(np.ceil(ns*qb/8)) + + iqSamples = buff[k:k+sz] + + def decode_nav(self, buff, k=10): + """ decode NAV message """ + trx = st.unpack_from(' nsig_max: + nsig_max = len(self.sig_tab[s][uTYP.L]) + + self.nsig[uTYP.C] = nsig_max + self.nsig[uTYP.L] = nsig_max + self.nsig[uTYP.D] = nsig_max + self.nsig[uTYP.S] = nsig_max + + obs.sat = np.empty(0, dtype=np.int32) + + trx, nm = st.unpack_from(' 1: + print("skip code={:}".format(sig)) + continue + idx = self.sig_tab[sys][sig.typ].index(sig) + wavelength = rCST.CLIGHT/sig.frequency() + + if sat not in pr.keys(): + pr[sat] = {} + cp[sat] = {} + dp[sat] = {} + lli[sat] = {} + cn[sat] = {} + + lli_ = 0 + + pr[sat][idx] = pr_ + cp[sat][idx] = cp_/wavelength # [m]->[cycle] + dp[sat][idx] = dop_ + lli[sat][idx] = lli_ + cn[sat][idx] = cn_ + + if sat not in obs.sat: + obs.sat = np.append(obs.sat, sat) + + # print(f"OBS {gnss}:{prn:3d} {svid} {sigid:2d}") + + nsat = len(obs.sat) + if nsat == 0: + return None + + obs.sat.sort() + obs.P = np.zeros((nsat, self.nsig[uTYP.C]), dtype=np.float64) + obs.L = np.zeros((nsat, self.nsig[uTYP.L]), dtype=np.float64) + obs.D = np.zeros((nsat, self.nsig[uTYP.D]), dtype=np.float64) + obs.S = np.zeros((nsat, self.nsig[uTYP.S]), dtype=np.float64) + obs.lli = np.zeros((nsat, self.nsig[uTYP.L]), dtype=np.int32) + + for k, sat in enumerate(obs.sat): + for i in pr[sat].keys(): + obs.P[k][i] = pr[sat][i] + obs.L[k][i] = cp[sat][i] + obs.D[k][i] = dp[sat][i] + obs.S[k][i] = cn[sat][i] + obs.lli[k][i] = lli[sat][i] + + return obs + + def decode(self, buff, len_, sys=[], prn=[]): + mt = buff[1:4] + k = 4 + sender, dlen = st.unpack_from('= maxlen: + break + + if not lgrdec.checksum(msg, k): + k += 1 + continue + + lgrdec.decode(msg[k:k+len_], len_) + k += len_ + + nep += 1 + if nep_max > 0 and nep >= nep_max: + break + + 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 + # + parser = argparse.ArgumentParser(description="LuGRE data converter") + + # Input file and folder + # + parser.add_argument( + "inpFileName", help="Input BIN file(s) (wildcards allowed)") + + parser.add_argument("--receiver", default='unknown', + help="Receiver type [unknown]") + parser.add_argument("--antenna", default='unknown', + help="Antenna type [unknown]") + + parser.add_argument("-g", "--gnss", default='GE', + help="GNSS [GE]") + + parser.add_argument("-j", "--jobs", default=int(mp.cpu_count() / 2), + type=int, help='Max. number of parallel processes') + + # Retrieve all command line arguments + # + args = parser.parse_args() + + opt = rcvOpt() + + opt.flg_rnxobs = True + + # args.inpFileName = '../data/doy2025-074/TLM_RAW_20250315_130709_26H_S_OP76_0.bin' + # args.inpFileName = '../data/doy2025-074/TLM_NAV_20250315_130709_26H_S_OP76_0.bin' + + # decode(args.inpFileName, opt, args) + # Start processing pool + # + with mp.Pool(processes=args.jobs) as pool: + pool.starmap(decode, [(f, opt, args) for f in glob(args.inpFileName)]) + + +# Call main function +# +if __name__ == "__main__": + main() diff --git a/receiver/decode_nov.py b/receiver/decode_nov.py index c1a17cc..d25468a 100644 --- a/receiver/decode_nov.py +++ b/receiver/decode_nov.py @@ -181,7 +181,7 @@ def decode_obs(self, buff): # prn -= 37 sat = prn2sat(sys, prn) - gfrq, pr_, sig_pr, cp_, sig_cp, dop_, cn0, lockt = \ + gfrq, pr_, sig_pr, adr_, sig_adr, dop_, cn0, lockt = \ st.unpack_from(' 0: - print(f"week={week} tow={tow:6.1f} id={id_:2d}") if id_ == 7: # GPS L1 C/A ephemeris pass + + elif id_ == 8: # IONUTC + pass + elif id_ == 25: # RAWGPSSUBFRAME if self.flg_gpslnav: k = self.head_len @@ -312,9 +316,22 @@ def decode(self, buff, len_, sys=[], prn=[]): if obs is not None: self.re.rnx_obs_header(obs.time, self.fh_rnxobs) self.re.rnx_obs_body(obs, self.fh_rnxobs) + + elif id_ == 47: # PSRPOS: Pseudorange position + pass + elif id_ == 93: # RXSTATUS pass + elif id_ == 128: # RXCONFIG + pass + + elif id_ == 243: # PSRXYZ: Pseudorange Cartesian position and velocity + pass + + elif id_ == 719: # GLOCLOCK: GLONASS L1 C/A clock information + pass + elif id_ == 722: # GLORAWSTRING if self.flg_gloca: k = self.head_len @@ -348,35 +365,46 @@ def decode(self, buff, len_, sys=[], prn=[]): elif id_ == 1336: # GLONASS Eph pass elif id_ == 1413: # GALFNAVRAWPAGE + if uGNSS.GAL not in self.sys_t: + return 0 + k = self.head_len ch, prn = st.unpack_from(' 5 and prn < 59: eph = self.rn.decode_bds_d1( self.week, self.tow, sat, msg) @@ -433,7 +461,10 @@ def decode(self, buff, len_, sys=[], prn=[]): elif id_ == 1696: # BDS Eph pass elif id_ == 2105: # NAVICRAWSUBFRAME - if self.flg_irnnav: + if uGNSS.IRN not in self.sys_t: + return 0 + + if self.flg_rnxnav: k = self.head_len ch, prn, sid = st.unpack_from(' 0: - self.sbas_frm = {} - self.time_p = self.time - k = self.head_len - prn, ch, src, pre, _, frm = st.unpack_from(' 0: + self.sbas_frm = {} - itype = 0 if src == 1 else 1 # 0:L1, 1:L5 + self.time_p = self.time + k = self.head_len + prn, ch, src, pre, _, frm = st.unpack_from(' 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)) @@ -158,52 +185,66 @@ 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, scanmode=scanmode) - _, obs, eph, geph, seph = self.rtcm.decode(buff, len_) + 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 - self.re.rnx_obs_header(obs.time, self.fh_rnxobs) if timediff(self.time, self.time_p) != 0.0: + + if self.time_p.time > 0: + self.re.rnx_obs_header(self.obs.time, self.fh_rnxobs) + if self.obs is not None: - self.re.rnx_obs_body(self.obs, self.fh_rnxobs) + if self.re.rnx_obs_header_sent: + self.re.rnx_obs_body(self.obs, self.fh_rnxobs) + self.init_obs(obs.time) + if not self.re.rnx_obs_header_sent: + for sys in obs.sig: + 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 = 1 - rtcmdec.rtcm.week = args.weekref +def rtcm_decode(rtcmdec, path, blen, scanmode=False): - path = str(Path(bdir) / fname) if bdir else fname - blen = os.path.getsize(path) 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: @@ -217,11 +258,33 @@ 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 + def main(): @@ -248,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() @@ -256,7 +322,20 @@ 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.weekref = 2380 + # args.inpFileName = '../data/doy2025-233/233h_qzsl6.rtcm3' + # 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: diff --git a/receiver/decode_sbf.py b/receiver/decode_sbf.py index 0bb2b74..45fa2f2 100755 --- a/receiver/decode_sbf.py +++ b/receiver/decode_sbf.py @@ -3,11 +3,14 @@ """ Septentrio Receiver SBF messages decoder - [1] mosaic-X5 Reference Guide, Applicable to version 4.14.10.1 - of the Firmware, 2024 + [1] mosaic-X5 Reference Guide, Applicable to version 4.15.0 + of the Firmware, July, 2025 - [1] PolaRX5 Reference Guide, Applicable to version 5.6.0 - of the Firmware, 2025 + [2] mosaic-G5 Reference Guide, Applicable to version 1.0.1 + of the Firmware, Novemrber, 2025 + + [3] PolaRX5 Reference Guide, Applicable to version 5.7.0 + of the Firmware, December, 2025 @author Rui Hirokawa """ @@ -131,37 +134,40 @@ def svid2prn(self, svid): elif svid <= 61: # R1-R24 sys = uGNSS.GLO prn = svid-37 - elif svid <= 62: # GLONASS, unknown slot + elif svid == 62: # GLONASS, unknown slot sys = uGNSS.GLO prn = 0 elif svid <= 68: # R25-R30 sys = uGNSS.GLO prn = svid-38 - elif svid >= 71 and svid <= 106: # E1-E36 + elif svid <= 70: # reserved (69-70) + sys = uGNSS.NONE + prn = 0 + elif svid <= 106: # E1-E36 sys = uGNSS.GAL prn = svid-70 - elif svid >= 71 and svid <= 119: # L-Band(MSS) + elif svid <= 119: # L-Band(MSS) sys = uGNSS.NONE prn = 0 - elif svid >= 120 and svid <= 140: # S120-S140 + elif svid <= 140: # S120-S140 sys = uGNSS.SBS prn = svid - elif svid >= 141 and svid <= 180: # C1-C40 + elif svid <= 180: # C1-C40 sys = uGNSS.BDS prn = svid-140 - elif svid >= 181 and svid <= 190: # J1-J10 + elif svid <= 190: # J1-J10 sys = uGNSS.QZS prn = svid-180+192 - elif svid >= 191 and svid <= 197: # I1-I7 + elif svid <= 197: # I1-I7 sys = uGNSS.IRN prn = svid-190 - elif svid >= 198 and svid <= 215: # S141-S158 + elif svid <= 215: # S141-S158 sys = uGNSS.SBS prn = svid-57 - elif svid >= 216 and svid <= 222: # I8-I14 + elif svid <= 222: # I8-I14 sys = uGNSS.IRN prn = svid-208 - elif svid >= 223 and svid <= 245: # C41-C63 + elif svid <= 245: # C41-C63 sys = uGNSS.BDS prn = svid-182 else: # reserved (246-255) @@ -260,22 +266,25 @@ def decode_obs(self, buff, k=8): else: S1 = cn0*0.25 + 10.0 - if code not in self.sig_tab[sys][code.typ]: - if self.monlevel > 1: - print("skip code={:}".format(code)) - k += nb2*sb2len - continue - idx = self.sig_tab[sys][code.typ].index(code) + if code in self.sig_tab[sys][code.typ]: + idx = self.sig_tab[sys][code.typ].index(code) - pr[idx] = P1 - cp[idx] = L1 - dp[idx] = D1 - ll[idx] = lli - cn[idx] = S1 + if code == rSigRnx('JL1E'): # for testing + idx = 0 - if self.monlevel >= 2: - print("{:6d} {:3d} {:} {:14.3f} {:14.3f} {:10.3f} {:4.2f}". - format(int(self.tow), prn, code, P1, L1, D1, S1)) + pr[idx] = P1 + cp[idx] = L1 + dp[idx] = D1 + ll[idx] = lli + cn[idx] = S1 + + if self.monlevel >= 2: + print("{:6d} {:3d} {:} {:14.3f} {:14.3f} {:10.3f} {:4.2f}". + format(int(self.tow), prn, code, P1, L1, D1, S1)) + + else: + if self.monlevel > 1: + print(f"skip code={code} for sys={sys}") for j in range(nb2): typ, ltime, cn0, ofst1, cp1, info, cofst0, cp0, dop0 = \ @@ -350,6 +359,8 @@ def decode_obs(self, buff, k=8): obs.S = obs.S.reshape(len(obs.sat), self.nsig[uTYP.S]) obs.lli = obs.lli.reshape(len(obs.sat), self.nsig[uTYP.L]) + obs.sort() + return obs def decode_gpsnav(self, buff, k=8): @@ -467,7 +478,8 @@ def decode(self, buff, len_, sys=[], prn=[]): if self.monlevel > 1 and blk_num not in \ (4002, 4004, 4006, 4007, 4017, 4018, 4019, 4020, 4021, 4022, 4023, 4024, 4026, 4027, 4036, 4047, 4066, 4067, 4068, 4069, 4081, 4093, - 4095, 4218, 4219, 4228, 4242, 4246, 4262, 5891, 5894, 5896): + 4095, 4218, 4219, 4227, 4228, 4242, 4246, 4262, 4270, 4271, 5891, + 5894, 5896): print("block_num = {:d} rev={:d} len={:d}".format( blk_num, blk_rev, len_)) @@ -497,65 +509,82 @@ def decode(self, buff, len_, sys=[], prn=[]): self.re.pos = pos2ecef([x, y, z]) elif blk_num in (4017, 4066): # GPSRawCA, QZSRawCA + if (blk_num == 4017 and uGNSS.GPS not in self.sys_t) or \ + (blk_num == 4066 and uGNSS.QZS not in self.sys_t): + return 0 + sys, prn = self.decode_head(buff, k) sat = prn2sat(sys, prn) k += 7 crcpass, _, src, _, ch = st.unpack_from(' 2: + print("crc error in GPSRawCA/QZSRawCA " + + "{:6d}\t{:2d}\t{:1d}\t{:2d}". + format(int(self.tow), prn, crcpass, src)) + return -1 + + msg = bytearray(40) + for i in range(10): + d = st.unpack_from('L', msg, i*4, d) + k += 4 + msg = bytes(msg) + if (sys == uGNSS.GPS and self.flg_gpslnav) or \ (sys == uGNSS.QZS and self.flg_qzslnav): - if crcpass != 1: - if self.monlevel > 2: - print("crc error in GPSRawCA/QZSRawCA " + - "{:6d}\t{:2d}\t{:1d}\t{:2d}". - format(int(self.tow), prn, crcpass, src)) - return -1 - - fh_ = self.fh_gpslnav if sys == uGNSS.GPS else self.fh_qzslnav + fh_ = self.fh_gpslnav blen = (300+7)//8 - fh_.write("{:4d}\t{:6d}\t{:3d}\t{:1d}\t{:3d}\t". - format(self.week, int(self.tow), prn, src, blen)) + fh_.write("{:4d}\t{:6d}\t{:3d}\t{:3d}\t". + format(self.week, int(self.tow), prn, blen)) - msg = bytearray(40) - for i in range(10): - d = st.unpack_from('L', msg, i*4, d) - k += 4 + for i in range(blen): + fh_.write("{:02x}".format(msg[i])) fh_.write("\n") - msg = bytes(msg) + + if self.flg_rnxnav: eph = self.rn.decode_gps_lnav(self.week, self.tow, sat, msg) if eph is not None: self.re.rnx_nav_body(eph, self.fh_rnxnav) elif blk_num in (4018, 4019, 4067, 4068): # GPSRawL2C/L5, QZSRawL2C/L5 + if (blk_num in [4018, 4019] and uGNSS.GPS not in self.sys_t) or \ + (blk_num in [4067, 4068] and uGNSS.QZS not in self.sys_t): + return 0 + sys, prn = self.decode_head(buff, k) k += 7 crcpass, cnt, src, freq, ch = st.unpack_from(' 2: + print("crc error in GPSRawL2C/L5, QZSRawL2C/L5 " + + "{:6d}\t{:2d}\t{:1d}\t{:1d}\t{:2d}". + format(int(self.tow), prn, crcpass, cnt, src)) + return -1 + + type_ = 3 if blk_num in [4018, 4067] else 4 + msg = bytearray(40) + for i in range(10): + st.pack_into('>L', msg, i*4, st.unpack_from(' 2: - print("crc error in GPSRawL2C/L5, QZSRawL2C/L5 " + - "{:6d}\t{:2d}\t{:1d}\t{:1d}\t{:2d}". - format(int(self.tow), prn, crcpass, cnt, src)) - return -1 - - fh_ = self.fh_gpscnav if sys == uGNSS.GPS else self.fh_qzscnav + fh_ = self.fh_gpscnav blen = (300+7)//8 fh_.write("{:4d}\t{:6d}\t{:3d}\t{:1d}\t{:3d}\t". - format(self.week, int(self.tow), prn, src, blen)) - msg = bytearray(40) - for i in range(10): - d = st.unpack_from('L', msg, i*4, d) - k += 4 + format(self.week, int(self.tow), prn, type_, blen)) + for i in range(blen): + fh_.write("{:02x}".format(msg[i])) fh_.write("\n") + if self.flg_rnxnav: sat = prn2sat(sys, prn) eph = self.rn.decode_gps_cnav( self.week, self.tow, sat, msg) @@ -567,27 +596,31 @@ def decode(self, buff, len_, sys=[], prn=[]): k += 7 crcpass, cnt, src, freq, ch = st.unpack_from(' 2: + print("crc error in GEORawL1/5 " + + "{:6d}\t{:2d}\t{:1d}\t{:1d}\t{:2d}". + format(int(self.tow), prn, crcpass, cnt, src)) + return -1 + + msg = bytearray(32) + for i in range(8): + d = st.unpack_from('L', msg, i*4, d) + k += 4 + if self.flg_sbas: itype = src-24 # 0:L1, 1:L5 - if crcpass != 1: - if self.monlevel > 2: - print("crc error in GEORawL1/5 " + - "{:6d}\t{:2d}\t{:1d}\t{:1d}\t{:2d}". - format(int(self.tow), prn, crcpass, cnt, src)) - return -1 - if prn < 120 or prn > 158: - return 0 - + if self.prn_ref > 0 and prn != self.prn_ref: return 0 - - msg = bytearray(32) - for i in range(8): - d = st.unpack_from('L', msg, i*4, d) - k += 4 - - self.output_sbas(prn, msg, self.fh_sbas, itype) + + self.output_sbas(prn, msg, self.fh_sbas, itype) + + if self.flg_rnxnav: + if prn < 120 or prn > 158: + return 0 sat = prn2sat(uGNSS.SBS, prn) seph = None @@ -597,47 +630,82 @@ def decode(self, buff, len_, sys=[], prn=[]): self.re.rnx_snav_body(seph, self.fh_rnxnav) elif blk_num == 4022: # GalRawFNAV + if uGNSS.GAL not in self.sys_t: + return 0 + sys, prn = self.decode_head(buff, k) sat = prn2sat(sys, prn) k += 7 crcpass, cnt, src, freq, ch = st.unpack_from(' 2: + print("crc error in GALRawFNAV " + + "{:6d}\t{:2d}\t{:1d}\t{:1d}\t{:2d}". + format(int(self.tow), prn, crcpass, cnt, + src & 0x1f)) + return -1 + + msg = bytearray(32) + for i in range(8): + d = st.unpack_from('L', msg, i*4, d) + k += 4 + if self.flg_galfnav: if src & 0x1f == 20: # E5a type_ = 1 else: return -1 - if crcpass != 1: - if self.monlevel > 2: - print("crc error in GALRawFNAV " + - "{:6d}\t{:2d}\t{:1d}\t{:1d}\t{:2d}". - format(int(self.tow), prn, crcpass, cnt, - src & 0x1f)) - return -1 - + blen = 32 self.fh_galfnav.write("{:4d}\t{:6d}\t{:3d}\t{:1d}\t{:3d}\t". format(self.week, int(self.tow), prn, - type_, 32)) - msg = bytearray(32) - for i in range(8): - d = st.unpack_from('L', msg, i*4, d) - k += 4 - + type_, blen)) + for i in range(blen): + self.fh_galfnav.write("{:02x}".format(msg[i])) self.fh_galfnav.write("\n") + if self.flg_rnxnav: eph = self.rn.decode_gal_fnav(self.week, self.tow, sat, 1, msg) if eph is not None: self.re.rnx_nav_body(eph, self.fh_rnxnav) elif blk_num == 4023: # GALRawINAV + if uGNSS.GAL not in self.sys_t: + return 0 + sys, prn = self.decode_head(buff, k) sat = prn2sat(sys, prn) k += 7 crcpass, cnt, src, freq, ch = st.unpack_from(' 1: + print("crc error in GALRawINAV " + + "{:6d}\t{:2d}\t{:1d}\t{:1d}\t{:2d}". + format(int(self.tow), prn, crcpass, cnt, + src & 0x1f)) + return -1 + + msg = bytearray(32) + for i in range(8): + d = st.unpack_from('L', msg, i*4, d) + k += 4 + + # GALRawINAV is missing tail bit (6) of even page + # add 6 bits offset for odd page + msg_ = bytearray(30) + msg_[0:15] = msg[0:15] # even page + k = 114 + for i in range(15): + d = bs.unpack_from('u8', bytes(msg), k)[0] + bs.pack_into('u8', msg_, 120+i*8, d) + k += 8 + if self.flg_galinav: if src & 0x1f == 17: # E1B type_ = 0 @@ -648,60 +716,43 @@ def decode(self, buff, len_, sys=[], prn=[]): print(f"unknown src: {src}") return -1 - if crcpass != 1: - if self.monlevel > 1: - print("crc error in GALRawINAV " + - "{:6d}\t{:2d}\t{:1d}\t{:1d}\t{:2d}". - format(int(self.tow), prn, crcpass, cnt, - src & 0x1f)) - return -1 - + blen = 30 self.fh_galinav.write("{:4d}\t{:6d}\t{:3d}\t{:1d}\t{:3d}\t". format(self.week, int(self.tow), prn, - type_, 30)) - msg = bytearray(32) - for i in range(8): - d = st.unpack_from('L', msg, i*4, d) - k += 4 + type_, blen)) - # GALRawINAV is missing tail bit (6) of even page - # add 6 bits offset for odd page - msg_ = bytearray(30) - msg_[0:15] = msg[0:15] # even page - k = 114 - for i in range(15): - d = bs.unpack_from('u8', bytes(msg), k)[0] - bs.pack_into('u8', msg_, 120+i*8, d) - k += 8 - - for i in range(30): + for i in range(blen): self.fh_galinav.write("{:02x}".format(msg_[i])) self.fh_galinav.write("\n") + if self.flg_rnxnav: eph = self.rn.decode_gal_inav(self.week, self.tow, sat, 2, msg_) if self.mode_galinav == 0 and eph is not None: self.re.rnx_nav_body(eph, self.fh_rnxnav) elif blk_num == 4024: # GALRawCNAV + if uGNSS.GAL not in self.sys_t: + return 0 + sys, prn = self.decode_head(buff, k) k += 7 crcpass, cnt, src, freq, ch = st.unpack_from(' 2: + print("crc error in GALRawCNAV " + + "{:6d}\t{:2d}\t{:1d}\t{:1d}\t{:2d}". + format(int(self.tow), prn, crcpass, cnt, src)) + return -1 + if self.flg_gale6: if src & 0x1f == 19: type_ = 6 else: return -1 - if crcpass != 1: - if self.monlevel > 2: - print("crc error in GALRawCNAV " + - "{:6d}\t{:2d}\t{:1d}\t{:1d}\t{:2d}". - format(int(self.tow), prn, crcpass, cnt, src)) - return -1 - blen = (492+7)//8 self.fh_gale6.write("{:4d}\t{:6d}\t{:3d}\t{:1d}\t{:3d}\t". format(self.week, int(self.tow), prn, @@ -713,25 +764,29 @@ def decode(self, buff, len_, sys=[], prn=[]): self.fh_gale6.write("\n") elif blk_num == 4026: # GLORawCA + if uGNSS.GLO not in self.sys_t: + return 0 + sys, prn = self.decode_head(buff, k) k += 7 crcpass, cnt, src, freq, ch = st.unpack_from(' 2: - print("crc error in GLORawCA " + - "{:6d}\t{:2d}\t{:1d}\t{:1d}\t{:2d}". - format(int(self.tow), prn, crcpass, cnt, src)) - return -1 - msg = bytearray(12) - for i in range(3): - d = st.unpack_from('L', msg, i*4, d) - k += 4 + if crcpass != 1: + if self.monlevel > 2: + print("crc error in GLORawCA " + + "{:6d}\t{:2d}\t{:1d}\t{:1d}\t{:2d}". + format(int(self.tow), prn, crcpass, cnt, src)) + return -1 + + msg = bytearray(12) + for i in range(3): + d = st.unpack_from('L', msg, i*4, d) + k += 4 + if self.flg_rnxnav: sat = prn2sat(sys, prn) geph = self.rn.decode_glo_fdma(self.week, self.tow, sat, msg, freq) @@ -745,18 +800,22 @@ def decode(self, buff, len_, sys=[], prn=[]): self.re.rnx_obs_body(obs, self.fh_rnxobs) elif blk_num == 4047: # BDSRaw + if uGNSS.BDS not in self.sys_t: + return 0 sys, prn = self.decode_head(buff, k) k += 7 # src 28: B1I (2I), 29: B2I (7I), 30: B3I (6I) crcpass, cnt, src, _, ch = st.unpack_from(' 2: - print("crc error in BDSRaw " + - "{:6d}\t{:2d}\t{:1d}\t{:1d}\t{:2d}". - format(int(self.tow), prn, crcpass, cnt, src)) - return -1 + + if crcpass != 1: + if self.monlevel > 2: + print("crc error in BDSRaw " + + "{:6d}\t{:2d}\t{:1d}\t{:1d}\t{:2d}". + format(int(self.tow), prn, crcpass, cnt, src)) + return -1 + + if self.flg_rnxnav and src == 28: # only D1 is supported msg = bytearray(40) for i in range(10): @@ -777,7 +836,7 @@ def decode(self, buff, len_, sys=[], prn=[]): if eph is not None: self.re.rnx_nav_body(eph, self.fh_rnxnav) - elif blk_num == 4069: # QZSRawL6 + elif blk_num in [4069, 4270, 4271]: # QZSRawL6, QZSRawL6D, QZSRawL6E sys, prn = self.decode_head(buff, k) k += 7 parity, rscnt, src, res, ch = st.unpack_from('>16)) + else: + self.fh_qzsl6.write("{:08x}".format(d)) k += 4 self.fh_qzsl6.write("\n") elif blk_num == 4093: # NAVICRaw + if uGNSS.IRN not in self.sys_t: + return 0 sys, prn = self.decode_head(buff, k) k += 7 crcpass, cnt, src, _, ch = st.unpack_from(' 2: - print("crc error in NAVICRaw " + - "{:6d}\t{:2d}\t{:1d}\t{:1d}\t{:2d}". - format(int(self.tow), prn, crcpass, cnt, src)) - return -1 + if crcpass != 1: + if self.monlevel > 2: + print("crc error in NAVICRaw " + + "{:6d}\t{:2d}\t{:1d}\t{:1d}\t{:2d}". + format(int(self.tow), prn, crcpass, cnt, src)) + return -1 + + if self.flg_rnxnav: msg = bytearray(40) for i in range(10): d = st.unpack_from(' 2: + print("crc error in BDSRawB1C " + + "{:6d}\t{:2d}\t{:1d}\t{:1d}\t{:2d}". + format(int(self.tow), prn, (crcsf2<<1)|crcsf3, + cnt, src)) + return -1 + + if self.flg_rnxnav: + #self.fh_bdsb1c.write("{:4d}\t{:6d}\t{:3d}\t{:3d}\t". + # format(self.week, int(self.tow), prn,225)) # 1800 deinterleaved symbols of a BeiDou B1C # (B-CNAV1) navigation frame # 24 unused bits in NAVBits[56] @@ -841,8 +916,8 @@ def decode(self, buff, len_, sys=[], prn=[]): d = st.unpack_from('L', v, i*4, d) - # self.fh_bdsb1c.write("{:08x}".format(d)) - # self.fh_bdsb1c.write("\n") + #self.fh_bdsb1c.write("{:08x}".format(d)) + #self.fh_bdsb1c.write("\n") v = bytes(v) prn_ = bs.unpack_from('u6', v, 0)[0] if prn != prn_: @@ -862,18 +937,21 @@ def decode(self, buff, len_, sys=[], prn=[]): self.re.rnx_nav_body(eph, self.fh_rnxnav) elif blk_num == 4219: # BDSRawB2a + if uGNSS.BDS not in self.sys_t: + return 0 sys, prn = self.decode_head(buff, k) k += 7 crcpass, cnt, src, _, ch = st.unpack_from(' 2: - print("crc error in BDSRawB2a " + - "{:6d}\t{:2d}\t{:1d}\t{:1d}\t{:2d}". - format(int(self.tow), prn, crcpass, cnt, src)) - return -1 - + + if crcpass != 1: + if self.monlevel > 2: + print("crc error in BDSRawB2a " + + "{:6d}\t{:2d}\t{:1d}\t{:2d}". + format(int(self.tow), prn, crcpass, cnt)) + return -1 + + if self.flg_rnxnav: msg = bytearray(40) for i in range(10): d = st.unpack_from(' 2: + print("crc error in GPSRawL1C, QZSRawL1C " + + "{:6d}\t{:2d}\t{:1d}\t{:1d}\t{:2d}". + format(int(self.tow), prn, + crcsf2, crcsf3, src)) + return -1 + + blen = (1800+7)//8 + msg = bytearray(228) + for i in range(57): + d = st.unpack_from('L', msg, i*4, d) + k += 4 + if (sys == uGNSS.GPS and self.flg_gpscnav2) or \ (sys == uGNSS.QZS and self.flg_qzscnav2): - if crcsf2 != 1 or crcsf3 != 1: - if self.monlevel > 2: - print("crc error in GPSRawL1C, QZSRawL1C " + - "{:6d}\t{:2d}\t{:1d}\t{:1d}\t{:2d}". - format(int(self.tow), prn, - crcsf2, crcsf3, src)) - return -1 - - fh_ = self.fh_gpscnav2 if sys == uGNSS.GPS \ - else self.fh_qzscnav2 - - blen = (1800+7)//8 + fh_ = self.fh_gpscnav2 fh_.write("{:4d}\t{:6d}\t{:3d}\t{:1d}\t{:3d}\t". format(self.week, int(self.tow), prn, src, blen)) - msg = bytearray(228) - for i in range(57): - d = st.unpack_from('L', msg, i*4, d) - k += 4 + for i in range(blen): + fh_.write("{:02x}".format(msg[i])) fh_.write("\n") + if self.flg_rnxnav: sat = prn2sat(sys, prn) eph = self.rn.decode_gps_cnav2(self.week, self.tow, sat, msg) if eph is not None: self.re.rnx_nav_body(eph, self.fh_rnxnav) elif blk_num in (4228, 4246): # QZSRawL1S, QZSRawL5S + if uGNSS.QZS not in self.sys_t: + return 0 + src_t = {24: 0, 25: 1, 33: 2, 39: 3} # L1C/A, L5, L1S, L5S sys, prn = self.decode_head(buff, k) k += 7 crc, cnt, src, freq, ch = st.unpack_from(' 2: - print("crc error in QZSRawL1S, QZSRawL5S " + - "{:6d}\t{:2d}\t{:1d}\t{:1d}\t{:2d}". - format(int(self.tow), prn, crc, cnt, src)) - return -1 + + if crc != 1: + if self.monlevel > 2: + print("crc error in QZSRawL1S, QZSRawL5S " + + "{:6d}\t{:2d}\t{:1d}\t{:1d}\t{:2d}". + format(int(self.tow), prn, crc, cnt, src)) + return -1 + if self.flg_qzsl1s or self.flg_qzsl5s: if src not in src_t.keys(): if self.monlevel > 0: print("src not recognized in QZSRawL1S/QZSRawL5S " + @@ -949,64 +1036,67 @@ def decode(self, buff, len_, sys=[], prn=[]): self.output_sbas(prn, msg, self.fh_sbas, itype) elif blk_num == 4242: # BDSRawB2b + if uGNSS.BDS not in self.sys_t: + return 0 + sys, prn = self.decode_head(buff, k) k += 7 crcpass, _, src, _, ch = st.unpack_from('= 59 else False - - if crcpass != 1: - if self.monlevel > 2: - print("crc error in BDSRawB2b " + - "{:6d}\t{:2d}\t{:1d}\t{:2d}". - format(int(self.tow), prn, crcpass, src)) - return -1 - - if flg_bdsppp: - self.fh_bdsb2b.write("{:4d}\t{:6d}\t{:3d}\t{:1d}\t{:3d}\t". - format(self.week, int(self.tow), prn, - src, 64)) - # 984 symbols of a BeiDou B2b navigation frame - # 8 unused bits in NAVBits[30] - msg = bytearray(64) - for i in range(16): - d = st.unpack_from('L', msg, i*4, d) - k += 4 + if crcpass != 1: + if self.monlevel > 2: + print("crc error in BDSRawB2b " + + "{:6d}\t{:2d}\t{:1d}\t{:2d}". + format(int(self.tow), prn, crcpass, src)) + return -1 - if flg_bdsppp: - if i == 0: - self.fh_bdsb2b.write("{:05x}".format(d & 0xfffff)) - elif i == 15: - self.fh_bdsb2b.write( - "{:05x}{:06x}".format((d >> 12) & 0xffffc, 0)) - else: - self.fh_bdsb2b.write("{:08x}".format(d)) - - if flg_bdsppp: - self.fh_bdsb2b.write("\n") - else: # B2b-CNAV3 - sat = prn2sat(sys, prn) - eph = self.rn.decode_bds_b2b(self.week, self.tow, sat, msg) - if eph is not None: - self.re.rnx_nav_body(eph, self.fh_rnxnav) + # 984 symbols of a BeiDou B2b navigation frame + # 8 unused bits in NAVBits[30] + msg = bytearray(64) + for i in range(16): + d = st.unpack_from('L', msg, i*4, d) + k += 4 + + + if self.flg_bdsb2b and prn >= 59: + self.fh_bdsb2b.write("{:4d}\t{:6d}\t{:3d}\t{:3d}\t". + format(self.week, int(self.tow), prn, 64)) + + blen = 64 + for i in range(blen): # skip MSB 12bit + if i==0: + continue + elif i==1: + self.fh_bdsb2b.write("{:01x}".format(msg[i]&0xf)) + else: + self.fh_bdsb2b.write("{:02x}".format(msg[i])) + self.fh_bdsb2b.write("\n") + + if self.flg_rnxnav and prn<59: + sat = prn2sat(sys, prn) + eph = self.rn.decode_bds_b2b(self.week, self.tow, sat, msg) + if eph is not None: + self.re.rnx_nav_body(eph, self.fh_rnxnav) elif blk_num == 4262: # NAVICRawL1 + if uGNSS.IRN not in self.sys_t: + return 0 sys, prn = self.decode_head(buff, k, sysref=uGNSS.IRN) k += 7 crc_sf2, crc_sf3, src, _, ch = st.unpack_from(' 2: - print("crc error in NAVICRawL1 " + - "{:6d}\t{:2d}\t{:1d}\t{:1d}\t{:2d}". - format(int(self.tow), prn, crc_sf2, crc_sf3, - src)) - return -1 + if crc_sf2 != 1 or crc_sf3 != 1: + if self.monlevel > 2: + print("crc error in NAVICRawL1 " + + "{:6d}\t{:2d}\t{:1d}\t{:1d}\t{:2d}". + format(int(self.tow), prn, crc_sf2, crc_sf3, + src)) + return -1 + + if self.flg_rnxnav: msg = bytearray(228) for i in range(57): d = st.unpack_from(' 1: @@ -432,7 +438,7 @@ def decode(f, opt, args): bdir, fname = os.path.split(f) - prefix = fname.removesuffix('.ubx')[-4:]+'_' + prefix = fname.removesuffix('.ubx') prefix = str(Path(bdir) / prefix) if bdir else prefix ubxdec = ubx(opt, prefix=prefix, gnss_t=args.gnss) ubxdec.monlevel = 1 diff --git a/samples/config.yml b/samples/config.yml new file mode 100644 index 0000000..556f1c0 --- /dev/null +++ b/samples/config.yml @@ -0,0 +1,22 @@ +# config.yml + +elmin: 10.0 # minimum elevation [deg] +atxfile: '../data/antex/igs20.atx' # ANTEX file + +nav: + monlevel: 1 + pmode: 1 # positioning mode 0:static, 1:kinematic + csmooth: True # carrier smoothing is enabled + + ephopt: 2 # SSR-APC + parmode: 2 # Partial AR 1: normal, 2: PAR + trop_opt: 1 + iono_opt: 1 + phw_opt: 1 + cmooth: False + trop_model: 0 + iono_model: 0 + +cs: + monlevel: 1 + diff --git a/samples/config_ppp.yml b/samples/config_ppp.yml new file mode 100644 index 0000000..783f8ae --- /dev/null +++ b/samples/config_ppp.yml @@ -0,0 +1,46 @@ + # config.yml for PPP (static) + +elmin: 10.0 # minimum elevation [deg] +atxfile: '../data/antex/igs20.atx' # ANTEX file + +nav: + monlevel: 1 + pmode: 0 # positioning mode 0:static, 1:kinematic + csmooth: False # carrier smoothing is enabled + ephopt: 2 # SSR-APC + trop_opt: 1 + iono_opt: 1 + phw_opt: 1 + trop_model: 'SAAST' + iono_model: 'KLOBUCHAR' + + eratio: [50.0, 50.0] + err: [0, 0.007, 0.0035] + # excl_sat: ['G02'] + + # initial value of covariance + sig_p0: 100.0 # [m] + sig_v0: 1.0 # [m/s] + sig_ztd0: 0.1 # [m] + sig_ion0: 10.0 # [m] + sig_n0: 30.0 # [cyc] + + # Process noise sigma + # sig_qp: 100.0 + sig_qp: 0.01 # [m/sqrt(s)] + sig_qv: 1.0 # [m] + sig_qztd: 0.0008 # [m] + sig_qion: 10.0 # [m] + sig_qb: 0.0001 # [m] + + # AR parameters + armode: 0 # 0:float-ppp,1:continuous,2:instantaneous,3:fix-and-hold + thresar: 3.0 # AR acceptance threshold + elmaskar: 15.0 + + parmode: 1 # Partial AR 1: normal, 2: PAR + par_P0: 0.995 # probability of sussefull AR + +cs: + monlevel: 1 + diff --git a/samples/config_ppprtk.yml b/samples/config_ppprtk.yml new file mode 100644 index 0000000..85152c4 --- /dev/null +++ b/samples/config_ppprtk.yml @@ -0,0 +1,47 @@ + # config.yml for PPP-RTK (static) + +elmin: 10.0 # minimum elevation [deg] +atxfile: '../data/antex/igs20.atx' # ANTEX file +griddef: '../data/clas_grid.def' # grid defintion for local correction + +nav: + monlevel: 1 + pmode: 0 # positioning mode 0:static, 1:kinematic + csmooth: False # carrier smoothing is enabled + ephopt: 2 # SSR-APC + trop_opt: 2 + iono_opt: 2 + phw_opt: 2 + trop_model: 'SAAST' + iono_model: 'KLOBUCHAR' + + eratio: [50.0, 50.0] + err: [0, 0.007, 0.0035] + + # initial value of covariance + sig_p0: 100.0 # [m] + sig_v0: 1.0 # [m/s] + sig_ztd0: 0.1 # [m] + sig_ion0: 10.0 # [m] + sig_n0: 30.0 # [cyc] + + # Process noise sigma + # sig_qp: 100.0 + sig_qp: 0.01 # [m/sqrt(s)] + sig_qv: 1.0 # [m] + sig_qztd: 0.0008 # [m] + sig_qion: 10.0 # [m] + sig_qb: 0.0001 # [m] + + # AR parameters + armode: 3 # 0:float-ppp,1:continuous,2:instantaneous,3:fix-and-hold + thresar: 2.0 # AR acceptance threshold + elmaskar: 20.0 + elmin: 10.0 + + parmode: 1 # Partial AR 1: normal, 2: PAR + par_P0: 0.995 # probability of sussefull AR + +cs: + monlevel: 1 + diff --git a/samples/igs_download.py b/samples/igs_download.py index 6bfe385..83b9bc4 100644 --- a/samples/igs_download.py +++ b/samples/igs_download.py @@ -12,7 +12,7 @@ # Configuration # LOCAL_BASE_DIR = "../data" # Data directory -FILE_LIST = "../data/igs_files.txt" # File containing destination & FTP URLs +DEFAULT_FILE_LIST = "../data/igs_files.txt" # File containing destination & FTP URLs def gps_to_datetime(gps_week, gps_day_of_week): @@ -74,14 +74,14 @@ def connect_ftp(host): return ftp -def read_file_list(): +def read_file_list(file_list_path=DEFAULT_FILE_LIST): """Read the list of destination folders and FTP URLs from a file.""" - if not os.path.exists(FILE_LIST): - print(f"Error: {FILE_LIST} not found!") + if not os.path.exists(file_list_path): + print(f"Error: {file_list_path} not found!") return [] entries = [] - with open(FILE_LIST, "r") as f: + with open(file_list_path, "r") as f: for line in f: parts = line.strip().split(maxsplit=1) if line.startswith('#') or len(parts) == 0: @@ -161,9 +161,9 @@ def download_http(url, filename, local_filename): print(f"Warning: Failed to download {filename} from {url} - {e}") -def download_files(): +def download_files(file_list_path=DEFAULT_FILE_LIST): """Download files from the provided list of FTP URLs.""" - entries = read_file_list() + entries = read_file_list(file_list_path) if not entries: print("No valid entries found. Exiting.") return @@ -244,4 +244,14 @@ def download_files(): if __name__ == "__main__": - download_files() + import argparse + parser = argparse.ArgumentParser(description="Download IGS files from a list.") + parser.add_argument( + "file_list", + nargs="?", + default=DEFAULT_FILE_LIST, + help=f"Path to the file containing destination folders & FTP URLs (default: {DEFAULT_FILE_LIST})" + ) + args = parser.parse_args() + + download_files(args.file_list) diff --git a/samples/test_dgps.py b/samples/test_dgps.py index a71206b..ea2b7dd 100644 --- a/samples/test_dgps.py +++ b/samples/test_dgps.py @@ -2,314 +2,121 @@ static test for DGPS (QZSS SLAS) """ from binascii import unhexlify -from copy import deepcopy -import matplotlib.pyplot as plt -import matplotlib.dates as md import numpy as np -from sys import exit as sys_exit -from sys import stdout - -import cssrlib.gnss as gn -from cssrlib.gnss import ecef2pos, Nav -from cssrlib.gnss import time2gpst, time2doy, time2str, timediff, epoch2time -from cssrlib.gnss import rSigRnx, sys2str, uIonoModel -from cssrlib.peph import atxdec, searchpcv +from cssrlib.gnss import Nav, time2gpst, time2doy, epoch2time +from cssrlib.gnss import rSigRnx, uIonoModel, load_config from cssrlib.pntpos import stdpos from cssrlib.dgps import dgpsDec from cssrlib.rinex import rnxdec +from cssrlib.utils import process + + +def decode_msg(v, tow, prn_ref): + """ find valid correction message """ + + vi = v[(v['tow'] == tow) & (v['prn'] == prn_ref) + & ((v['type'] >= 47) & (v['type'] <= 51))] + if len(vi) > 0: + msg = unhexlify(vi['nav'][0]) + else: + msg = None + + return msg +config = load_config('config.yml') + # Select test case # dataset = 0 +nep = 3600 +ttl = 'test_dgps' # title for case + +navfile = None # Start epoch and number of epochs # -if dataset == 0: # MSAS, L1 SBAS - ep = [2023, 8, 11, 21, 0, 0] - # navfile = '../data/doy2023-308/308c_rnx.nav' - navfile = '../data/brdc/BRD400DLR_S_20232230000_01D_MN.rnx' - obsfile = '../data/doy2023-223/SEPT223Y.23O' # PolaRX5 - file_sbas = '../data/doy2023-223/223v_sbas.txt' - xyz_ref = [-3962108.6726, 3381309.4719, 3668678.6264] - prn_ref = 189 # satellite PRN for SLAS - sbas_type = 0 # L1: 0, L5: 1 - nf = 1 +if dataset == 0: # QZSS SLAS + ep = [2025, 8, 21, 7, 0, 0] + xyz_ref = [-3962108.6836, 3381309.5672, 3668678.6720] # Kamakura -elif dataset == 1: # MSAS, L1 SBAS - ep = [2025, 2, 15, 12, 0, 0] - navfile = '../data/doy2025-046/046r_rnx.nav' - obsfile = '../data/doy2025-046/046r_rnx.obs' # PolaRX5 - file_sbas = '../data/doy2025-046/046m_sbas.txt' - xyz_ref = [-3962108.6726, 3381309.4719, 3668678.6264] - prn_ref = 189 # satellite PRN for SBAS correction + prn_ref = 199 # satellite PRN for SBAS correction sbas_type = 0 # L1: 0, L5: 1 nf = 1 -elif dataset == 2: # SouthPAN L5 - print("ERROR: datset not yet available!") - sys_exit(1) - - time = epoch2time(ep) year = ep[0] doy = int(time2doy(time)) +ses = chr(ord('a')+ep[3]) -nep = 3600 -# nep = 360 - -pos_ref = ecef2pos(xyz_ref) +if navfile is None: + navfile = f'../data/doy{year}-{doy:03d}/{doy:03d}{ses}_rnx.nav' + obsfile = f'../data/doy{year}-{doy:03d}/{doy:03d}{ses}_rnx.obs' + file_sbas = f'../data/doy{year}-{doy:03d}/{doy:03d}{ses}_sbas.txt' -# Define signals to be processed -# -sigs = [] - -if sbas_type == 0: # single frequency L1 SBAS - nf = 1 - gnss = "GJ" - if 'G' in gnss: - sigs.extend([rSigRnx("GC1C"), rSigRnx("GL1C"), rSigRnx("GS1C")]) - if 'J' in gnss: - sigs.extend([rSigRnx("JC1C"), rSigRnx("JL1C"), rSigRnx("JS1C")]) +if sbas_type == 0: # DGNSS + sig_t = {'G': ['1C'], 'J': ['1C']} else: # dual frequency - nf = 2 - gnss = "GE" - if 'G' in gnss: - sigs.extend([rSigRnx("GC1C"), rSigRnx("GC5Q"), - rSigRnx("GL1C"), rSigRnx("GL5Q"), - rSigRnx("GS1C"), rSigRnx("GS5Q")]) - if 'E' in gnss: - sigs.extend([rSigRnx("EC1C"), rSigRnx("EC5Q"), - rSigRnx("EL1C"), rSigRnx("EL5Q"), - rSigRnx("ES1C"), rSigRnx("ES5Q")]) + sig_t = {'G': ['1C', '5Q'], 'E': ['1C', '5Q']} rnx = rnxdec() -rnx.setSignals(sigs) - nav = Nav() +proc = process(nav, rnx, config, nep=nep, xyz_ref=xyz_ref) -# Positioning mode -# 0:static, 1:kinematic +# Define signals to be processed # -nav.pmode = 0 +sigs, nav.nf = proc.init_sig(sig_t) +rnx.setSignals(sigs) # Decode RINEX NAV data # nav = rnx.decode_nav(navfile, nav) -cs = dgpsDec('test_dgps_cs.log') -cs.monlevel = 0 - -# Load ANTEX data for satellites and stations -# -atx = atxdec() -atx.readpcv('../data/antex/igs20.atx') - -# Initialize data structures for results -# -t = np.zeros(nep) -enu = np.ones((nep, 3))*np.nan -sol = np.zeros((nep, 4)) -ztd = np.zeros((nep, 1)) -smode = np.zeros(nep, dtype=int) - -# Logging level -# -nav.monlevel = 1 # TODO: enabled for testing! +cs = dgpsDec(f'{ttl}_cs.log') +cs.monlevel = config['cs']['monlevel'] # Load RINEX OBS file header # -if rnx.decode_obsh(obsfile) >= 0: - - # Auto-substitute signals - # - rnx.autoSubstituteSignals() - - # Initialize position - # - std = stdpos(nav, rnx.pos, 'test_dgps.log', trop_opt=2, iono_opt=2) - nav.elmin = np.deg2rad(5.0) - - std.ionoModel = uIonoModel.SBAS +rnx.decode_obsh(obsfile) - nav.nf = nf - - # Get equipment information - # - nav.fout.write("FileName: {}\n".format(obsfile)) - nav.fout.write("Start : {}\n".format(time2str(rnx.ts))) - if rnx.te is not None: - nav.fout.write("End : {}\n".format(time2str(rnx.te))) - nav.fout.write("Receiver: {}\n".format(rnx.rcv)) - nav.fout.write("Antenna : {}\n".format(rnx.ant)) - nav.fout.write("\n") +std = stdpos(nav, rnx.pos, f'{ttl}.log', trop_opt=2, iono_opt=2) +std.ionoModel = uIonoModel.SBAS - if 'UNKNOWN' in rnx.ant or rnx.ant.strip() == "": - nav.fout.write("ERROR: missing antenna type in RINEX OBS header!\n") +proc.prepare_signal(obsfile) - # Set PCO/PCV information - # - nav.sat_ant = atx.pcvs - nav.rcv_ant = searchpcv(atx.pcvr, rnx.ant, rnx.ts) - if nav.rcv_ant is None: - nav.fout.write("ERROR: missing antenna type <{}> in ANTEX file!\n" - .format(rnx.ant)) - - # Print available signals - # - nav.fout.write("Available signals\n") - for sys, sigs in rnx.sig_map.items(): - txt = "{:7s} {}\n".format(sys2str(sys), ' '. - join([sig.str() for sig in sigs.values()])) - nav.fout.write(txt) - nav.fout.write("\n") - - nav.fout.write("Selected signals\n") - for sys, tmp in rnx.sig_tab.items(): - txt = "{:7s} ".format(sys2str(sys)) - for _, sigs in tmp.items(): - txt += "{} ".format(' '.join([sig.str() for sig in sigs])) - nav.fout.write(txt+"\n") - nav.fout.write("\n") - - # Skip epochs until start time - # +# Skip epochs until start time +# +obs = rnx.decode_obs() +while time > obs.t and obs.t.time != 0: obs = rnx.decode_obs() - while time > obs.t and obs.t.time != 0: - obs = rnx.decode_obs() - - dtype = [('wn', 'int'), ('tow', 'float'), ('prn', 'int'), - ('type', 'int'), ('len', 'int'), ('nav', 'S124')] - v = np.genfromtxt(file_sbas, dtype=dtype) - - # Loop over number of epoch from file start - # - for ne in range(nep): - - week, tow = time2gpst(obs.t) - cs.week = week - cs.tow0 = tow//86400*86400 - cs.time = obs.t - - # Set initial epoch - # - if ne == 0: - nav.t = deepcopy(obs.t) - t0 = deepcopy(obs.t) - t0.time = t0.time//30*30 - nav.time_p = t0 - - vi = v[(v['tow'] == tow) & (v['prn'] == prn_ref) - & (v['type'] == sbas_type)] - if len(vi) > 0: - buff = unhexlify(vi['nav'][0]) - cs.decode_cssr(buff, 0) - - # cs.check_validity(obs.t) - - # Call PPP module with PVS corrections - # - if (cs.lc[0].cstat & 0x7) == 0x7: - std.process(obs, cs=cs) - # std.process(obs) - - # Save output - # - t[ne] = timediff(nav.t, t0)/86400.0 - - sol = nav.xa[0:3] if nav.smode == 4 else nav.x[0:3] - enu[ne, :] = gn.ecef2enu(pos_ref, sol-xyz_ref) - # ztd[ne] = nav.xa[std.IT(nav.na)] \ - # if nav.smode == 4 else nav.x[std.IT(nav.na)] - smode[ne] = nav.smode +dtype = [('wn', 'int'), ('tow', 'float'), ('prn', 'int'), + ('type', 'int'), ('marker', 'S2'), ('nav', 'S124')] +v = np.genfromtxt(file_sbas, dtype=dtype) - nav.fout.write("{} {:14.4f} {:14.4f} {:14.4f} " - "ENU {:7.3f} {:7.3f} {:7.3f}, 2D {:6.3f}, mode {:1d}\n" - .format(time2str(obs.t), - sol[0], sol[1], sol[2], - enu[ne, 0], enu[ne, 1], enu[ne, 2], - np.sqrt(enu[ne, 0]**2+enu[ne, 1]**2), - smode[ne])) - - # Log to standard output - # - stdout.write('\r {} ENU {:7.3f} {:7.3f} {:7.3f}, ' - '2D {:6.3f}, mode {:1d}' - .format(time2str(obs.t), - enu[ne, 0], enu[ne, 1], enu[ne, 2], - np.sqrt(enu[ne, 0]**2+enu[ne, 1]**2), - smode[ne])) - - # Get new epoch, exit after last epoch - # - obs = rnx.decode_obs() - if obs.t.time == 0: - break - - # Send line-break to stdout - # - stdout.write('\n') - - # Close RINEX observation file - # - rnx.fobs.close() +# Loop over number of epoch from file start +# +for ne in range(nep): + _, tow = cs.set_time(obs.t) # set time for reference - # Close output file + # Set initial epoch # - if nav.fout is not None: - nav.fout.close() - -fig_type = 1 -ylim_h = 2.0 -ylim_v = 6.0 - -idx2 = np.where(smode == 2)[0] -idx1 = np.where(smode == 1)[0] -idx0 = np.where(smode == 0)[0] - -fig = plt.figure(figsize=[7, 9]) -fig.set_rasterized(True) - -fmt = '%H:%M' - -if fig_type == 1: - - lbl_t = ['East [m]', 'North [m]', 'Up [m]'] - - for k in range(3): - ylim = ylim_h if k < 2 else ylim_v - plt.subplot(3, 1, k+1) - plt.plot(t[idx0], enu[idx0, k], 'r.', label='none') - plt.plot(t[idx2], enu[idx2, k], 'y.', label='SBAS/DGPS') - plt.plot(t[idx1], enu[idx1, k], 'g.', label='standalone') - - plt.ylabel(lbl_t[k]) - plt.grid() - plt.ylim([-ylim, ylim]) - plt.gca().xaxis.set_major_formatter(md.DateFormatter(fmt)) - - plt.xlabel('Time [HH:MM]') - plt.legend() - -elif fig_type == 2: + if ne == 0: + proc.init_time(obs.t) - ax = fig.add_subplot(111) + msg = decode_msg(v, tow, prn_ref) + if msg is not None: + cs.decode_cssr(msg, 0) # decode DGPS correction - plt.plot(enu[idx0, 0], enu[idx0, 1], 'r.', label='none') - plt.plot(enu[idx2, 0], enu[idx2, 1], 'y.', label='SBAS/DGPS') - plt.plot(enu[idx1, 0], enu[idx1, 1], 'g.', label='standalone') + if (cs.lc[0].cstat & 0x7) == 0x7: # wait for mask/clock/orbit + std.process(obs, cs=cs) # standalone positioning - plt.xlabel('Easting [m]') - plt.ylabel('Northing [m]') - plt.grid() - plt.axis('equal') - plt.legend() - # ax.set(xlim=(-ylim, ylim), ylim=(-ylim, ylim)) + proc.save_output(obs.t, ne) # save output -plotFileFormat = 'eps' -plotFileName = '.'.join(('test_dgps', plotFileFormat)) + obs = rnx.decode_obs() # get new epoch + if obs.t.time == 0: # exit after last epoch + break -plt.savefig(plotFileName, format=plotFileFormat, bbox_inches='tight', dpi=300) -# plt.show() +proc.close() +proc.plot(ttl, fig_type=1, ylim=2, ylim_v=4) diff --git a/samples/test_integ.py b/samples/test_integ.py index 350b6ab..242ff45 100644 --- a/samples/test_integ.py +++ b/samples/test_integ.py @@ -1,20 +1,16 @@ """ - interopetabity test for SSR Integrity messages MT11,12,13 + interopetabity test for SSR Integrity messages MT11,12,13,MT54.* @author Rui Hirokawa """ import os import copy -from cssrlib.gnss import uGNSS, prn2sat +from cssrlib.gnss import uGNSS, prn2sat, gpst2time from cssrlib.rtcm import rtcm, rtcme, Integrity from random import randint, seed, sample from binascii import unhexlify -# Validity Period DFi065 -vp_tbl = [1, 2, 5, 10, 15, 30, 60, 120, 240, - 300, 600, 900, 1800, 3600, 7200, 10800] - def read_asc(file): b = bytearray() @@ -25,6 +21,12 @@ def read_asc(file): return b +def read_bin(file): + with open(file, 'rb') as fh: + return fh.read() + return None + + def gen_data(mt, sys_t, svid_t): """ generate random message data """ intr = Integrity() @@ -43,7 +45,7 @@ def gen_data(mt, sys_t, svid_t): # issue of GNSS satellite mask DFi010 intr.pid = randint(0, 4095) # provider id DFi027 (0-4095) intr.tow = randint(0, 604799)*1e-3 # tow - intr.vp = vp_tbl[randint(0, 15)] # validity period DFi065 (0-15) + intr.vp = intr.vp_tbl[randint(0, 15)] # validity period DFi065 (0-15) intr.uri = randint(0, 65535)*0.1 # update rate interval DFi067 intr.pidssr = randint(0, 65535) # SSR Provider ID DFi078 (0-65535) @@ -100,9 +102,14 @@ def write_rtcm(file_rtcm, msg_t, intr, nep=1): return msg[:k] -def decode_rtcm(msg, intr=None, nep=1, logfile=None, maxlen=1024): +def decode_rtcm(msg, intr=None, nep=1, logfile=None, maxlen=1024, + mt_skip=None, weekref = -1): cs = rtcm(foutname=logfile) cs.monlevel = 2 + cs.week = weekref + cs.time = gpst2time(cs.week, 0) + if mt_skip is not None: + cs.mt_skip = mt_skip k = 0 for ne in range(nep): @@ -136,7 +143,8 @@ def decode_rtcm(msg, intr=None, nep=1, logfile=None, maxlen=1024): return cs -def read_rtcm(file_rtcm, intr, nep=1, logfile=None): +def read_rtcm(file_rtcm, intr, nep=1, logfile=None, mt_skip=None): + """ read test script for SC-134 messages """ fc = open(file_rtcm, 'rb') if not fc: @@ -147,45 +155,72 @@ def read_rtcm(file_rtcm, intr, nep=1, logfile=None): maxlen = len(msg)-5 fc.close() - return decode_rtcm(msg, intr, nep, logfile, maxlen) + return decode_rtcm(msg, intr, nep, logfile, maxlen, mt_skip) if __name__ == "__main__": - file_rtcm = '../data/sample.rtcm' - file_log = '../data/sample.log' - - file_asc = '../data/sc134/MT05_DFi56=00.txt' - - nep = 1 - maxlen = 1024 - nsatmax = 10 - - seed_ = 1 - # parameters - # msg_t = [11, 12, 13] - msg_t = [11] - # msg_t = [12] - # msg_t = [13] - - # 0:GPS,1:GLO,2:GAL,3:BDS,4:QZS,5:IRN - sys_t = [uGNSS.GPS, uGNSS.GLO, uGNSS.GAL, uGNSS.QZS] - - # GNSS satellite mask DFi009 - prn_rng_t = {uGNSS.GPS: [1, 32], # Table 8.5-1 - uGNSS.GLO: [1, 27], # Table 8.5-3 - uGNSS.GAL: [1, 36], # Table 8.5-5 - uGNSS.BDS: [1, 63], - uGNSS.QZS: [193, 209], - uGNSS.IRN: [1, 14]} - - mt = msg_t[0] - - seed(seed_) - prn_t = gen_sat_list(sys_t, prn_rng_t) # generate random sat list - intr = gen_data(mt, sys_t, prn_t) # generate random message data - msg = write_rtcm(file_rtcm, msg_t, intr, nep) - cs = read_rtcm(file_rtcm, intr, nep, logfile=file_log) - - # msg = read_asc(file_asc) - # cs = rtcm(foutname=file_log) - # decode_rtcm(msg) + bdir = '../data/sc134/msg/' + flg_sim = False + + mt_skip = [] + #mt_skip = [1267] # work-around for SC134 SSR interop-test + + if flg_sim: # generate test data + file_rtcm = bdir+'test.rtcm' + file_log = bdir+'test.log' + nep = 1 + maxlen = 1024 + nsatmax = 10 + + seed_ = 1 + # parameters + # msg_t = [11, 12, 13] + msg_t = [11] + # msg_t = [12] + # msg_t = [13] + + # 0:GPS,1:GLO,2:GAL,3:BDS,4:QZS,5:IRN + sys_t = [uGNSS.GPS, uGNSS.GLO, uGNSS.GAL, uGNSS.QZS] + + # GNSS satellite mask DFi009 + prn_rng_t = {uGNSS.GPS: [1, 32], # Table 8.5-1 + uGNSS.GLO: [1, 27], # Table 8.5-3 + uGNSS.GAL: [1, 36], # Table 8.5-5 + uGNSS.BDS: [1, 63], + uGNSS.QZS: [193, 209], + uGNSS.IRN: [1, 14]} + + mt = msg_t[0] + + seed(seed_) + prn_t = gen_sat_list(sys_t, prn_rng_t) # generate random sat list + intr = gen_data(mt, sys_t, prn_t) # generate random message data + msg = write_rtcm(file_rtcm, msg_t, intr, nep) + cs = read_rtcm(file_rtcm, intr, nep, logfile=file_log) + + else: # decode using sample dataset (*.bin) + + icase = 2 + + if icase == 1: + + file_rtcm = ['MT54_9', 'MT54_10_DFi209=0', + 'MT54_10_DFi209=1', 'MT54_10_DFi209=2'] + weekref = 2403 + + #file_rtcm = ['RTCM134test_21012026'] + + # file_rtcm = ['sampledataMT03-04-05-06-07'] + + elif icase == 2: + file_rtcm = ['ssr/SSRTEST_20260206_CORR_v123', + 'ssr/ROVRMSG', + 'ssr/ROMAMSG'] + weekref = 2404 + + # msg = read_asc(file_asc) + for f in file_rtcm: + file_log = bdir+f+'.dlg' + msg = read_bin(bdir+f+'.bin') + decode_rtcm(msg, logfile=file_log, maxlen=len(msg),mt_skip=mt_skip, + weekref=weekref) diff --git a/samples/test_pppbds.py b/samples/test_pppbds.py index fd16b87..601bb85 100644 --- a/samples/test_pppbds.py +++ b/samples/test_pppbds.py @@ -2,90 +2,70 @@ static test for PPP (BeiDou PPP) """ from binascii import unhexlify -from copy import deepcopy -import matplotlib.pyplot as plt import numpy as np -from sys import exit as sys_exit -from sys import stdout -import cssrlib.gnss as gn -from cssrlib.gnss import ecef2pos, Nav -from cssrlib.gnss import time2gpst, time2doy, time2str, timediff, epoch2time -from cssrlib.gnss import rSigRnx -from cssrlib.gnss import sys2str -from cssrlib.peph import atxdec, searchpcv +from cssrlib.gnss import Nav, load_config, time2doy, epoch2time from cssrlib.cssr_bds import cssr_bds from cssrlib.pppssr import pppos from cssrlib.rinex import rnxdec -from cssrlib.plot import plot_enu +from cssrlib.utils import process + + +def decode_msg(v, tow, prn_ref): + """ find valid correction message """ + msg = None + vi = v[(v['tow'] == tow) & (v['prn'] == prn_ref)] + if len(vi) > 0: + msg = unhexlify(vi['nav'][0]) + return msg + + +config = load_config('config_ppp.yml') # Select test case # -dataset = 3 +nep = 900*4 +prn_ref = 59 # satellite PRN to receive BDS PPP collection +dataset = 1 +navfile = None +file_bds = None +ttl = 'test_pppbds' # Start epoch and number of epochs # if dataset == 0: - ep = [2023, 7, 8, 4, 0, 0] - xyz_ref = [-3962108.7007, 3381309.5532, 3668678.6648] - navfile = '../data/brdc/BRD400DLR_S_20231890000_01D_MN.rnx' - obsfile = '../data/doy2023-189/SEPT1890.23O' - file_bds = '../data/doy2023-189/bdsb2b_189e.txt' -elif dataset == 1: - ep = [2023, 8, 11, 21, 0, 0] - xyz_ref = [-3962108.7007, 3381309.5532, 3668678.6648] - navfile = '../data/brdc/BRD400DLR_S_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 - file_bds = '../data/doy2023-223/223v_bdsb2b.txt' -elif dataset == 2: ep = [2025, 2, 15, 17, 0, 0] xyz_ref = [-3962108.6836, 3381309.5672, 3668678.6720] - navfile = '../data/doy2025-046/046r_rnx.nav' - obsfile = '../data/doy2025-046/046r_rnx.obs' # SEPT MOSAIC-X5 - file_bds = '../data/doy2025-046/046r_bdsb2b.txt' -elif dataset == 3: +elif dataset == 1: ep = [2025, 8, 21, 7, 0, 0] - navfile = '../data/doy2025-233/233h_rnx.nav' - # navfile = '../data/brdc/BRD400DLR_S_20252330000_01D_MN.rnx' - obsfile = '../data/doy2025-233/233h_rnx.obs' # SEPT MOSAIC-X5 - file_bds = '../data/doy2025-233/233h_bdsb2b.txt' xyz_ref = [-3962108.6836, 3381309.5672, 3668678.6720] # Kamakura time = epoch2time(ep) year = ep[0] doy = int(time2doy(time)) +ses = chr(ord('a')+ep[3]) + +if navfile is None: + navfile = f'../data/doy{year}-{doy:03d}/{doy:03d}{ses}_rnx.nav' + obsfile = f'../data/doy{year}-{doy:03d}/{doy:03d}{ses}_rnx.obs' +if file_bds is None: + file_bds = f'../data/doy{year}-{doy:03d}/{doy:03d}{ses}_bdsb2b.txt' -nep = 900*4 dtype = [('wn', 'int'), ('tow', 'int'), ('prn', 'int'), ('type', 'int'), ('len', 'int'), ('nav', 'S124')] v = np.genfromtxt(file_bds, dtype=dtype) -prn_ref = 59 # satellite PRN to receive BDS PPP collection +rnx = rnxdec() +nav = Nav() +proc = process(nav, rnx, config, nep=nep, xyz_ref=xyz_ref) -pos_ref = ecef2pos(xyz_ref) +sig_t = {'G': ['1C', '2W'], 'C': ['1P', '5P']} # GC -# Define signals to be processed -# -gnss = "GC" -sigs = [] -if 'G' in gnss: - sigs.extend([rSigRnx("GC1C"), rSigRnx("GC2W"), - rSigRnx("GL1C"), rSigRnx("GL2W"), - rSigRnx("GS1C"), rSigRnx("GS2W")]) -if 'C' in gnss: - sigs.extend([rSigRnx("CC1P"), rSigRnx("CC5P"), - rSigRnx("CL1P"), rSigRnx("CL5P"), - rSigRnx("CS1P"), rSigRnx("CS5P")]) - -rnx = rnxdec() +sigs, nav.nf = proc.init_sig(sig_t) rnx.setSignals(sigs) -nav = Nav() - # Positioning mode # 0:static, 1:kinematic # @@ -98,186 +78,55 @@ cs = cssr_bds() cs.monlevel = 0 """ -cs = cssr_bds('test_pppbds_ssr.log') +cs = cssr_bds(f'{ttl}_ssr.log') cs.monlevel = 2 """ -# Load ANTEX data for satellites and stations +# Load RINEX OBS file header # +rnx.decode_obsh(obsfile) -if time > epoch2time([2022, 11, 27, 0, 0, 0]): - atxfile = '../data/antex/igs20.atx' -else: - atxfile = '../data/antex/igs14.atx' +# Initialize position +# +ppp = pppos(nav, rnx.pos, f'{ttl}.log') -atx = atxdec() -atx.readpcv(atxfile) +if time < epoch2time([2022, 11, 27, 0, 0, 0]): + config['atxfile'] = '../data/antex/igs14.atx' -# Initialize data structures for results -# -t = np.zeros(nep) -enu = np.ones((nep, 3))*np.nan -sol = np.zeros((nep, 4)) -ztd = np.zeros((nep, 1)) -smode = np.zeros(nep, dtype=int) -nsat = np.zeros((nep, 3), dtype=int) +proc.prepare_signal(obsfile) -# Logging level +# Skip epochs until start time # -nav.monlevel = 1 # TODO: enabled for testing! +obs = rnx.decode_obs() +while time > obs.t and obs.t.time != 0: + obs = rnx.decode_obs() -# Load RINEX OBS file header +# Loop over number of epoch from file start # -if rnx.decode_obsh(obsfile) >= 0: +for ne in range(nep): + week, tow = cs.set_time(obs.t) # set time for reference - # Auto-substitute signals + # Set initial epoch # - rnx.autoSubstituteSignals() + if ne == 0: + proc.init_time(obs.t) - # Initialize position - # - ppp = pppos(nav, rnx.pos, 'test_pppbds.log') - nav.elmin = np.deg2rad(5.0) + msg = decode_msg(v, tow, prn_ref) + if msg is not None: + cs.decode_cssr(msg, 0) - # Get equipment information + # Call PPP module with BDS-PPP corrections # - nav.fout.write("FileName: {}\n".format(obsfile)) - nav.fout.write("Start : {}\n".format(time2str(rnx.ts))) - if rnx.te is not None: - nav.fout.write("End : {}\n".format(time2str(rnx.te))) - nav.fout.write("Receiver: {}\n".format(rnx.rcv)) - nav.fout.write("Antenna : {}\n".format(rnx.ant)) - nav.fout.write("\n") + if (cs.lc[0].cstat & 0xf) == 0xf: # wait for mask/orb/clk/cbias + ppp.process(obs, cs=cs) - # Set satellite PCO/PCV information - # - nav.sat_ant = atx.pcvs + proc.save_output(obs.t, ne, ppp) # save output - # Set receiver PCO/PCV information, check antenna name and exit if unknown - # - # NOTE: comment out the line with 'sys_exit(1)' to continue with zero - # receiver antenna corrections! - # - if 'UNKNOWN' in rnx.ant or rnx.ant.strip() == "": - nav.fout.write("ERROR: missing antenna type in RINEX OBS header!\n") - sys_exit(1) - else: - nav.rcv_ant = searchpcv(atx.pcvr, rnx.ant, rnx.ts) - if nav.rcv_ant is None: - nav.fout.write("ERROR: missing antenna type <{}> in ANTEX file!\n" - .format(rnx.ant)) - sys_exit(1) - - if nav.rcv_ant is None: - nav.fout.write("WARNING: no receiver antenna corrections applied!\n") - nav.fout.write("\n") - - # Print available signals - # - nav.fout.write("Available signals\n") - for sys, sigs in rnx.sig_map.items(): - txt = "{:7s} {}\n".format(sys2str(sys), - ' '.join([sig.str() for sig in sigs.values()])) - nav.fout.write(txt) - nav.fout.write("\n") - - nav.fout.write("Selected signals\n") - for sys, tmp in rnx.sig_tab.items(): - txt = "{:7s} ".format(sys2str(sys)) - for _, sigs in tmp.items(): - txt += "{} ".format(' '.join([sig.str() for sig in sigs])) - nav.fout.write(txt+"\n") - nav.fout.write("\n") - - # Skip epochs until start time + # Get new epoch, exit after last epoch # obs = rnx.decode_obs() - while time > obs.t and obs.t.time != 0: - obs = rnx.decode_obs() - - # Loop over number of epoch from file start - # - for ne in range(nep): - - week, tow = time2gpst(obs.t) - cs.week = week - cs.tow0 = tow//86400*86400 - - # Set initial epoch - # - if ne == 0: - nav.t = deepcopy(obs.t) - t0 = deepcopy(obs.t) - t0.time = t0.time//30*30 - nav.time_p = t0 - - vi = v[(v['tow'] == tow) & (v['prn'] == prn_ref)] - if len(vi) > 0: - buff = unhexlify(vi['nav'][0]) - # prn, rev = bs.unpack_from('u6u6', buff, 0) - cs.decode_cssr(buff, 0) - - # Call PPP module with BDS-PPP corrections - # - if (cs.lc[0].cstat & 0xf) == 0xf: - ppp.process(obs, cs=cs) - - # Save output - # - t[ne] = timediff(nav.t, t0)/86400.0 - - sol = nav.xa[0:3] if nav.smode == 4 else nav.x[0:3] - enu[ne, :] = gn.ecef2enu(pos_ref, sol-xyz_ref) - - ztd[ne] = nav.xa[ppp.IT(nav.na)] \ - if nav.smode == 4 else nav.x[ppp.IT(nav.na)] - smode[ne] = nav.smode - nsat[ne, :] = nav.nsat - - nav.fout.write("{} {:14.4f} {:14.4f} {:14.4f} " - "ENU {:7.3f} {:7.3f} {:7.3f}, 2D {:6.3f}, mode {:1d}\n" - .format(time2str(obs.t), - sol[0], sol[1], sol[2], - enu[ne, 0], enu[ne, 1], enu[ne, 2], - np.sqrt(enu[ne, 0]**2+enu[ne, 1]**2), - smode[ne])) - - # Log to standard output - # - stdout.write('\r {} ENU {:7.3f} {:7.3f} {:7.3f}, 2D {:6.3f}, mode {:1d}' - .format(time2str(obs.t), - enu[ne, 0], enu[ne, 1], enu[ne, 2], - np.sqrt(enu[ne, 0]**2+enu[ne, 1]**2), - smode[ne])) - - # Get new epoch, exit after last epoch - # - obs = rnx.decode_obs() - 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 - -if fig_type == 1: - plot_enu(t, enu, smode, ztd) -elif fig_type == 2: - plot_enu(t, enu, smode, figtype=fig_type) - -plotFileFormat = 'eps' -plotFileName = '.'.join(('test_pppbds', plotFileFormat)) + if obs.t.time == 0: + break -plt.savefig(plotFileName, format=plotFileFormat, bbox_inches='tight', dpi=300) -# plt.show() +proc.close() +proc.plot(ttl, fig_type=1) diff --git a/samples/test_ppphas.py b/samples/test_ppphas.py index d1f92c0..725265d 100644 --- a/samples/test_ppphas.py +++ b/samples/test_ppphas.py @@ -2,28 +2,25 @@ static test for PPP (Galileo HAS SIS) """ -from copy import deepcopy -import matplotlib.pyplot as plt import numpy as np -from sys import exit as sys_exit -from sys import stdout - -import cssrlib.gnss as gn -from cssrlib.gnss import ecef2pos, Nav -from cssrlib.gnss import time2gpst, time2doy, time2str, timediff, epoch2time -from cssrlib.gnss import rSigRnx, prn2sat, uGNSS -from cssrlib.gnss import sys2str -from cssrlib.peph import atxdec, searchpcv +from cssrlib.gnss import Nav, time2doy, epoch2time +from cssrlib.gnss import prn2sat, uGNSS, load_config from cssrlib.cssr_has import cssr_has, cnav_msg from cssrlib.pppssr import pppos from cssrlib.rinex import rnxdec -from cssrlib.plot import plot_enu +from cssrlib.utils import process +config = load_config('config_ppp.yml') # Select test case # +ttl = 'test_ppphas' dataset = 3 +nep = 900*4 + excl_sat = [] +navfile = None +file_has = None # Start epoch and number of epochs # fromSbfConvert = False @@ -43,30 +40,27 @@ elif dataset == 2: ep = [2025, 2, 15, 17, 0, 0] xyz_ref = [-3962108.6836, 3381309.5672, 3668678.6720] - navfile = '../data/doy2025-046/046r_rnx.nav' # Mosaic-X5 - obsfile = '../data/doy2025-046/046r_rnx.obs' # Mosaic-X5 - file_has = '../data/doy2025-046/046r_gale6.txt' excl_sat = [prn2sat(uGNSS.GAL, 29)] # E29 elif dataset == 3: ep = [2025, 8, 21, 7, 0, 0] - navfile = '../data/doy2025-233/233h_rnx.nav' - # navfile = '../data/brdc/BRD400DLR_S_20252330000_01D_MN.rnx' - obsfile = '../data/doy2025-233/233h_rnx.obs' # SEPT MOSAIC-X5 # obsfile = '../data/doy2025-233/sept233h_rnx.obs' # SEPT POLARX5 # obsfile = '../data/doy2025-233/ux2233h_rnx.obs' # u-blox X20P # obsfile = '../data/doy2025-233/jav3233h_rnx.obs' # javad DELTA-3S - file_has = '../data/doy2025-233/233h_gale6.txt' xyz_ref = [-3962108.6836, 3381309.5672, 3668678.6720] # Kamakura # Convert epoch and user reference position # -pos_ref = ecef2pos(xyz_ref) time = epoch2time(ep) year = ep[0] doy = int(time2doy(time)) +ses = chr(ord('a')+ep[3]) -nep = 900*4 +if navfile is None: + navfile = f'../data/doy{year}-{doy:03d}/{doy:03d}{ses}_rnx.nav' + obsfile = f'../data/doy{year}-{doy:03d}/{doy:03d}{ses}_rnx.obs' +if file_has is None: + file_has = f'../data/doy{year}-{doy:03d}/{doy:03d}{ses}_gale6.txt' # Load SSR correction file # @@ -89,26 +83,17 @@ # Define signals to be processed # -gnss = "GE" -sigs = [] -if 'G' in gnss: - if dataset == 3: - sigs.extend([rSigRnx("GC1C"), rSigRnx("GC2L"), - rSigRnx("GL1C"), rSigRnx("GL2L"), - rSigRnx("GS1C"), rSigRnx("GS2L")]) - else: - sigs.extend([rSigRnx("GC1C"), rSigRnx("GC2W"), - rSigRnx("GL1C"), rSigRnx("GL2W"), - rSigRnx("GS1C"), rSigRnx("GS2W")]) -if 'E' in gnss: - sigs.extend([rSigRnx("EC1C"), rSigRnx("EC5Q"), - rSigRnx("EL1C"), rSigRnx("EL5Q"), - rSigRnx("ES1C"), rSigRnx("ES5Q")]) +if dataset == 3: + sig_t = {'G': ['1C', '2L'], 'E': ['1C', '5Q']} # GE +else: + sig_t = {'G': ['1C', '2W'], 'E': ['1C', '5Q']} # GE rnx = rnxdec() -rnx.setSignals(sigs) - nav = Nav() +proc = process(nav, rnx, config, nep=nep, xyz_ref=xyz_ref) + +sigs, nav.nf = proc.init_sig(sig_t) +rnx.setSignals(sigs) # Positioning mode # 0:static, 1:kinematic @@ -136,181 +121,54 @@ cnav = cnav_msg() cnav.load_gmat(file_gm) -# Load ANTEX data for satellites and stations +# Load RINEX OBS file header # -atx = atxdec() -if time > epoch2time([2025, 5, 15, 17, 18, 0]): - atx.readpcv('../data/antex/igs20.atx') -else: - atx.readpcv('../data/antex/has14_2345.atx') -# atx.readpcv('../data/antex/igs20.atx', onlyReceiver=True) +rnx.decode_obsh(obsfile) -# Initialize data structures for results +# Initialize position # -t = np.zeros(nep) -enu = np.ones((nep, 3))*np.nan -sol = np.zeros((nep, 4)) -ztd = np.zeros((nep, 1)) -smode = np.zeros(nep, dtype=int) -nsat = np.zeros((nep, 3), dtype=int) - -# Logging level -# -nav.monlevel = 1 # TODO: enabled for testing! +ppp = pppos(nav, rnx.pos, f'{ttl}.log') -# Load RINEX OBS file header -# -if rnx.decode_obsh(obsfile) >= 0: +if time < epoch2time([2025, 5, 15, 17, 18, 0]): + config['atxfile'] = '../data/antex/has14_2345.atx' - # Auto-substitute signals - # - rnx.autoSubstituteSignals() +proc.prepare_signal(obsfile) - # Initialize position - # - ppp = pppos(nav, rnx.pos, 'test_ppphas.log') - nav.elmin = np.deg2rad(5.0) +# Skip epochs until start time +# +obs = rnx.decode_obs() +while time > obs.t and obs.t.time != 0: + obs = rnx.decode_obs() - # Get equipment information - # - nav.fout.write("FileName: {}\n".format(obsfile)) - nav.fout.write("Start : {}\n".format(time2str(rnx.ts))) - if rnx.te is not None: - nav.fout.write("End : {}\n".format(time2str(rnx.te))) - nav.fout.write("Receiver: {}\n".format(rnx.rcv)) - nav.fout.write("Antenna : {}\n".format(rnx.ant)) - nav.fout.write("\n") - - # Set satellite PCO/PCV information - # - nav.sat_ant = atx.pcvs +# Loop over number of epoch from file start +# +for ne in range(nep): + week, tow = cs.set_time(obs.t) # set time for reference - # Set receiver PCO/PCV information, check antenna name and exit if unknown - # - # NOTE: comment out the line with 'sys_exit(1)' to continue with zero - # receiver antenna corrections! - # - if 'UNKNOWN' in rnx.ant or rnx.ant.strip() == "": - nav.fout.write("ERROR: missing antenna type in RINEX OBS header!\n") - sys_exit(1) - else: - nav.rcv_ant = searchpcv(atx.pcvr, rnx.ant, rnx.ts) - if nav.rcv_ant is None: - nav.fout.write("ERROR: missing antenna type <{}> in ANTEX file!\n" - .format(rnx.ant)) - sys_exit(1) - - if nav.rcv_ant is None: - nav.fout.write("WARNING: no receiver antenna corrections applied!\n") - nav.fout.write("\n") - - # Print available signals + # Set initial epoch # - nav.fout.write("Available signals\n") - for sys, sigs in rnx.sig_map.items(): - txt = "{:7s} {}\n".format(sys2str(sys), - ' '.join([sig.str() for sig in sigs.values()])) - nav.fout.write(txt) - nav.fout.write("\n") - - nav.fout.write("Selected signals\n") - for sys, tmp in rnx.sig_tab.items(): - txt = "{:7s} ".format(sys2str(sys)) - for _, sigs in tmp.items(): - txt += "{} ".format(' '.join([sig.str() for sig in sigs])) - nav.fout.write(txt+"\n") - nav.fout.write("\n") - - # Skip epochs until start time - # - obs = rnx.decode_obs() - while time > obs.t and obs.t.time != 0: - obs = rnx.decode_obs() + if ne == 0: + proc.init_time(obs.t) - # Loop over number of epoch from file start - # - for ne in range(nep): - - week, tow = time2gpst(obs.t) - tow = (tow+0.05)//1 - cs.week = week - cs.tow0 = tow//3600*3600 - - # Set initial epoch - # - if ne == 0: - nav.t = deepcopy(obs.t) - t0 = deepcopy(obs.t) - t0.time = t0.time//30*30 - nav.time_p = t0 - - vi = v[v['tow'] == tow] - - HASmsg = cnav.decode_cnav(tow, vi) # decode CNAV pages - if HASmsg is not None: - cs.msgtype = cnav.msgtype - cs.decode_cssr(HASmsg) # decode HAS messages - - # Call PPP module with HAS corrections - # - if (cs.lc[0].cstat & 0xf) == 0xf: - ppp.process(obs, cs=cs) - - # Save output - # - t[ne] = timediff(nav.t, t0)/86400.0 - - sol = nav.xa[0:3] if nav.smode == 4 else nav.x[0:3] - enu[ne, :] = gn.ecef2enu(pos_ref, sol-xyz_ref) - - ztd[ne] = nav.xa[ppp.IT(nav.na)] \ - if nav.smode == 4 else nav.x[ppp.IT(nav.na)] - smode[ne] = nav.smode - nsat[ne, :] = nav.nsat - - nav.fout.write("{} {:14.4f} {:14.4f} {:14.4f} " - "ENU {:7.3f} {:7.3f} {:7.3f}, 2D {:6.3f}, mode {:1d}\n" - .format(time2str(obs.t), - sol[0], sol[1], sol[2], - enu[ne, 0], enu[ne, 1], enu[ne, 2], - np.sqrt(enu[ne, 0]**2+enu[ne, 1]**2), - smode[ne])) - - # Log to standard output - # - stdout.write('\r {} ENU {:7.3f} {:7.3f} {:7.3f}, 2D {:6.3f}, mode {:1d}' - .format(time2str(obs.t), - enu[ne, 0], enu[ne, 1], enu[ne, 2], - np.sqrt(enu[ne, 0]**2+enu[ne, 1]**2), - smode[ne])) - - # Get new epoch, exit after last epoch - # - obs = rnx.decode_obs() - if obs.t.time == 0: - break - - # Send line-break to stdout - # - stdout.write('\n') + vi = v[v['tow'] == tow] - # Close RINEX observation file - # - rnx.fobs.close() + HASmsg = cnav.decode_cnav(tow, vi) # decode CNAV pages + if HASmsg is not None: + cs.msgtype = cnav.msgtype + cs.decode_cssr(HASmsg) # decode HAS messages - # Close output file + # Call PPP module with HAS corrections # - if nav.fout is not None: - nav.fout.close() + if (cs.lc[0].cstat & 0xf) == 0xf: + ppp.process(obs, cs=cs) -fig_type = 1 + proc.save_output(obs.t, ne, ppp) # save output -if fig_type == 1: - plot_enu(t, enu, smode, ztd) -elif fig_type == 2: - plot_enu(t, enu, smode, figtype=fig_type) + # Get new epoch, exit after last epoch + # + obs = rnx.decode_obs() + if obs.t.time == 0: + break -plotFileFormat = 'eps' # 'eps' or 'png' -plotFileName = '.'.join(('test_ppphas', plotFileFormat)) -plt.savefig(plotFileName, format=plotFileFormat, bbox_inches='tight', dpi=300) -# plt.show() +proc.close() +proc.plot(ttl, fig_type=1) diff --git a/samples/test_pppigs.py b/samples/test_pppigs.py index 587caf9..8516e7a 100644 --- a/samples/test_pppigs.py +++ b/samples/test_pppigs.py @@ -2,25 +2,19 @@ static test for PPP (IGS) """ from copy import deepcopy -import matplotlib.pyplot as plt -import numpy as np -from sys import exit as sys_exit -from sys import stdout - -import cssrlib.gnss as gn -from cssrlib.gnss import ecef2pos, Nav -from cssrlib.gnss import time2doy, time2str, timediff, epoch2time -from cssrlib.gnss import rSigRnx -from cssrlib.gnss import sys2str -from cssrlib.peph import atxdec, searchpcv +from cssrlib.gnss import ecef2pos, Nav, load_config +from cssrlib.gnss import time2doy, epoch2time from cssrlib.peph import peph, biasdec from cssrlib.pppssr import pppos from cssrlib.rinex import rnxdec -from cssrlib.plot import plot_enu +from cssrlib.utils import process + +config = load_config('config_ppp.yml') # Start epoch and number of epochs # dataset = 4 +ttl = 'test_pppigs' if dataset == 0: # SETP078M.21O ep = [2021, 3, 19, 12, 0, 0] @@ -82,35 +76,17 @@ # Define signals to be processed # -gnss = "GEJ" -sigs = [] -if 'G' in gnss: - sigs.extend([rSigRnx("GC1C"), rSigRnx("GC2W"), - rSigRnx("GL1C"), rSigRnx("GL2W"), - rSigRnx("GS1C"), rSigRnx("GS2W")]) -if 'R' in gnss: - sigs.extend([rSigRnx("RC1C"), rSigRnx("RC2P"), - rSigRnx("RL1C"), rSigRnx("RL2P"), - rSigRnx("RS1C"), rSigRnx("RS2P")]) -if 'E' in gnss: - sigs.extend([rSigRnx("EC1C"), rSigRnx("EC5Q"), - rSigRnx("EL1C"), rSigRnx("EL5Q"), - rSigRnx("ES1C"), rSigRnx("ES5Q")]) -if 'C' in gnss: - sigs.extend([rSigRnx("CC2I"), rSigRnx("CC6I"), - rSigRnx("CL2I"), rSigRnx("CL6I"), - rSigRnx("CS2I"), rSigRnx("CS6I")]) -if 'J' in gnss: - sigs.extend([rSigRnx("JC1C"), rSigRnx("JC5Q"), - rSigRnx("JL1C"), rSigRnx("JL5Q"), - rSigRnx("JS1C"), rSigRnx("JS5Q")]) +sig_t = {'G': ['1C', '2W'], 'E': ['1C', '5Q'], 'J': ['1C', '5Q']} # GEJ rnx = rnxdec() -rnx.setSignals(sigs) - nav = Nav() orb = peph() +proc = process(nav, rnx, config, nep=nep, xyz_ref=xyz_ref) + +sigs, nav.nf = proc.init_sig(sig_t) +rnx.setSignals(sigs) + # Set site code for GLONASS code biases # site = "SUWN" @@ -144,164 +120,50 @@ else: atxfile += 'M14.ATX' if 'COD0MGXFIN' in ac else 'igs14.atx' -atx = atxdec() -atx.readpcv(atxfile) - -# Initialize data structures for results -# -t = np.zeros(nep) -enu = np.ones((nep, 3))*np.nan -sol = np.zeros((nep, 4)) -ztd = np.zeros((nep, 1)) -smode = np.zeros(nep, dtype=int) +config['atxfile'] = atxfile -# Logging level -# -nav.monlevel = 1 # TODO: enabled for testing! # Load RINEX OBS file header # -if rnx.decode_obsh(obsfile) >= 0: - - # Auto-substitute signals - # - rnx.autoSubstituteSignals() +rnx.decode_obsh(obsfile) - # Initialize position - # - ppp = pppos(nav, rnx.pos, 'test_pppigs.log') - nav.ephopt = 4 # IGS - nav.armode = 3 # 1: continuous, 3: fix-and-hold - nav.parmode = 1 # 1: normal, 2: partial ambiguity resolution - nav.thresar = 2.0 +# Initialize position +# +ppp = pppos(nav, rnx.pos, f'{ttl}.log') +nav.ephopt = 4 # IGS +nav.armode = 3 # 1: continuous, 3: fix-and-hold +nav.parmode = 1 # 1: normal, 2: partial ambiguity resolution +nav.thresar = 2.0 - nav.elmin = np.deg2rad(10.0) +proc.prepare_signal(obsfile) - # Get equipment information - # - nav.fout.write("FileName: {}\n".format(obsfile)) - nav.fout.write("Start : {}\n".format(time2str(rnx.ts))) - if rnx.te is not None: - nav.fout.write("End : {}\n".format(time2str(rnx.te))) - nav.fout.write("Receiver: {}\n".format(rnx.rcv)) - nav.fout.write("Antenna : {}\n".format(rnx.ant)) - nav.fout.write("\n") - - # Set satellite PCO/PCV information - # - nav.sat_ant = atx.pcvs - - # Set receiver PCO/PCV information, check antenna name and exit if unknown - # - # NOTE: comment out the line with 'sys_exit(1)' to continue with zero - # receiver antenna corrections! - # - if 'UNKNOWN' in rnx.ant or rnx.ant.strip() == "": - nav.fout.write("ERROR: missing antenna type in RINEX OBS header!\n") - sys_exit(1) - else: - nav.rcv_ant = searchpcv(atx.pcvr, rnx.ant, rnx.ts) - if nav.rcv_ant is None: - nav.fout.write("ERROR: missing antenna type <{}> in ANTEX file!\n" - .format(rnx.ant)) - sys_exit(1) - - if nav.rcv_ant is None: - nav.fout.write("WARNING: no receiver antenna corrections applied!\n") - nav.fout.write("\n") - - # Print available signals - # - nav.fout.write("Available signals\n") - for sys, sigs in rnx.sig_map.items(): - txt = "{:7s} {}\n".format(sys2str(sys), - ' '.join([sig.str() for sig in sigs.values()])) - nav.fout.write(txt) - nav.fout.write("\n") - - nav.fout.write("Selected signals\n") - for sys, tmp in rnx.sig_tab.items(): - txt = "{:7s} ".format(sys2str(sys)) - for _, sigs in tmp.items(): - txt += "{} ".format(' '.join([sig.str() for sig in sigs])) - nav.fout.write(txt+"\n") - nav.fout.write("\n") - - # Skip epochs until start time - # +# Skip epochs until start time +# +obs = rnx.decode_obs() +while time > obs.t and obs.t.time != 0: obs = rnx.decode_obs() - while time > obs.t and obs.t.time != 0: - obs = rnx.decode_obs() - # Loop over number of epoch from file start - # - for ne in range(nep): - - # Set initial epoch - # - if ne == 0: - nav.t = deepcopy(obs.t) - t0 = deepcopy(obs.t) - - # Call PPP module with IGS products - # - ppp.process(obs, orb=orb, bsx=bsx) - - # Save output - # - t[ne] = timediff(nav.t, t0)/86400.0 - - sol = nav.xa[0:3] if nav.smode == 4 else nav.x[0:3] - enu[ne, :] = gn.ecef2enu(pos_ref, sol-xyz_ref) - - ztd[ne] = nav.xa[ppp.IT(nav.na)] \ - if nav.smode == 4 else nav.x[ppp.IT(nav.na)] - smode[ne] = nav.smode - - nav.fout.write("{} {:14.4f} {:14.4f} {:14.4f} " - "ENU {:7.3f} {:7.3f} {:7.3f}, 2D {:6.3f}, mode {:1d}\n" - .format(time2str(obs.t), - sol[0], sol[1], sol[2], - enu[ne, 0], enu[ne, 1], enu[ne, 2], - np.sqrt(enu[ne, 0]**2+enu[ne, 1]**2), - smode[ne])) - - # Log to standard output - # - stdout.write('\r {} ENU {:7.3f} {:7.3f} {:7.3f}, 2D {:6.3f}, mode {:1d}' - .format(time2str(obs.t), - enu[ne, 0], enu[ne, 1], enu[ne, 2], - np.sqrt(enu[ne, 0]**2+enu[ne, 1]**2), - smode[ne])) - - # Get new epoch, exit after last epoch - # - obs = rnx.decode_obs() - if obs.t.time == 0: - break - - # Send line-break to stdout - # - stdout.write('\n') +# Loop over number of epoch from file start +# +for ne in range(nep): - # Close RINEX observation file + # Set initial epoch # - rnx.fobs.close() + if ne == 0: + nav.t = deepcopy(obs.t) + t0 = deepcopy(obs.t) - # Close output file + # Call PPP module with IGS products # - if nav.fout is not None: - nav.fout.close() + ppp.process(obs, orb=orb, bsx=bsx) -fig_type = 1 + proc.save_output(obs.t, ne, ppp) # save output -if fig_type == 1: - plot_enu(t, enu, smode, ztd) -elif fig_type == 2: - plot_enu(t, enu, smode, figtype=fig_type) - -plotFileFormat = 'eps' -plotFileName = '.'.join(('test_pppigs', plotFileFormat)) + # Get new epoch, exit after last epoch + # + obs = rnx.decode_obs() + if obs.t.time == 0: + break -plt.savefig(plotFileName, format=plotFileFormat, bbox_inches='tight', dpi=300) -# plt.show() +proc.close() +proc.plot(ttl, fig_type=1) diff --git a/samples/test_pppmdc.py b/samples/test_pppmdc.py index a782e6d..c5eed9d 100644 --- a/samples/test_pppmdc.py +++ b/samples/test_pppmdc.py @@ -2,67 +2,84 @@ static test for PPP (MADOCA PPP) """ from binascii import unhexlify -from copy import deepcopy -import matplotlib.pyplot as plt import numpy as np from sys import exit as sys_exit -from sys import stdout - -from cssrlib.gnss import ecef2pos, ecef2enu, Nav, rSigRnx, sys2str -from cssrlib.gnss import time2gpst, time2doy, time2str, timediff, epoch2time -from cssrlib.peph import atxdec, searchpcv +from cssrlib.gnss import ecef2pos, Nav, load_config +from cssrlib.gnss import time2doy, epoch2time from cssrlib.cssrlib import sCType as sc from cssrlib.cssr_mdc import cssr_mdc from cssrlib.pppssr import pppos from cssrlib.rinex import rnxdec -from cssrlib.plot import plot_enu +from cssrlib.utils import process + +config = load_config('config_ppp.yml') + + +def decode_msg(v, tow, prn_ref, l6_ch=0, prn_ref_ext=0, l6_ch_ext=0): + """ find valid correction message """ + + msg, msg_e = None, None + + vi_ = v[v['tow'] == tow] + vi = vi_[(vi_['type'] == l6_ch) & (vi_['prn'] == prn_ref)] + if len(vi) > 0: + msg = unhexlify(vi['nav'][0]) + + # load regional STEC info (experimental) + if prn_ref_ext > 0: + vi = vi_[(vi_['type'] == l6_ch_ext) & + (vi_['prn'] == prn_ref_ext)] + if len(vi) > 0: + msg_e = unhexlify(vi['nav'][0]) + + return msg, msg_e + # Select test case # +ttl = 'test_pppmdc' l6_mode = 0 # 0: from receiver log, 1: from archive on QZSS -dataset = 3 - +dataset = 1 +navfile = None +file_l6 = None file_stec = None +prn_ref = 199 # QZSS PRN +l6_ch = 1 # 0:L6D, 1:L6E + +prn_ref_ext = -1 +l6_ch_ext = 0 # 0:L6D,1:L6E + # Start epoch and number of epochs # if dataset == 0: - ep = [2023, 7, 8, 4, 0, 0] - xyz_ref = [-3962108.6836, 3381309.5672, 3668678.6720] - navfile = '../data/doy2023-189/SEPT1890.23P' - obsfile = '../data/doy2023-189/SEPT1890.23O' - file_l6 = '../data/doy2023-189/qzsl6_189e.txt' -elif dataset == 1: - ep = [2023, 8, 11, 21, 0, 0] - xyz_ref = [-3962108.6836, 3381309.5672, 3668678.6720] - navfile = '../data/doy2023-223/NAV223.23p' - # navfile = '../data/brdc/BRD400DLR_S_20232230000_01D_MN.rnx' - # obsfile = '../data/doy2023-223/SEPT223Z.23O' # MOSAIC-CLAS - obsfile = '../data/doy2023-223/SEPT223Y.23O' # PolaRX5 - file_l6 = '../data/doy2023-223/223v_qzsl6.txt' -elif dataset == 2: ep = [2025, 2, 15, 17, 0, 0] xyz_ref = [-3962108.6836, 3381309.5672, 3668678.6720] - navfile = '../data/doy2025-046/046r_rnx.nav' # - obsfile = '../data/doy2025-046/046r_rnx.obs' # SEPT MOSAIC-X5 - if l6_mode == 0: - file_l6 = '../data/doy2025-046/046r_qzsl6.txt' - elif l6_mode == 1: + if l6_mode == 1: file_l6 = '../data/doy2025-046/2025046R.l6' -elif dataset == 3: +elif dataset == 1: ep = [2025, 8, 21, 7, 0, 0] - navfile = '../data/doy2025-233/233h_rnx.nav' - # navfile = '../data/brdc/BRD400DLR_S_20252330000_01D_MN.rnx' - obsfile = '../data/doy2025-233/233h_rnx.obs' # SEPT MOSAIC-X5 - file_l6 = '../data/doy2025-233/233h_qzsl6.txt' xyz_ref = [-3962108.6836, 3381309.5672, 3668678.6720] # Kamakura - # file_stec = '../data/qzsl6/2025233H.201.l6' # STEC correction +elif dataset == 2: # MADOCA-PPP with iono correction + ep = [2025, 8, 21, 7, 0, 0] + xyz_ref = [-3962108.6836, 3381309.5672, 3668678.6720] # Kamakura + file_stec = '../data/qzsl6/2025233H.201.l6' # STEC correction + # prn_ref_ext = 201 # QZSS PRN for iono-correction + time = epoch2time(ep) year = ep[0] doy = int(time2doy(time)) +ses = chr(ord('a')+ep[3]) + +if navfile is None: + navfile = f'../data/doy{year}-{doy:03d}/{doy:03d}{ses}_rnx.nav' + obsfile = f'../data/doy{year}-{doy:03d}/{doy:03d}{ses}_rnx.obs' +if file_l6 is None: + file_l6 = f'../data/doy{year}-{doy:03d}/{doy:03d}{ses}_qzsl6.txt' -nep = 900*4-5 +nep = 3600-5 +# nep = 120 if l6_mode == 1: fc = open(file_l6, 'rb') @@ -74,13 +91,6 @@ ('type', 'int'), ('len', 'int'), ('nav', 'S500')] v = np.genfromtxt(file_l6, dtype=dtype) -prn_ref = 199 # QZSS PRN -l6_ch = 1 # 0:L6D, 1:L6E - -prn_ref_ext = -1 -# prn_ref_ext = 201 # QZSS PRN for iono-correction -l6_ch_ext = 0 # 0:L6D,1:L6E - if file_stec is not None: fc = open(file_stec, 'rb') if not fc: @@ -89,40 +99,18 @@ iono_opt = 2 if prn_ref_ext > 0 or (file_stec is not None) else 1 -pos_ref = ecef2pos(xyz_ref) +rnx = rnxdec() +nav = Nav() +proc = process(nav, rnx, config, nep=nep, xyz_ref=xyz_ref) # Define signals to be processed -# -# gnss = "GE" -gnss = "GEJRC" -sigs = [] -if 'G' in gnss: - sigs.extend([rSigRnx("GC1C"), rSigRnx("GC2W"), - rSigRnx("GL1C"), rSigRnx("GL2W"), - rSigRnx("GS1C"), rSigRnx("GS2W")]) -if 'E' in gnss: - sigs.extend([rSigRnx("EC1C"), rSigRnx("EC5Q"), - rSigRnx("EL1C"), rSigRnx("EL5Q"), - rSigRnx("ES1C"), rSigRnx("ES5Q")]) -if 'J' in gnss: - sigs.extend([rSigRnx("JC1C"), rSigRnx("JC2L"), - rSigRnx("JL1C"), rSigRnx("JL2L"), - rSigRnx("JS1C"), rSigRnx("JS2L")]) - -if 'R' in gnss: - sigs.extend([rSigRnx("RC1C"), rSigRnx("RC2C"), - rSigRnx("RL1C"), rSigRnx("RL2C"), - rSigRnx("RS1C"), rSigRnx("RS2C")]) - -if 'C' in gnss: - sigs.extend([rSigRnx("CC2I"), rSigRnx("CC5P"), - rSigRnx("CL2I"), rSigRnx("CL5P"), - rSigRnx("CS2I"), rSigRnx("CS5P")]) -rnx = rnxdec() -rnx.setSignals(sigs) +# sig_t = {'G': ['1C', '2W'], 'E': ['1C', '5Q']} # GE +sig_t = {'G': ['1C', '2W'], 'E': ['1C', '5Q'], 'J': ['1C', '5Q'], + 'R': ['1C', '2C'], 'C': ['2I', '5P']} # GEJRC -nav = Nav() +sigs, nav.nf = proc.init_sig(sig_t) +rnx.setSignals(sigs) # Positioning mode # 0:static, 1:kinematic @@ -149,217 +137,77 @@ cs_.monlevel = 2 """ -# Load ANTEX data for satellites and stations +# Load RINEX OBS file header # -atxfile = '../data/antex/' -if time > epoch2time([2022, 11, 27, 0, 0, 0]): - atxfile += 'igs20.atx' -else: - atxfile += 'igs14.atx' - -atx = atxdec() -atx.readpcv(atxfile) +rnx.decode_obsh(obsfile) -# Initialize data structures for results +# Initialize position # -t = np.zeros(nep) -enu = np.ones((nep, 3))*np.nan -sol = np.zeros((nep, 4)) -ztd = np.zeros((nep, 1)) -smode = np.zeros(nep, dtype=int) -nsat = np.zeros((nep, 3), dtype=int) - -# Logging level +ppp = pppos(nav, rnx.pos, f'{ttl}.log', iono_opt=iono_opt) +# nav.armode = 3 +# nav.thresar = 2.0 + +proc.prepare_signal(obsfile) + +# Skip epochs until start time # -nav.monlevel = 1 # TODO: enabled for testing! +obs = rnx.decode_obs() +while time > obs.t and obs.t.time != 0: + obs = rnx.decode_obs() -# Load RINEX OBS file header +# Loop over number of epoch from file start # -if rnx.decode_obsh(obsfile) >= 0: +for ne in range(nep): + week, tow = cs.set_time(obs.t) # set time for reference - # Auto-substitute signals + # Set initial epoch # - rnx.autoSubstituteSignals() + if ne == 0: + proc.init_time(obs.t) + + if l6_mode == 1: # from log file + cs.decode_l6msg(fc.read(250), 0) + if cs.fcnt == 5: # end of sub-frame + cs.week = week + cs.decode_cssr(cs.buff, 0) + else: # from L6 log + msg, msg_e = decode_msg(v, tow, prn_ref, l6_ch, prn_ref_ext, l6_ch_ext) + if msg is not None: + cs.decode_l6msg(msg, 0) + if cs.fcnt == 5: # end of sub-frame + cs.decode_cssr(bytes(cs.buff), 0) - # Initialize position - # - ppp = pppos(nav, rnx.pos, 'test_pppmdc.log', iono_opt=iono_opt) - # nav.armode = 3 - # nav.thresar = 2.0 + # load regional STEC info (experimental) + if prn_ref_ext > 0 and msg_e is not None: + cs_.decode_l6msg(msg_e, 0) + if cs_.sid == 1: # end of sub-frame + cs.decode_cssr(bytes(cs_.buff_p), 0) - nav.elmin = np.deg2rad(5.0) + if file_stec is not None: # STEC read from file + cs_.decode_l6msg(fc.read(250), 0) + if cs_.sid == 1: # end of sub-frame + cs.decode_cssr(bytes(cs_.buff_p), 0) - nav.glo_ch = rnx.glo_ch + cs.inet = cs.find_grid_index(ecef2pos(nav.x[:3])) - # Get equipment information + # Call PPP module # - nav.fout.write("FileName: {}\n".format(obsfile)) - nav.fout.write("Start : {}\n".format(time2str(rnx.ts))) - if rnx.te is not None: - nav.fout.write("End : {}\n".format(time2str(rnx.te))) - nav.fout.write("Receiver: {}\n".format(rnx.rcv)) - nav.fout.write("Antenna : {}\n".format(rnx.ant)) - nav.fout.write("\n") - - # Set satellite PCO/PCV information - # - nav.sat_ant = atx.pcvs - - # Set receiver PCO/PCV information, check antenna name and exit if unknown - # - # NOTE: comment out the line with 'sys_exit(1)' to continue with zero - # receiver antenna corrections! - # - if 'UNKNOWN' in rnx.ant or rnx.ant.strip() == "": - nav.fout.write("ERROR: missing antenna type in RINEX OBS header!\n") - sys_exit(1) - else: - nav.rcv_ant = searchpcv(atx.pcvr, rnx.ant, rnx.ts) - if nav.rcv_ant is None: - nav.fout.write("ERROR: missing antenna type <{}> in ANTEX file!\n" - .format(rnx.ant)) - sys_exit(1) - - if nav.rcv_ant is None: - nav.fout.write("WARNING: no receiver antenna corrections applied!\n") - nav.fout.write("\n") - - # Print available signals - # - nav.fout.write("Available signals\n") - for sys_, sigs in rnx.sig_map.items(): - txt = "{:7s} {}\n".format(sys2str(sys_), - ' '.join([sig.str() for sig in sigs.values()])) - nav.fout.write(txt) - nav.fout.write("\n") - - nav.fout.write("Selected signals\n") - for sys_, tmp in rnx.sig_tab.items(): - txt = "{:7s} ".format(sys2str(sys_)) - for _, sigs in tmp.items(): - txt += "{} ".format(' '.join([sig.str() for sig in sigs])) - nav.fout.write(txt+"\n") - nav.fout.write("\n") - - # Skip epochs until start time - # - obs = rnx.decode_obs() - while time > obs.t and obs.t.time != 0: - obs = rnx.decode_obs() - - # Loop over number of epoch from file start - # - for ne in range(nep): - - week, tow = time2gpst(obs.t) - cs.week = week - cs.tow0 = tow//3600*3600 - - # Set initial epoch - # - if ne == 0: - nav.t = deepcopy(obs.t) - t0 = deepcopy(obs.t) - t0.time = t0.time//30*30 - nav.time_p = t0 - - if l6_mode == 1: - cs.decode_l6msg(fc.read(250), 0) - if cs.fcnt == 5: # end of sub-frame - cs.week = week - cs.decode_cssr(cs.buff, 0) - else: # multi-channel mode - vi_ = v[v['tow'] == tow] - - vi = vi_[(vi_['type'] == l6_ch) & (vi_['prn'] == prn_ref)] - if len(vi) > 0: - cs.decode_l6msg(unhexlify(vi['nav'][0]), 0) - if cs.fcnt == 5: # end of sub-frame - cs.decode_cssr(bytes(cs.buff), 0) - - # load regional STEC info (experimental) - if prn_ref_ext > 0: - vi = vi_[(vi_['type'] == l6_ch_ext) & - (vi_['prn'] == prn_ref_ext)] - if len(vi) > 0: - cs_.decode_l6msg(unhexlify(vi['nav'][0]), 0) - if cs_.sid == 1: # end of sub-frame - cs.decode_cssr(bytes(cs_.buff_p), 0) - if file_stec is not None: # STEC read from file - cs_.decode_l6msg(fc.read(250), 0) - if cs_.sid == 1: # end of sub-frame - cs.decode_cssr(bytes(cs_.buff_p), 0) - - cs.inet = cs.find_grid_index(ecef2pos(nav.x[:3])) - - # Call PPP module - # - if (cs.lc[0].cstat & 0xf) == 0xf: - if iono_opt == 2: # STEC is available - mask_s = 1 << sc.STEC - if cs.inet > 0 and \ - cs.lc[cs.inet].cstat & mask_s == mask_s: - ppp.process(obs, cs=cs) - else: + if (cs.lc[0].cstat & 0xf) == 0xf: + if iono_opt == 2: # STEC is available + mask_s = 1 << sc.STEC + if cs.inet > 0 and \ + cs.lc[cs.inet].cstat & mask_s == mask_s: ppp.process(obs, cs=cs) + else: + ppp.process(obs, cs=cs) - # Save output - # - t[ne] = timediff(nav.t, t0)/86400.0 - - sol = nav.xa[0:3] if nav.smode == 4 else nav.x[0:3] - enu[ne, :] = ecef2enu(pos_ref, sol-xyz_ref) - - ztd[ne] = nav.xa[ppp.IT(nav.na)] \ - if nav.smode == 4 else nav.x[ppp.IT(nav.na)] - smode[ne] = nav.smode - nsat[ne, :] = nav.nsat - - nav.fout.write("{} {:14.4f} {:14.4f} {:14.4f} " - "ENU {:7.3f} {:7.3f} {:7.3f}, 2D {:6.3f}, mode {:1d}\n" - .format(time2str(obs.t), - sol[0], sol[1], sol[2], - enu[ne, 0], enu[ne, 1], enu[ne, 2], - np.sqrt(enu[ne, 0]**2+enu[ne, 1]**2), - smode[ne])) - - # Log to standard output - # - stdout.write('\r {} ENU {:7.3f} {:7.3f} {:7.3f}, 2D {:6.3f}, mode {:1d}' - .format(time2str(obs.t), - enu[ne, 0], enu[ne, 1], enu[ne, 2], - np.sqrt(enu[ne, 0]**2+enu[ne, 1]**2), - smode[ne])) - - # Get new epoch, exit after last epoch - # - obs = rnx.decode_obs() - if obs.t.time == 0: - break - - # Send line-break to stdout - # - stdout.write('\n') - - # Close RINEX observation file - # - rnx.fobs.close() + proc.save_output(obs.t, ne, ppp) # save output - # Close output file + # Get new epoch, exit after last epoch # - if nav.fout is not None: - nav.fout.close() - -fig_type = 1 - -if fig_type == 1: - plot_enu(t, enu, smode, ztd) -elif fig_type == 2: - plot_enu(t, enu, smode, figtype=fig_type) - -plotFileFormat = 'eps' -plotFileName = '.'.join(('test_pppmdc', plotFileFormat)) + obs = rnx.decode_obs() + if obs.t.time == 0: + break -plt.savefig(plotFileName, format=plotFileFormat, - bbox_inches='tight', dpi=300) -# plt.show() +proc.close() +proc.plot(ttl, fig_type=1) diff --git a/samples/test_ppppvs.py b/samples/test_ppppvs.py index db72b1f..2f400da 100644 --- a/samples/test_ppppvs.py +++ b/samples/test_ppppvs.py @@ -2,39 +2,44 @@ static test for PPP (PVS PPP) """ from binascii import unhexlify -from copy import deepcopy -import matplotlib.pyplot as plt import numpy as np from sys import exit as sys_exit -from sys import stdout - -import cssrlib.gnss as gn -from cssrlib.gnss import ecef2pos, Nav -from cssrlib.gnss import time2gpst, time2doy, time2str, timediff, epoch2time -from cssrlib.gnss import rSigRnx -from cssrlib.gnss import sys2str -from cssrlib.peph import atxdec, searchpcv + +from cssrlib.gnss import Nav, time2doy, timediff, epoch2time, load_config from cssrlib.cssr_pvs import cssr_pvs from cssrlib.pppssr import pppos from cssrlib.rinex import rnxdec from cssrlib.cssr_pvs import decode_sinca_line -from cssrlib.plot import plot_enu +from cssrlib.utils import process + +config = load_config('config_ppp.yml') # Select test case # -dataset = 3 +dataset = 2 +navfile = None +file_pvs = None +ttl = 'test_ppppvs' -# Start epoch and input files -# -if dataset == 0: # SIS - ep = [2023, 11, 4, 2, 0, 0] - navfile = '../data/brdc/BRD400DLR_S_20233080000_01D_MN.rnx' - obsfile = '../data/doy2023-308/308c_rnx.obs' # Mosaic-X5 - file_pvs = '../data/doy2023-308/308c_sbas.txt' - xyz_ref = [-3962108.7007, 3381309.5532, 3668678.6648] +def decode_msg(v, tow, prn_ref, sbas_type=1): + """ find valid correction message """ + msg = None + vi = v[(v['tow'] == tow) & (v['prn'] == prn_ref)] + if sbas_type == 0: # L1 + vi = vi[vi['type'] <= 30] + else: # DFMC L5 + vi = vi[vi['type'] > 30] + + if len(vi) > 0: + msg = unhexlify(vi['nav'][0]) -elif dataset == 1: # DAS + return msg + + +# Start epoch and input files +# +if dataset == 0: # DAS ep = [2025, 4, 20, 5, 0, 0] navfile = '../data/doy2025-110/BRD400DLR_S_20251100000_01D_MN.rnx' @@ -42,59 +47,47 @@ file_pvs = '../data/doy2025-110/DAS2025110f.txt' xyz_ref = [-4052052.9320, 4212835.9496, -2545104.3074] -elif dataset == 2: # SIS +elif dataset == 1: # SIS ep = [2025, 2, 15, 17, 0, 0] - navfile = '../data/doy2025-046/046r_rnx.nav' - obsfile = '../data/doy2025-046/046r_rnx.obs' # SEPT MOSAIC-X5 - file_pvs = '../data/doy2025-046/046r_sbas.txt' xyz_ref = [-3962108.6836, 3381309.5672, 3668678.6720] -elif dataset == 3: # SIS +elif dataset == 2: # SIS ep = [2025, 8, 21, 7, 0, 0] - navfile = '../data/doy2025-233/233h_rnx.nav' - obsfile = '../data/doy2025-233/233h_rnx.obs' # SEPT MOSAIC-X5 - file_pvs = '../data/doy2025-233/233h_sbas.txt' xyz_ref = [-3962108.6836, 3381309.5672, 3668678.6720] # Kamakura -elif dataset == 4: # SIS +elif dataset == 3: # SIS (Australia) ep = [2025, 8, 21, 7, 0, 0] navfile = '../data/doy2025-233/alby233h_rnx.nav' obsfile = '../data/doy2025-233/alby233h_rnx.obs' # SEPT POLARX5 - file_pvs = '../data/doy2025-233/233h_sbas.txt' xyz_ref = [-2441715.2741, 4629128.6896, -3633362.5218] # Albany, AUSTRALIA time = epoch2time(ep) year = ep[0] doy = int(time2doy(time)) +ses = chr(ord('a')+ep[3]) -nep = 900*4 - +if navfile is None: + navfile = f'../data/doy{year}-{doy:03d}/{doy:03d}{ses}_rnx.nav' + obsfile = f'../data/doy{year}-{doy:03d}/{doy:03d}{ses}_rnx.obs' +if file_pvs is None: + file_pvs = f'../data/doy{year}-{doy:03d}/{doy:03d}{ses}_sbas.txt' +nep = 900*4 prn_ref = 122 # satellite PRN for PRN122 sbas_type = 1 # L1: 0, L5: 1 -pos_ref = ecef2pos(xyz_ref) +rnx = rnxdec() +nav = Nav() +proc = process(nav, rnx, config, nep=nep, xyz_ref=xyz_ref) -# Define signals to be processed -# -gnss = "GE" -sigs = [] -if 'G' in gnss: - sigs.extend([rSigRnx("GC1C"), rSigRnx("GC5Q"), - rSigRnx("GL1C"), rSigRnx("GL5Q"), - rSigRnx("GS1C"), rSigRnx("GS5Q")]) -if 'E' in gnss: - sigs.extend([rSigRnx("EC1C"), rSigRnx("EC5Q"), - rSigRnx("EL1C"), rSigRnx("EL5Q"), - rSigRnx("ES1C"), rSigRnx("ES5Q")]) +sig_t = {'G': ['1C', '5Q'], 'E': ['1C', '5Q']} # GE -rnx = rnxdec() +sigs, nav.nf = proc.init_sig(sig_t) rnx.setSignals(sigs) -nav = Nav() # Positioning mode # 0:static, 1:kinematic @@ -107,206 +100,72 @@ # cs = cssr_pvs() # cs.monlevel = 0 -cs = cssr_pvs('test_ppppvs_ssr.log') +cs = cssr_pvs(f'{ttl}_ssr.log') cs.monlevel = 2 -# Load ANTEX data for satellites and stations +# Load RINEX OBS file header # -atx = atxdec() -atx.readpcv('../data/antex/igs20.atx') +rnx.decode_obsh(obsfile) -# Initialize data structures for results +# Initialize position # -t = np.zeros(nep) -enu = np.ones((nep, 3))*np.nan -sol = np.zeros((nep, 4)) -ztd = np.zeros((nep, 1)) -smode = np.zeros(nep, dtype=int) -nsat = np.zeros((nep, 3), dtype=int) - -# Logging level +ppp = pppos(nav, rnx.pos, f'{ttl}.log') + +proc.prepare_signal(obsfile) + +# Skip epochs until start time # -nav.monlevel = 1 # TODO: enabled for testing! +obs = rnx.decode_obs() +while time > obs.t and obs.t.time != 0: + obs = rnx.decode_obs() -# Load RINEX OBS file header +if 'sbas' in file_pvs: # SIS + dtype = [('wn', 'int'), ('tow', 'float'), ('prn', 'int'), + ('type', 'int'), ('marker', 'S2'), ('nav', 'S124')] + v = np.genfromtxt(file_pvs, dtype=dtype) +elif 'DAS' in file_pvs: # DAS + fc = open(file_pvs, 'rt') +else: + print("ERROR: unknown file format for correction data") + sys_exit(1) + +# Loop over number of epoch from file start # -if rnx.decode_obsh(obsfile) >= 0: +for ne in range(nep): - # Auto-substitute signals - # - rnx.autoSubstituteSignals() + week, tow = cs.set_time(obs.t) # set time for reference - # Initialize position + # Set initial epoch # - ppp = pppos(nav, rnx.pos, 'test_ppppvs.log') - nav.elmin = np.deg2rad(10.0) + if ne == 0: + proc.init_time(obs.t) - # nav.q[0:3] = 0.0 # use zero process noise on position + if 'sbas' in file_pvs: # SIS + msg = decode_msg(v, tow, prn_ref, sbas_type) + if msg is not None: + cs.decode_cssr(msg, 0) - # Get equipment information - # - nav.fout.write("FileName: {}\n".format(obsfile)) - nav.fout.write("Start : {}\n".format(time2str(rnx.ts))) - if rnx.te is not None: - nav.fout.write("End : {}\n".format(time2str(rnx.te))) - nav.fout.write("Receiver: {}\n".format(rnx.rcv)) - nav.fout.write("Antenna : {}\n".format(rnx.ant)) - nav.fout.write("\n") - - # Set satellite PCO/PCV information - # - nav.sat_ant = atx.pcvs + else: # DAS + for line in fc: + tc, buff = decode_sinca_line(line) + cs.decode_cssr(buff, 0) + if timediff(obs.t, tc) >= 0.0: + break - # Set receiver PCO/PCV information, check antenna name and exit if unknown - # - # NOTE: comment out the line with 'sys_exit(1)' to continue with zero - # receiver antenna corrections! - # - if 'UNKNOWN' in rnx.ant or rnx.ant.strip() == "": - nav.fout.write("ERROR: missing antenna type in RINEX OBS header!\n") - sys_exit(1) - else: - nav.rcv_ant = searchpcv(atx.pcvr, rnx.ant, rnx.ts) - if nav.rcv_ant is None: - nav.fout.write("ERROR: missing antenna type <{}> in ANTEX file!\n" - .format(rnx.ant)) - sys_exit(1) - - if nav.rcv_ant is None: - nav.fout.write("WARNING: no receiver antenna corrections applied!\n") - nav.fout.write("\n") - - # Print available signals - # - nav.fout.write("Available signals\n") - for sys, sigs in rnx.sig_map.items(): - txt = "{:7s} {}\n".format(sys2str(sys), - ' '.join([sig.str() for sig in sigs.values()])) - nav.fout.write(txt) - nav.fout.write("\n") - - nav.fout.write("Selected signals\n") - for sys, tmp in rnx.sig_tab.items(): - txt = "{:7s} ".format(sys2str(sys)) - for _, sigs in tmp.items(): - txt += "{} ".format(' '.join([sig.str() for sig in sigs])) - nav.fout.write(txt+"\n") - nav.fout.write("\n") - - # Skip epochs until start time - # - obs = rnx.decode_obs() - while time > obs.t and obs.t.time != 0: - obs = rnx.decode_obs() + cs.check_validity(obs.t) - if 'sbas' in file_pvs: # SIS - dtype = [('wn', 'int'), ('tow', 'float'), ('prn', 'int'), - ('type', 'int'), ('marker', 'S2'), ('nav', 'S124')] - v = np.genfromtxt(file_pvs, dtype=dtype) - elif 'DAS' in file_pvs: # DAS - fc = open(file_pvs, 'rt') - else: - print("ERROR: unknown file format for correction data") - sys_exit(1) - - # Loop over number of epoch from file start - # - for ne in range(nep): - - week, tow = time2gpst(obs.t) - cs.week = week - cs.tow0 = tow//86400*86400 - cs.time0 = obs.t - - # Set initial epoch - # - if ne == 0: - nav.t = deepcopy(obs.t) - t0 = deepcopy(obs.t) - t0.time = t0.time//30*30 - nav.time_p = t0 - - if 'sbas' in file_pvs: # SIS - - vi = v[(v['tow'] == tow) & (v['prn'] == prn_ref)] - if sbas_type == 0: # L1 - vi = vi[vi['type'] <= 30] - else: # DFMC L5 - vi = vi[vi['type'] > 30] - - if len(vi) > 0: - buff = unhexlify(vi['nav'][0]) - cs.decode_cssr(buff, 0) - else: # DAS - for line in fc: - tc, buff = decode_sinca_line(line) - cs.decode_cssr(buff, 0) - if timediff(obs.t, tc) >= 0.0: - break - - cs.check_validity(obs.t) - - # Call PPP module with PVS corrections - # - if (cs.lc[0].cstat & 0x6) == 0x6: - ppp.process(obs, cs=cs) - - # Save output - # - t[ne] = timediff(nav.t, t0)/86400.0 - - sol = nav.xa[0:3] if nav.smode == 4 else nav.x[0:3] - enu[ne, :] = gn.ecef2enu(pos_ref, sol-xyz_ref) - - ztd[ne] = nav.xa[ppp.IT(nav.na)] \ - if nav.smode == 4 else nav.x[ppp.IT(nav.na)] - smode[ne] = nav.smode - nsat[ne, :] = nav.nsat - - nav.fout.write("{} {:14.4f} {:14.4f} {:14.4f} " - "ENU {:7.3f} {:7.3f} {:7.3f}, 2D {:6.3f}, mode {:1d}\n" - .format(time2str(obs.t), - sol[0], sol[1], sol[2], - enu[ne, 0], enu[ne, 1], enu[ne, 2], - np.sqrt(enu[ne, 0]**2+enu[ne, 1]**2), - smode[ne])) - - # Log to standard output - # - stdout.write('\r {} ENU {:7.3f} {:7.3f} {:7.3f}, 2D {:6.3f}, mode {:1d}' - .format(time2str(obs.t), - enu[ne, 0], enu[ne, 1], enu[ne, 2], - np.sqrt(enu[ne, 0]**2+enu[ne, 1]**2), - smode[ne])) - - # Get new epoch, exit after last epoch - # - obs = rnx.decode_obs() - if obs.t.time == 0: - break - - # Send line-break to stdout + # Call PPP module with PVS corrections # - stdout.write('\n') + if (cs.lc[0].cstat & 0x6) == 0x6: + ppp.process(obs, cs=cs) - # Close RINEX observation file - # - rnx.fobs.close() + proc.save_output(obs.t, ne, ppp) # save output - # Close output file + # Get new epoch, exit after last epoch # - if nav.fout is not None: - nav.fout.close() - -fig_type = 1 - -if fig_type == 1: - plot_enu(t, enu, smode, ztd) -elif fig_type == 2: - plot_enu(t, enu, smode, figtype=fig_type) - -plotFileFormat = 'png' -plotFileName = '.'.join(('test_ppppvs', plotFileFormat)) + obs = rnx.decode_obs() + if obs.t.time == 0: + break -plt.savefig(plotFileName, format=plotFileFormat, bbox_inches='tight', dpi=300) -# plt.show() +proc.close() +proc.plot(ttl, fig_type=1) diff --git a/samples/test_ppprtcm.py b/samples/test_ppprtcm.py index 607e644..1c8a74b 100644 --- a/samples/test_ppprtcm.py +++ b/samples/test_ppprtcm.py @@ -1,33 +1,30 @@ """ - static test for PPP (Galileo HAS IDD) + static test for PPP (Galileo HAS IDD, JPL GDGPS HAS) """ import os -from copy import deepcopy -import matplotlib.pyplot as plt -import matplotlib.dates as md -import numpy as np -from sys import exit as sys_exit -from sys import stdout - -import cssrlib.gnss as gn -from cssrlib.gnss import ecef2pos, Nav -from cssrlib.gnss import time2gpst, time2doy, time2str, timediff, epoch2time -from cssrlib.gnss import rSigRnx, sys2str -from cssrlib.cssrlib import sCSSRTYPE -from cssrlib.peph import atxdec, searchpcv +from cssrlib.gnss import Nav, load_config, time2doy, timediff, epoch2time, timeadd +from cssrlib.cssrlib import sCSSRTYPE, sCType from cssrlib.rtcm import rtcm from cssrlib.pppssr import pppos from cssrlib.rinex import rnxdec -from cssrlib.cssrlib import sCType +from cssrlib.utils import process +config = load_config('config_ppp.yml') # Select test case # icase = 3 +nep = 900*4-5 +# nep = 60*16 +ttl = 'test_ppprtcm' +navfile = None +file_rtcm = None # Start epoch and number of epochs # +cs_mask = 1 << sCType.CLOCK | 1 << sCType.ORBIT | 1 << sCType.CBIAS + if icase == 1: # Galileo HAS IDD ep = [2023, 8, 17, 2, 0, 0] @@ -37,75 +34,60 @@ xyz_ref = [4186704.2262, 834903.7677, 4723664.9337] file_rtcm = '../data/doy2023-229/idd2023229c.rtc' file_rtcm_log = '../data/doy2023-229/idd2023229c.log' - gnss = "GE" - cs_mask = 1 << sCType.CLOCK | 1 << sCType.ORBIT | 1 << sCType.CBIAS elif icase == 2: # JPL GDGPS Mosaic-X5 ep = [2024, 2, 12, 7, 0, 0] - navfile = '../data/doy2024-043/043h_rnx.nav' - # navfile = '../data/brdc/BRD400DLR_S_20240430000_01D_MN.rnx' - obsfile = '../data/doy2024-043/043h_rnx.obs' xyz_ref = [-3962108.7007, 3381309.5532, 3668678.6648] file_rtcm = '../data/doy2024-043/JPL32T2043h.rtcm3' file_rtcm_log = '../data/doy2024-043/JPL32T2043h.log' - gnss = "GE" - cs_mask = 1 << sCType.CLOCK | 1 << sCType.ORBIT | 1 << sCType.CBIAS elif icase == 3: # JPL GDGPS (w/o code bias) JAVAD DELTA-3S ep = [2025, 8, 21, 7, 0, 0] - navfile = '../data/doy2025-233/233h_rnx.nav' - obsfile = '../data/doy2025-233/233h_rnx.obs' xyz_ref = [-3962108.6836, 3381309.5672, 3668678.6720] # SSRA11JPL0 GPS+GAL orbit+clock corrs file_rtcm = '../data/doy2025-233/jpl233h.rtcm3' file_rtcm_log = '../data/doy2025-233/jpl233h.log' - gnss = "GE" cs_mask = 1 << sCType.CLOCK | 1 << sCType.ORBIT +elif icase == 4: # RTCM-SC134 + + ep = [2026, 2, 6,20, 0, 39] + # xyz_ref = [4649341.5804, 1029778.2736, 4229073.4732] + xyz_ref = [4649331.6995, 1029774.1231, 4229062.7567] + # SSRA11JPL0 GPS+GAL orbit+clock corrs + bdir = '../data/sc134/msg/ssr/' + file_rtcm = bdir+'SSRTEST_20260206_CORR_v123.bin' + file_rtcm_log = bdir+'SSRTEST_20260206_CORR_v123.log' + #navfile = bdir+'ROMA037.nav' + navfile = bdir+'M0SE00ITA_R_20260370000_01D_MN.rnx' + obsfile = bdir+'ROVR037.obs' + time = epoch2time(ep) year = ep[0] doy = int(time2doy(time)) +ses = chr(ord('a')+ep[3]) -nep = 900*4 - - -# Set user reference position -# -pos_ref = ecef2pos(xyz_ref) +if navfile is None: + navfile = f'../data/doy{year}-{doy:03d}/{doy:03d}{ses}_rnx.nav' + obsfile = f'../data/doy{year}-{doy:03d}/{doy:03d}{ses}_rnx.obs' # Define signals to be processed # - -sigs = [] - if icase in [1, 2]: + sig_t = {'G': ['1C', '2W'], 'E': ['1C', '7Q']} # GE - if 'G' in gnss: - sigs.extend([rSigRnx("GC1C"), rSigRnx("GC2W"), - rSigRnx("GL1C"), rSigRnx("GL2W"), - rSigRnx("GS1C"), rSigRnx("GS2W")]) - if 'E' in gnss: - sigs.extend([rSigRnx("EC1C"), rSigRnx("EC7Q"), - rSigRnx("EL1C"), rSigRnx("EL7Q"), - rSigRnx("ES1C"), rSigRnx("ES7Q")]) elif icase in [3, 4]: - - if 'G' in gnss: - sigs.extend([rSigRnx("GC1W"), rSigRnx("GC2W"), - rSigRnx("GL1W"), rSigRnx("GL2W"), - rSigRnx("GS1W"), rSigRnx("GS2W")]) - if 'E' in gnss: - sigs.extend([rSigRnx("EC1C"), rSigRnx("EC7Q"), - rSigRnx("EL1C"), rSigRnx("EL7Q"), - rSigRnx("ES1C"), rSigRnx("ES7Q")]) + sig_t = {'G': ['1W', '2W'], 'E': ['1C', '7Q']} # GE rnx = rnxdec() -rnx.setSignals(sigs) - nav = Nav() +proc = process(nav, rnx, config, nep=nep, xyz_ref=xyz_ref) + +sigs, nav.nf = proc.init_sig(sig_t) +rnx.setSignals(sigs) # Positioning mode # 0:static, 1:kinematic @@ -121,6 +103,8 @@ cs.cssrmode = sCSSRTYPE.RTCM3_SSR cs.inet = 0 +# cs.workaround_sc134 = True + if icase in [2, 3]: # mask phase-bias for JPL GDGPS cs.mask_pbias = True @@ -134,252 +118,70 @@ maxlen = len(msg)-5 fc.close() -# Load ANTEX data for satellites and stations -# -atxfile = '../data/antex/igs20.atx' -atx = atxdec() -atx.readpcv(atxfile) - -# Initialize data structures for results -# -t = np.zeros(nep) -enu = np.ones((nep, 3))*np.nan -sol = np.zeros((nep, 4)) -ztd = np.zeros((nep, 1)) -smode = np.zeros(nep, dtype=int) -nsat = np.zeros((nep, 3), dtype=int) - -# Logging level -# -nav.monlevel = 1 # TODO: enabled for testing! - # Load RINEX OBS file header # -if rnx.decode_obsh(obsfile) >= 0: +rnx.decode_obsh(obsfile) - # Auto-substitute signals - # - rnx.autoSubstituteSignals() - # Initialize position - # - ppp = pppos(nav, rnx.pos, 'test_ppprtcm.log') - nav.elmin = np.deg2rad(5.0) +# Initialize position +# +ppp = pppos(nav, rnx.pos, f'{ttl}.log') - # Get equipment information - # - nav.fout.write("FileName: {}\n".format(obsfile)) - nav.fout.write("Start : {}\n".format(time2str(rnx.ts))) - if rnx.te is not None: - nav.fout.write("End : {}\n".format(time2str(rnx.te))) - nav.fout.write("Receiver: {}\n".format(rnx.rcv)) - nav.fout.write("Antenna : {}\n".format(rnx.ant)) - nav.fout.write("\n") - - # Set satellite PCO/PCV information - # - nav.sat_ant = atx.pcvs +proc.prepare_signal(obsfile) - # Set receiver PCO/PCV information, check antenna name and exit if unknown - # - # NOTE: comment out the line with 'sys_exit(1)' to continue with zero - # receiver antenna corrections! - # - if 'UNKNOWN' in rnx.ant or rnx.ant.strip() == "": - nav.fout.write("ERROR: missing antenna type in RINEX OBS header!\n") - sys_exit(1) - else: - nav.rcv_ant = searchpcv(atx.pcvr, rnx.ant, rnx.ts) - if nav.rcv_ant is None: - nav.fout.write("ERROR: missing antenna type <{}> in ANTEX file!\n" - .format(rnx.ant)) - sys_exit(1) - - if nav.rcv_ant is None: - nav.fout.write("WARNING: no receiver antenna corrections applied!\n") - nav.fout.write("\n") - - # Print available signals - # - nav.fout.write("Available signals\n") - for sys, sigs in rnx.sig_map.items(): - txt = "{:7s} {}\n".format(sys2str(sys), - ' '.join([sig.str() for sig in sigs.values()])) - nav.fout.write(txt) - nav.fout.write("\n") - - nav.fout.write("Selected signals\n") - for sys, tmp in rnx.sig_tab.items(): - txt = "{:7s} ".format(sys2str(sys)) - for _, sigs in tmp.items(): - txt += "{} ".format(' '.join([sig.str() for sig in sigs])) - nav.fout.write(txt+"\n") - nav.fout.write("\n") - - # Skip epochs until start time - # +# Skip epochs until start time +# +obs = rnx.decode_obs() +while time > obs.t and obs.t.time != 0: obs = rnx.decode_obs() - while time > obs.t and obs.t.time != 0: - obs = rnx.decode_obs() - k = 0 - # Loop over number of epoch from file start +cs.time = obs.t +k = 0 +# Loop over number of epoch from file start +# +for ne in range(nep): + week, tow = cs.set_time(obs.t) # set time for reference + + # Set initial epoch # - for ne in range(nep): - - week, tow = time2gpst(obs.t) - cs.week = week - cs.tow0 = tow//3600*3600 - - # Set initial epoch - # - if ne == 0: - nav.t = deepcopy(obs.t) - t0 = deepcopy(obs.t) - t0.time = t0.time//30*30 - nav.time_p = t0 - - while True: - stat = cs.sync(msg, k) - if stat is False: - k += 1 - continue - if not cs.checksum(msg, k, maxlen): - k += 1 - continue - - tc = cs.decode_time(msg[k:k+cs.len+3]) - if (tc is not False) and timediff(tc, obs.t) > 0: - break - - _, _, eph, geph, seph = cs.decode(msg[k:k+cs.len+3]) + if ne == 0: + # nav.t = deepcopy(obs.t) + # t0 = deepcopy(obs.t) + # t0.time = t0.time//30*30 + # nav.time_p = t0 + proc.init_time(obs.t) + + while k 0: break - # Send line-break to stdout - # - stdout.write('\n') + _, _, eph, geph, seph = cs.decode(msg[k:k+cs.len+3]) + k += cs.dlen - # Close RINEX observation file + if cs.msgtype in cs.eph_t.values(): + nav.eph.append(eph) + + # Call PPP module with HAS corrections # - rnx.fobs.close() + if (cs.lc[0].cstat & cs_mask) == cs_mask: + ppp.process(obs, cs=cs) + + proc.save_output(obs.t, ne, ppp) # save output - # Close output file + # Get new epoch, exit after last epoch # - if nav.fout is not None: - nav.fout.close() - -fig_type = 1 -ylim = 1.0 - -idx4 = np.where(smode == 4)[0] -idx5 = np.where(smode == 5)[0] -idx0 = np.where(smode == 0)[0] - -fig = plt.figure(figsize=[7, 9]) -fig.set_rasterized(True) - -fmt = '%H:%M' -col_t = ['r', 'y', 'g'] - -if fig_type == 1: - - lbl_t = ['East [m]', 'North [m]', 'Up [m]'] - # nm = 4 - nm = 3 - - for k in range(3): - plt.subplot(nm, 1, k+1) - plt.plot(t[idx0], enu[idx0, k], color=col_t[0], - marker='.', label=None if nm > 3 else 'none') - plt.plot(t[idx5], enu[idx5, k], color=col_t[1], - marker='.', label=None if nm > 3 else 'float') - plt.plot(t[idx4], enu[idx4, k], color=col_t[2], - marker='.', label=None if nm > 3 else 'fix') - - plt.ylabel(lbl_t[k]) - plt.grid() - plt.ylim([-ylim, ylim]) - plt.gca().xaxis.set_major_formatter(md.DateFormatter(fmt)) - if nm < 4: - plt.legend() - - if nm > 3: - plt.subplot(nm, 1, 4) - plt.plot(t[idx0], ztd[idx0]*1e2, color=col_t[0], - marker='.', markersize=8, label='none') - plt.plot(t[idx5], ztd[idx5]*1e2, color=col_t[1], - marker='.', markersize=8, label='float') - plt.plot(t[idx4], ztd[idx4]*1e2, color=col_t[2], - marker='.', markersize=8, label='fix') - plt.ylabel('ZTD [cm]') - plt.grid() - plt.gca().xaxis.set_major_formatter(md.DateFormatter(fmt)) - plt.legend() - - plt.xlabel('Time [HH:MM]') - -elif fig_type == 2: - - ax = fig.add_subplot(111) - - plt.plot(enu[idx0, 0], enu[idx0, 1], - color=col_t[0], marker='.', label='none') - plt.plot(enu[idx5, 0], enu[idx5, 1], - color=col_t[1], marker='.', label='float') - plt.plot(enu[idx4, 0], enu[idx4, 1], - color=col_t[2], marker='.', label='fix') - - plt.xlabel('Easting [m]') - plt.ylabel('Northing [m]') - plt.grid() - plt.axis('equal') - plt.legend() - # ax.set(xlim=(-ylim, ylim), ylim=(-ylim, ylim)) - -plotFileFormat = 'png' -plotFileName = '.'.join(('test_ppprtcm', plotFileFormat)) - -plt.savefig(plotFileName, format=plotFileFormat, bbox_inches='tight', dpi=300) -# plt.show() + obs = rnx.decode_obs() + if obs.t.time == 0: + break + +proc.close() +proc.plot(ttl, fig_type=1) diff --git a/samples/test_ppprtk.py b/samples/test_ppprtk.py index 2e08c50..83af557 100644 --- a/samples/test_ppprtk.py +++ b/samples/test_ppprtk.py @@ -1,24 +1,44 @@ """ static test for PPP-RTK (QZSS CLAS) """ -from copy import deepcopy -import matplotlib.pyplot as plt import numpy as np from sys import exit as sys_exit -from sys import stdout -import cssrlib.gnss as gn from cssrlib.cssrlib import cssr -from cssrlib.gnss import ecef2pos, Nav, time2gpst, timediff, time2str, time2doy -from cssrlib.gnss import rSigRnx, sys2str, epoch2time -from cssrlib.peph import atxdec, searchpcv +from cssrlib.gnss import ecef2pos, Nav, time2gpst, time2doy +from cssrlib.gnss import epoch2time, load_config from cssrlib.ppprtk import ppprtkpos from cssrlib.rinex import rnxdec from binascii import unhexlify -from cssrlib.plot import plot_enu +from cssrlib.utils import process + + +def decode_msg(v, tow, l6_ch, prn_ref): + """ find valid correction message """ + + msg = None + vi = v[(v['tow'] == tow) & (v['type'] == l6_ch) & (v['prn'] == prn_ref)] + if len(vi) > 0: + msg = unhexlify(vi['nav'][0]) + + return msg + + +config = load_config('config_ppprtk.yml') l6_mode = 0 # 0: from receiver log, 1: from archive on QZSS -dataset = 2 +dataset = 3 +nep = 900*4 +# nep = 60 + +navfile = None +file_l6 = None +ttl = 'test_ppprtk' +l6_ch = 0 # 0:L6D, 1:L6E +prn_p1 = 199 +prn_p2 = -1 + +sig_t = {'G': ['1C', '2W'], 'E': ['1C', '5Q'], 'J': ['1C', '2L']} # GEJ if l6_mode == 1: # from archive @@ -35,17 +55,6 @@ ep = [2025, 2, 15, 17, 0, 0] xyz_ref = [-3962108.6836, 3381309.5672, 3668678.6720] - time = epoch2time(ep) - year = ep[0] - doy = int(time2doy(time)) - let = chr(ord('a')+ep[3]) - - bdir = '../data/doy{:04d}-{:03d}/'.format(year, doy) - - navfile = bdir+'{:03d}{}_rnx.nav'.format(doy, let) - obsfile = bdir+'{:03d}{}_rnx.obs'.format(doy, let) # SEPT MOSAIC-X5 - l6file = bdir+'{:04d}{:03d}{}.l6'.format(year, doy, let.upper()) - else: # from receiver log if dataset == 0: @@ -54,256 +63,140 @@ xyz_ref = [-3962108.7007, 3381309.5532, 3668678.6648] navfile = '../data/doy2023-223/NAV223.23p' obsfile = '../data/doy2023-223/SEPT223Y.23O' # PolaRX5 - file_l6 = '../data/doy2023-223/223v_qzsl6.txt' - elif dataset == 2: + elif dataset == 1: # single channel + + ep = [2025, 8, 21, 7, 0, 0] + xyz_ref = [-3962108.6836, 3381309.5672, 3668678.6720] + + elif dataset == 2: # two channel ep = [2025, 8, 21, 7, 0, 0] xyz_ref = [-3962108.6836, 3381309.5672, 3668678.6720] - time = epoch2time(ep) - year = ep[0] - doy = int(time2doy(time)) - let = chr(ord('a')+ep[3]) + # from Tab 4.1.1-1 of IS-QZSS-L6 + prn_p1 = 199 # QZSS PRN pattern 1 (195, 197, 199) + prn_p2 = 194 # QZSS PRN pattern 2 (194, 196) - bdir = '../data/doy{:04d}-{:03d}/'.format(year, doy) + elif dataset == 3: # single channel - navfile = bdir+'{:03d}{}_rnx.nav'.format(doy, let) - obsfile = bdir+'{:03d}{}_rnx.obs'.format(doy, let) # SEPT MOSAIC-X5 - file_l6 = bdir+'{:03d}{}_qzsl6.txt'.format(doy, let) + ep = [2025, 8, 21, 7, 0, 0] + xyz_ref = [-3962108.6836, 3381309.5672, 3668678.6720] - # obsfile = bdir+'ux2233h_rnx.obs' # u-blox X20P - # navfile = '../data/brdc/BRD400DLR_S_20252330000_01D_MN.rnx' +time = epoch2time(ep) +year = ep[0] +doy = int(time2doy(time)) +ses = chr(ord('a')+ep[3]) + +if navfile is None: + navfile = f'../data/doy{year}-{doy:03d}/{doy:03d}{ses}_rnx.nav' + obsfile = f'../data/doy{year}-{doy:03d}/{doy:03d}{ses}_rnx.obs' +if file_l6 is None: + file_l6 = f'../data/doy{year}-{doy:03d}/{doy:03d}{ses}_qzsl6.txt' + if l6_mode == 1: + file_l6 = f'../data/doy{year}-{doy:03d}/{year}{doy:03d}{ses.upper()}.txt' - # from Tab 4.1.1-1 of IS-QZSS-L6 - prn_p1 = 199 # QZSS PRN pattern 1 (195, 197, 199) - prn_p2 = 194 # QZSS PRN pattern 2 (194, 196) - l6_ch = 0 # 0:L6D, 1:L6E time = epoch2time(ep) -atxfile = '../data/antex/' -if time > epoch2time([2022, 11, 27, 0, 0, 0]): - atxfile += 'igs20.atx' -else: - atxfile += 'igs14.atx' -griddef = '../data/clas_grid.def' +if time < epoch2time([2022, 11, 27, 0, 0, 0]): + config['atxfile'] = '../data/antex/igs14.atx' + +griddef = config['griddef'] pos_ref = ecef2pos(xyz_ref) cs = cssr() -cs.monlevel = 1 +cs.monlevel = config['cs']['monlevel'] cs.week = time2gpst(time)[0] cs.read_griddef(griddef) cs_ = cssr() -cs_.monlevel = 1 +cs_.monlevel = config['cs']['monlevel'] cs_.week = cs.week cs_.read_griddef(griddef) -# Define signals to be processed -# -gnss = "GEJ" # "GEJ" -sigs = [] -if 'G' in gnss: - sigs.extend([rSigRnx("GC1C"), rSigRnx("GC2W"), - rSigRnx("GL1C"), rSigRnx("GL2W"), - rSigRnx("GS1C"), rSigRnx("GS2W")]) -if 'E' in gnss: - sigs.extend([rSigRnx("EC1C"), rSigRnx("EC5Q"), - rSigRnx("EL1C"), rSigRnx("EL5Q"), - rSigRnx("ES1C"), rSigRnx("ES5Q")]) -if 'J' in gnss: - sigs.extend([rSigRnx("JC1C"), rSigRnx("JC2L"), - rSigRnx("JL1C"), rSigRnx("JL2L"), - rSigRnx("JS1C"), rSigRnx("JS2L")]) - rnx = rnxdec() +nav = Nav() +proc = process(nav, rnx, config, nep=nep, xyz_ref=xyz_ref) + +sigs, nav.nf = proc.init_sig(sig_t) rnx.setSignals(sigs) -nav = Nav() nav = rnx.decode_nav(navfile, nav) -nep = 900*4 +rnx.decode_obsh(obsfile) -# Load ANTEX data for satellites and stations +# Initialize position # -atx = atxdec() -atx.readpcv(atxfile) - -t = np.zeros(nep) -enu = np.ones((nep, 3))*np.nan -sol = np.zeros((nep, 4)) -smode = np.zeros(nep, dtype=int) -nsat = np.zeros((nep, 3), dtype=int) +ppprtk = ppprtkpos(nav, rnx.pos, logfile=f'{ttl}.log', config=config) -if rnx.decode_obsh(obsfile) >= 0: +proc.prepare_signal(obsfile) - # Auto-substitute signals - # - rnx.autoSubstituteSignals() +# Get grid location +# +pos = ecef2pos(rnx.pos) +inet = cs.find_grid_index(pos) + +if l6_mode == 1: + fc = open(file_l6, 'rb') + if not fc: + nav.fout.write(f"ERROR: cannot open L6 message file {file_l6}!") + sys_exit(-1) +else: + dtype = [('wn', 'int'), ('tow', 'int'), ('prn', 'int'), + ('type', 'int'), ('len', 'int'), ('nav', 'S500')] + v = np.genfromtxt(file_l6, dtype=dtype) - # Initialize position - # - ppprtk = ppprtkpos(nav, rnx.pos, 'test_ppprtk.log') +# Skip epoch until start time +# +obs = rnx.decode_obs() +while time > obs.t and obs.t.time != 0: + obs = rnx.decode_obs() - # Get equipment information - # - nav.fout.write("FileName: {}\n".format(obsfile)) - nav.fout.write("Start : {}\n".format(time2str(rnx.ts))) - if rnx.te is not None: - nav.fout.write("End : {}\n".format(time2str(rnx.te))) - nav.fout.write("Receiver: {}\n".format(rnx.rcv)) - nav.fout.write("Antenna : {}\n".format(rnx.ant)) - nav.fout.write("\n") - - # Set satellite PCO/PCV information - # - nav.sat_ant = atx.pcvs +msg, msg2 = None, None - # Set receiver PCO/PCV information, check antenna name and exit if unknown - # - # NOTE: comment out the line with 'sys_exit(1)' to continue with zero - # receiver antenna corrections! - # - if 'UNKNOWN' in rnx.ant or rnx.ant.strip() == "": - nav.fout.write("ERROR: missing antenna type in RINEX OBS header!\n") - sys_exit(1) - else: - nav.rcv_ant = searchpcv(atx.pcvr, rnx.ant, rnx.ts) - if nav.rcv_ant is None: - nav.fout.write("ERROR: missing antenna type <{}> in ANTEX file!\n" - .format(rnx.ant)) - sys_exit(1) +for ne in range(nep): + week, tow = cs.set_time(obs.t) # set time for reference - if nav.rcv_ant is None: - nav.fout.write("WARNING: no receiver antenna corrections applied!\n") - nav.fout.write("\n") + if ne == 0: + proc.init_time(obs.t) - # Print available signals - # - nav.fout.write("Available signals\n") - for sys, sigs in rnx.sig_map.items(): - txt = "{:7s} {}\n".format(sys2str(sys), ' '. - join([sig.str() for sig in sigs.values()])) - nav.fout.write(txt) - nav.fout.write("\n") - - nav.fout.write("Selected signals\n") - for sys, tmp in rnx.sig_tab.items(): - txt = "{:7s} ".format(sys2str(sys)) - for _, sigs in tmp.items(): - txt += "{} ".format(' '.join([sig.str() for sig in sigs])) - nav.fout.write(txt+"\n") - nav.fout.write("\n") - - # Get grid location - # - pos = ecef2pos(rnx.pos) - inet = cs.find_grid_index(pos) - - if l6_mode == 1: - fc = open(l6file, 'rb') - if not fc: - nav.fout.write("ERROR: cannot open L6 message file {}!" - .format(l6file)) - sys_exit(-1) + if l6_mode == 1: # from log file + cs.decode_l6msg(fc.read(250), 0) + if cs.fcnt == 5: # end of sub-frame + cs.week = week + cs.decode_cssr(cs.buff, 0) else: - dtype = [('wn', 'int'), ('tow', 'int'), ('prn', 'int'), - ('type', 'int'), ('len', 'int'), ('nav', 'S500')] - v = np.genfromtxt(file_l6, dtype=dtype) + msg = decode_msg(v, tow, l6_ch, prn_p1) - # Skip epoch until start time - # - obs = rnx.decode_obs() - while time > obs.t and obs.t.time != 0: - obs = rnx.decode_obs() + if msg is not None: + cs.decode_l6msg(msg, 0) + if cs.fcnt == 5: # end of sub-frame + cs.decode_cssr(bytes(cs.buff), 0) - for ne in range(nep): + if prn_p2 > 0: + msg2 = decode_msg(v, tow, l6_ch, prn_p2) - week, tow = time2gpst(obs.t) + if msg2 is not None: + cs_.decode_l6msg(msg2, 0) + if cs_.fcnt == 5: # end of sub-frame + cs_.decode_cssr(bytes(cs_.buff), 0) + cs.merge_cssr(cs_) - if l6_mode == 1: - cs.decode_l6msg(fc.read(250), 0) - if cs.fcnt == 5: # end of sub-frame - cs.week = week - cs.decode_cssr(cs.buff, 0) - else: - vi = v[(v['tow'] == tow) & (v['type'] == l6_ch)] - vi_p1 = vi[vi['prn'] == prn_p1] - vi_p2 = vi[vi['prn'] == prn_p2] - if len(vi_p1) > 0: - cs.decode_l6msg(unhexlify(vi_p1['nav'][0]), 0) - if cs.fcnt == 5: # end of sub-frame - cs.decode_cssr(bytes(cs.buff), 0) - if len(vi_p2) > 0: - cs_.decode_l6msg(unhexlify(vi_p2['nav'][0]), 0) - if cs_.fcnt == 5: # end of sub-frame - cs_.decode_cssr(bytes(cs_.buff), 0) - - if ne == 0: - nav.t = deepcopy(obs.t) - t0 = deepcopy(obs.t) - t0.time = t0.time//30*30 - cs.time = obs.t - nav.time_p = t0 - - cstat = cs.chk_stat() - if cstat: - ppprtk.process(obs, cs=cs) - - t[ne] = timediff(nav.t, t0)/60 - - sol = nav.xa[0:3] if nav.smode == 4 else nav.x[0:3] - enu[ne, :] = gn.ecef2enu(pos_ref, sol-xyz_ref) - smode[ne] = nav.smode - nsat[ne, :] = nav.nsat - - nav.fout.write("{} {:14.4f} {:14.4f} {:14.4f} " - "ENU {:7.3f} {:7.3f} {:7.3f}, 2D {:6.3f}, mode {:1d}\n" - .format(time2str(obs.t), - sol[0], sol[1], sol[2], - enu[ne, 0], enu[ne, 1], enu[ne, 2], - np.sqrt(enu[ne, 0]**2+enu[ne, 1]**2), - smode[ne])) - - # Log to standard output - # - stdout.write('\r {} ENU {:7.3f} {:7.3f} {:7.3f}, 2D {:6.3f}, mode {:1d}' - .format(time2str(obs.t), - enu[ne, 0], enu[ne, 1], enu[ne, 2], - np.sqrt(enu[ne, 0]**2+enu[ne, 1]**2), - smode[ne])) - - # Get new epoch, exit after last epoch - # - obs = rnx.decode_obs() - if obs.t.time == 0: - break - - # Send line-break to stdout - # - stdout.write('\n') + cstat = cs.chk_stat() + if cstat: + ppprtk.process(obs, cs=cs) - # Close RINEX observation and CLAS correction file - # - if l6_mode == 1: - fc.close() - rnx.fobs.close() + proc.save_output(obs.t, ne) # save output - # Close output file + # Get new epoch, exit after last epoch # - if nav.fout is not None: - nav.fout.close() - -fig_type = 1 - -if fig_type == 1: - plot_enu(t, enu, smode) -elif fig_type == 2: - plot_enu(t, enu, smode, figtype=fig_type) + obs = rnx.decode_obs() + if obs.t.time == 0: + break -plotFileFormat = 'eps' -plotFileName = '.'.join(('test_ppprtk', plotFileFormat)) +if l6_mode == 1: + fc.close() -plt.savefig(plotFileName, format=plotFileFormat, bbox_inches='tight', dpi=300) -# plt.show() +proc.close() +proc.plot(ttl, fig_type=1) diff --git a/samples/test_ppprtk2.py b/samples/test_ppprtk2.py index ff1718c..b1c7d3c 100644 --- a/samples/test_ppprtk2.py +++ b/samples/test_ppprtk2.py @@ -7,26 +7,33 @@ from cssrlib.cssrlib import cssr import cssrlib.rinex as rn -import cssrlib.gnss as gn -from cssrlib.gnss import rSigRnx, time2str, sys2str, epoch2time +from cssrlib.gnss import rSigRnx, time2str, sys2str, time2doy, Nav, time2gpst +from cssrlib.gnss import epoch2time, load_config, pos2ecef, ecef2pos, timediff +from cssrlib.gnss import ecef2enu from cssrlib.ppprtk import ppprtkpos from cssrlib.peph import atxdec, searchpcv +config = load_config('config_ppprtk.yml') + +ttl = 'test_ppprtk2' ep = [2021, 9, 22, 6, 30, 0] -time = gn.epoch2time(ep) -atxfile = '../data/antex/' -if time > epoch2time([2022, 11, 27, 0, 0, 0]): - atxfile += 'igs20.atx' -else: - atxfile += 'igs14.atx' -navfile = '../data/doy2021-265/SEPT2650.21P' -obsfile = '../data/doy2021-265/SEPT265G.21O' -l6file = '../data/doy2021-265/2021265G.l6' +time = epoch2time(ep) +year = ep[0] +doy = int(time2doy(time)) +ses = chr(ord('A')+ep[3]) + +if time < epoch2time([2022, 11, 27, 0, 0, 0]): + config['atxfile'] = '../data/antex/igs14.atx' + griddef = '../data/clas_grid.def' -xyz_ref = gn.pos2ecef([35.342058098, 139.521986657, 47.5515], True) +navfile = f'../data/doy{year}-{doy:03d}/SEPT{doy:03d}0.{year%100:02d}P' +obsfile = f'../data/doy{year}-{doy:03d}/SEPT{doy:03d}{ses}.{year%100:02d}O' +l6file = f'../data/doy{year}-{doy:03d}/{year:4d}{doy:03d}{ses}.l6' + +xyz_ref = pos2ecef([35.342058098, 139.521986657, 47.5515], True) # Initial position guess # @@ -34,7 +41,7 @@ cs = cssr() cs.monlevel = 1 -cs.week = 2176 # 2021/9/22 +cs.week,_ = time2gpst(time) cs.read_griddef(griddef) nep = 360 @@ -42,12 +49,12 @@ enu = np.zeros((nep, 3)) smode = np.zeros(nep, dtype=int) # rr0 = [-3961951.1326752, 3381198.11019757, 3668916.0417232] # from pntpos -pos_ref = gn.ecef2pos(xyz_ref) +pos_ref = ecef2pos(xyz_ref) # Load ANTEX data for satellites and stations # atx = atxdec() -atx.readpcv(atxfile) +atx.readpcv(config['atxfile']) # Define signals to be processed # @@ -70,7 +77,7 @@ rnx = rn.rnxdec() rnx.setSignals(sigs) -nav = gn.Nav() +nav = Nav() rnx.decode_nav(navfile, nav) if rnx.decode_obsh(obsfile) >= 0: @@ -81,7 +88,7 @@ # Initialize position # - ppprtk = ppprtkpos(nav, rnx.pos, 'test_ppprtk2.log') + ppprtk = ppprtkpos(nav, rnx.pos, logfile=f'{ttl}.log', config=config) nav.armode = 3 # Satellite exclusion @@ -138,7 +145,7 @@ nav.fout.write(txt+"\n") nav.fout.write("\n") - pos = gn.ecef2pos(rr0) + pos = ecef2pos(rr0) inet = cs.find_grid_index(pos) fc = open(l6file, 'rb') @@ -163,7 +170,7 @@ for ne in range(nep): - week, tow = gn.time2gpst(obs.t) + week, tow = time2gpst(obs.t) cs.decode_l6msg(fc.read(250), 0) if cs.fcnt == 5: # end of sub-frame @@ -179,9 +186,9 @@ if cstat: ppprtk.process(obs, cs=cs) - t[ne] = gn.timediff(nav.t, t0) + t[ne] = timediff(nav.t, t0) sol = nav.xa[0:3] if nav.smode == 4 else nav.x[0:3] - enu[ne, :] = gn.ecef2enu(pos_ref, sol-xyz_ref) + enu[ne, :] = ecef2enu(pos_ref, sol-xyz_ref) smode[ne] = nav.smode nav.fout.write("{} {:14.4f} {:14.4f} {:14.4f} {:14.4f} {:14.4f} {:14.4f} {:2d}\n" @@ -231,7 +238,7 @@ plt.legend() plt.grid() -plotFileName = '.'.join(('test_ppprtk2_1', plotFileFormat)) +plotFileName = '.'.join((f'{ttl}_1', plotFileFormat)) plt.savefig(plotFileName, format=plotFileFormat, bbox_inches='tight', dpi=300) # plt.show() @@ -250,7 +257,7 @@ plt.axis('equal') plt.legend() -plotFileName = '.'.join(('test_ppprtk2_2', plotFileFormat)) +plotFileName = '.'.join((f'{ttl}_2', plotFileFormat)) plt.savefig(plotFileName, format=plotFileFormat, bbox_inches='tight', dpi=300) # plt.show() diff --git a/samples/test_sbas.py b/samples/test_sbas.py index 6bf4cd1..484499d 100644 --- a/samples/test_sbas.py +++ b/samples/test_sbas.py @@ -2,43 +2,60 @@ static test for SBAS (L1 or DFMC) """ from binascii import unhexlify -from copy import deepcopy -import matplotlib.pyplot as plt -import matplotlib.dates as md import numpy as np -from sys import stdout from sys import exit as sys_exit -from cssrlib.gnss import ecef2pos, Nav, ecef2enu -from cssrlib.gnss import time2gpst, time2doy, time2str, timediff, epoch2time -from cssrlib.gnss import rSigRnx, sys2str, uIonoModel -from cssrlib.peph import atxdec, searchpcv +from cssrlib.gnss import Nav, load_config, uIonoModel +from cssrlib.gnss import time2doy, timediff, epoch2time from cssrlib.pntpos import stdpos from cssrlib.sbas import sbasDec from cssrlib.rinex import rnxdec from cssrlib.cssr_pvs import decode_sinca_line +from cssrlib.utils import process + + +def decode_msg(v, tow, prn_ref, sbas_type=0): + """ find valid correction message """ + + if len(prn_ref) == 1: + vi = v[(v['tow'] == tow) & (v['prn'] == prn_ref)] + else: + vi = v[(v['tow'] == tow) & (v['prn'] >= prn_ref[0]) & + (v['prn'] <= prn_ref[1])] + if sbas_type == 0: # L1 + vi = vi[vi['type'] <= 28] + else: # DFMC L5 + vi = vi[(vi['type'] == 31) | (vi['type'] == 32) | + ((vi['type'] >= 34) & (vi['type'] <= 37))] + if len(vi) > 0: + msg = {} + for vi_ in vi: + msg[vi_['prn']] = unhexlify(vi_['nav']) + else: + msg = None + + return msg # Select test case # -dataset = 4 +nep = 3600 +ttl = 'test_sbas' # title for case +dataset = 1 +navfile = None + +config = load_config('config.yml') # Start epoch and number of epochs # if dataset == 1: # MSAS, L1 SBAS - ep = [2025, 2, 15, 17, 0, 0] - navfile = '../data/doy2025-046/046r_rnx.nav' - obsfile = '../data/doy2025-046/046r_rnx.obs' # mosaic-X5 - file_sbas = '../data/doy2025-046/046r_sbas.txt' + ep = [2025, 8, 21, 7, 0, 0] xyz_ref = [-3962108.6836, 3381309.5672, 3668678.6720] prn_ref = [137] # satellite PRN for SBAS correction sbas_type = 0 # L1: 0, L5: 1 nf = 1 elif dataset == 2: # QZSS, L5 DFMC - ep = [2025, 2, 15, 17, 0, 0] - navfile = '../data/doy2025-046/046r_rnx.nav' - obsfile = '../data/doy2025-046/046r_rnx.obs' # mosaic-X5 - file_sbas = '../data/doy2025-046/046r_sbas.txt' + ep = [2025, 8, 21, 7, 0, 0] xyz_ref = [-3962108.6836, 3381309.5672, 3668678.6720] # prn_ref = [193, 202] # satellite PRN for SBAS correction prn_ref = [199] @@ -56,11 +73,7 @@ nf = 2 elif dataset == 4: # SouthPAN L5 DFMC (SIS) - ep = [2025, 8, 21, 7, 0, 0] - navfile = '../data/doy2025-233/233h_rnx.nav' - obsfile = '../data/doy2025-233/233h_rnx.obs' # SEPT MOSAIC-X5 - file_sbas = '../data/doy2025-233/233h_sbas.txt' xyz_ref = [-3962108.6836, 3381309.5672, 3668678.6720] # Kamakura prn_ref = [122] sbas_type = 1 # L1: 0, L5: 1 @@ -69,298 +82,99 @@ time = epoch2time(ep) year = ep[0] doy = int(time2doy(time)) +ses = chr(ord('a')+ep[3]) -nep = 900*4 -# nep = 360 - +if navfile is None: + navfile = f'../data/doy{year}-{doy:03d}/{doy:03d}{ses}_rnx.nav' + obsfile = f'../data/doy{year}-{doy:03d}/{doy:03d}{ses}_rnx.obs' + file_sbas = f'../data/doy{year}-{doy:03d}/{doy:03d}{ses}_sbas.txt' -pos_ref = ecef2pos(xyz_ref) - -# Define signals to be processed -# -sigs = [] -if sbas_type == 0: # single frequency L1 SBAS - nf = 1 - gnss = "G" - if 'G' in gnss: - sigs.extend([rSigRnx("GC1C"), rSigRnx("GL1C"), rSigRnx("GS1C")]) +if sbas_type == 0: # L1 SBAS + sig_t = {'G': ['1C']} -else: # dual frequency - nf = 2 - gnss = "GE" - if 'G' in gnss: - sigs.extend([rSigRnx("GC1C"), rSigRnx("GC5Q"), - rSigRnx("GL1C"), rSigRnx("GL5Q"), - rSigRnx("GS1C"), rSigRnx("GS5Q")]) - if 'E' in gnss: - sigs.extend([rSigRnx("EC1C"), rSigRnx("EC5Q"), - rSigRnx("EL1C"), rSigRnx("EL5Q"), - rSigRnx("ES1C"), rSigRnx("ES5Q")]) - if 'J' in gnss: - sigs.extend([rSigRnx("JC1C"), rSigRnx("JC5Q"), - rSigRnx("JL1C"), rSigRnx("JL5Q"), - rSigRnx("JS1C"), rSigRnx("JS5Q")]) - if 'S' in gnss: - sigs.extend([rSigRnx("SC1C"), rSigRnx("SC5Q"), - rSigRnx("SL1C"), rSigRnx("SL5Q"), - rSigRnx("SS1C"), rSigRnx("SS5Q")]) +else: # DFMC + sig_t = {'G': ['1C', '5Q'], 'E': ['1C', '5Q']} rnx = rnxdec() -rnx.setSignals(sigs) - -nav = Nav(nf=nf) +nav = Nav() +proc = process(nav, rnx, config, nep=nep, xyz_ref=xyz_ref) -# Positioning mode -# 0:static, 1:kinematic +# Define signals to be processed # -nav.pmode = 1 +sigs, nav.nf = proc.init_sig(sig_t) +rnx.setSignals(sigs) # Decode RINEX NAV data # nav = rnx.decode_nav(navfile, nav) -cs = sbasDec('test_sbas_cs.log') -cs.monlevel = 2 - -# Load ANTEX data for satellites and stations -# -atx = atxdec() -atx.readpcv('../data/antex/igs20.atx') +cs = sbasDec('test_sbas_cs.log', nf=nf) +cs.monlevel = config['cs']['monlevel'] -# Initialize data structures for results -# -t = np.zeros(nep) -enu = np.ones((nep, 3))*np.nan -sol = np.zeros((nep, 4)) -ztd = np.zeros((nep, 1)) -smode = np.zeros(nep, dtype=int) -nsat = np.zeros(nep, dtype=int) - -# Logging level -# -nav.monlevel = 1 # TODO: enabled for testing! +rnx.decode_obsh(obsfile) # Load RINEX OBS file header -# Load RINEX OBS file header +# Initialize position # -if rnx.decode_obsh(obsfile) >= 0: +std = stdpos(nav, rnx.pos, f'{ttl}.log') +std.monlevel = 1 +std.ionoModel = uIonoModel.SBAS - # Auto-substitute signals - # - rnx.autoSubstituteSignals() +proc.prepare_signal(obsfile) - # Initialize position - # - std = stdpos(nav, rnx.pos, 'test_sbas.log') - std.monlevel = 1 - nav.elmin = np.deg2rad(5.0) - - std.ionoModel = uIonoModel.SBAS +nav.rmode = 2 if nav.nf == 2 else 0 # L1/L5 iono-free combination - nav.nf = nf - nav.rmode = 2 if nav.nf == 2 else 0 # L1/L5 iono-free combination - nav.csmooth = True +# Skip epochs until start time +# +obs = rnx.decode_obs() +while time > obs.t and obs.t.time != 0: + obs = rnx.decode_obs() - # Get equipment information - # - nav.fout.write("FileName: {}\n".format(obsfile)) - nav.fout.write("Start : {}\n".format(time2str(rnx.ts))) - if rnx.te is not None: - nav.fout.write("End : {}\n".format(time2str(rnx.te))) - nav.fout.write("Receiver: {}\n".format(rnx.rcv)) - nav.fout.write("Antenna : {}\n".format(rnx.ant)) - nav.fout.write("\n") - - if 'UNKNOWN' in rnx.ant or rnx.ant.strip() == "": - nav.fout.write("ERROR: missing antenna type in RINEX OBS header!\n") - - # Set PCO/PCV information - # - nav.sat_ant = atx.pcvs - nav.rcv_ant = searchpcv(atx.pcvr, rnx.ant, rnx.ts) - if nav.rcv_ant is None: - nav.fout.write("ERROR: missing antenna type <{}> in ANTEX file!\n" - .format(rnx.ant)) +if 'sbas' in file_sbas: # SIS + dtype = [('wn', 'int'), ('tow', 'float'), ('prn', 'int'), + ('type', 'int'), ('marker', 'S2'), ('nav', 'S124')] + v = np.genfromtxt(file_sbas, dtype=dtype) +elif 'DAS' in file_sbas: # DAS + fc = open(file_sbas, 'rt') +else: + print("ERROR: unknown file format for correction data") + sys_exit(1) + +# Loop over number of epoch from file start +# +for ne in range(nep): + _, tow = cs.set_time(obs.t) # set time for reference - # Print available signals + # Set initial epoch # - nav.fout.write("Available signals\n") - for sys, sigs in rnx.sig_map.items(): - txt = "{:7s} {}\n".format(sys2str(sys), - ' '.join([s.str() for s in sigs.values()])) - nav.fout.write(txt) - nav.fout.write("\n") - - nav.fout.write("Selected signals\n") - for sys, tmp in rnx.sig_tab.items(): - txt = "{:7s} ".format(sys2str(sys)) - for _, sigs in tmp.items(): - txt += "{} ".format(' '.join([sig.str() for sig in sigs])) - nav.fout.write(txt+"\n") - nav.fout.write("\n") - - # Skip epochs until start time - # - obs = rnx.decode_obs() - while time > obs.t and obs.t.time != 0: - obs = rnx.decode_obs() + if ne == 0: + proc.init_time(obs.t) if 'sbas' in file_sbas: # SIS - dtype = [('wn', 'int'), ('tow', 'float'), ('prn', 'int'), - ('type', 'int'), ('marker', 'S2'), ('nav', 'S124')] - v = np.genfromtxt(file_sbas, dtype=dtype) - elif 'DAS' in file_sbas: # DAS - fc = open(file_sbas, 'rt') - else: - print("ERROR: unknown file format for correction data") - sys_exit(1) - - # Loop over number of epoch from file start + msgs = decode_msg(v, tow, prn_ref, sbas_type) + if msgs is not None: + for prn, msg in msgs.items(): + cs.decode_cssr(msg, 0, src=sbas_type, prn=prn) + + else: # DAS + for line in fc: + tc, buff = decode_sinca_line(line) + cs.decode_cssr(buff, 0, src=sbas_type) + if timediff(obs.t, tc) >= 0.0: + break + + # Call standard positioning module with SBAS corrections # - for ne in range(nep): - - week, tow = time2gpst(obs.t) - cs.week = week - cs.tow0 = tow//86400*86400 - cs.time = obs.t - - # Set initial epoch - # - if ne == 0: - nav.t = deepcopy(obs.t) - t0 = deepcopy(obs.t) - t0.time = t0.time//30*30 - nav.time_p = t0 - - if 'sbas' in file_sbas: # SIS - - if len(prn_ref) == 1: - vi = v[(v['tow'] == tow) & (v['prn'] == prn_ref)] - else: - vi = v[(v['tow'] == tow) & (v['prn'] >= prn_ref[0]) & - (v['prn'] <= prn_ref[1])] - if sbas_type == 0: # L1 - vi = vi[vi['type'] <= 28] - else: # DFMC L5 - vi = vi[(vi['type'] == 31) | (vi['type'] == 32) | - ((vi['type'] >= 34) & (vi['type'] <= 37))] - if len(vi) > 0: - for vi_ in vi: - buff = unhexlify(vi_['nav']) - cs.decode_cssr(buff, 0, src=sbas_type, prn=vi_['prn']) - - else: # DAS - for line in fc: - tc, buff = decode_sinca_line(line) - cs.decode_cssr(buff, 0, src=sbas_type) - if timediff(obs.t, tc) >= 0.0: - break - - # cs.check_validity(obs.t) - - # Call PPP module with PVS corrections - # + if (cs.lc[0].cstat & 0x6) == 0x6: # wait for orbit/clock correction std.process(obs, cs=cs) - # Save output - # - t[ne] = timediff(nav.t, t0)/86400.0 - - sol = nav.xa[0:3] if nav.smode == 4 else nav.x[0:3] - enu[ne, :] = ecef2enu(pos_ref, sol-xyz_ref) - - # ztd[ne] = nav.xa[std.IT(nav.na)] \ - # if nav.smode == 4 else nav.x[std.IT(nav.na)] - smode[ne] = nav.smode - # nsat[ne] = std.nsat - - nav.fout.write("{} {:14.4f} {:14.4f} {:14.4f} " - "ENU {:7.3f} {:7.3f} {:7.3f}, 2D {:6.3f}, mode {:1d}, " - "nsat {:1d}\n" - .format(time2str(obs.t), - sol[0], sol[1], sol[2], - enu[ne, 0], enu[ne, 1], enu[ne, 2], - np.sqrt(enu[ne, 0]**2+enu[ne, 1]**2), - smode[ne], std.nsat)) - - # Log to standard output - # - stdout.write("\r {} ENU {:7.3f} {:7.3f} {:7.3f}, 2D {:6.3f}, " - "mode {:1d}, nsat {:1d}" - .format(time2str(obs.t), - enu[ne, 0], enu[ne, 1], enu[ne, 2], - np.sqrt(enu[ne, 0]**2+enu[ne, 1]**2), - smode[ne], std.nsat)) - - # Get new epoch, exit after last epoch - # - obs = rnx.decode_obs() - if obs.t.time == 0: - break - - # Send line-break to stdout - # - stdout.write('\n') - - # Close RINEX observation file - # - rnx.fobs.close() + proc.save_output(obs.t, ne) # save output - # Close output file + # Get new epoch, exit after last epoch # - if nav.fout is not None: - nav.fout.close() - -fig_type = 1 -ylim_h = 2.0 -ylim_v = 3.0 - -idx2 = np.where(smode == 2)[0] -idx1 = np.where(smode == 1)[0] -idx0 = np.where(smode == 0)[0] - -fig = plt.figure(figsize=[7, 9]) -fig.set_rasterized(True) - -fmt = '%H:%M' - -if fig_type == 1: - - lbl_t = ['East [m]', 'North [m]', 'Up [m]'] - - for k in range(3): - ylim = ylim_h if k < 2 else ylim_v - - plt.subplot(3, 1, k+1) - plt.plot(t[idx0], enu[idx0, k], 'r.', label='none') - plt.plot(t[idx2], enu[idx2, k], 'y.', label='SBAS/DGPS') - plt.plot(t[idx1], enu[idx1, k], 'g.', label='standalone') - - plt.ylabel(lbl_t[k]) - plt.grid() - plt.ylim([-ylim, ylim]) - plt.gca().xaxis.set_major_formatter(md.DateFormatter(fmt)) - - plt.xlabel('Time [HH:MM]') - plt.legend() - -elif fig_type == 2: - - ax = fig.add_subplot(111) - - plt.plot(enu[idx0, 0], enu[idx0, 1], 'r.', label='none') - plt.plot(enu[idx2, 0], enu[idx2, 1], 'y.', label='SBAS/DGPS') - plt.plot(enu[idx1, 0], enu[idx1, 1], 'g.', label='standalone') - - plt.xlabel('Easting [m]') - plt.ylabel('Northing [m]') - plt.grid() - plt.axis('equal') - plt.legend() - # ax.set(xlim=(-ylim, ylim), ylim=(-ylim, ylim)) - -plotFileFormat = 'png' -plotFileName = '.'.join(('test_sbas', plotFileFormat)) + obs = rnx.decode_obs() + if obs.t.time == 0: + break -plt.savefig(plotFileName, format=plotFileFormat, bbox_inches='tight', dpi=300) -# plt.show() +proc.close() +proc.plot(ttl, fig_type=1, ylim=2, ylim_v=4)