はじめに
GNSSのRAWデータ(観測値)の共通フォーマットであるRINEX形式のファイルから値を読んでプロットします。
頑張って自分で読み書きするプログラムを書くのは好きですが（実際に以前に行いましたが）、今回はOSSを利用します。便利な世の中になったものだ。
内容
準備
georinex を使います。普通にインストールできました。
また、読み込むのに利用するRINEXファイルを用意します。自分で受信機を持って外で取得しても良いですが、電子基準点の値をダウンロードすることもできます。ここでは、こちらから取得したものを利用しました。
読んで描画
読み込みます。わりと時間がかかります。
データはload() で読み込めます。RINEXのversion はヘッダに書かれている内容を自動で読みます。
infile = './3019148c.23o'
rnxobs = gr.load(infile)
ここで読み込んだのはxarray というものらしく、下記のように表示されます。
In [3]: type(rnxobs)
Out[3]: xarray.core.dataset.Dataset
In [4]: print(rnxobs)
<xarray.Dataset>
Dimensions: (time: 119, sv: 32)
Coordinates:
* time (time) datetime64[ns] 2023-05-28T02:00:00 ... 2023-05-28T02:59:00
* sv (sv) <U3 'E10' 'E12' 'E14' 'E21' 'E24' ... 'R16' 'R17' 'R18' 'R19'
Data variables: (12/44)
C1C (time, sv) float64 nan nan nan nan ... nan 2.171e+07 2.016e+07
L1C (time, sv) float64 nan nan nan nan ... nan 1.159e+08 1.078e+08
D1C (time, sv) float64 nan nan nan nan ... nan -2.064e+03 -365.7
S1C (time, sv) float64 nan nan nan nan nan ... nan 46.2 nan 45.9 40.5
C1W (time, sv) float64 nan nan nan nan nan nan ... nan nan nan nan nan
L1W (time, sv) float64 nan nan nan nan nan nan ... nan nan nan nan nan
... ...
D7X (time, sv) float64 nan 1.826e+03 -3.091e+03 ... nan nan nan
S7X (time, sv) float64 nan 39.5 43.4 39.2 nan ... nan nan nan nan nan
C8X (time, sv) float64 nan 2.676e+07 2.011e+07 ... nan nan nan
L8X (time, sv) float64 nan 1.064e+08 7.995e+07 ... nan nan nan
D8X (time, sv) float64 nan 1.803e+03 -3.052e+03 ... nan nan nan
S8X (time, sv) float64 nan 40.0 43.6 39.0 nan ... nan nan nan nan nan
Attributes:
version: 3.02
interval: 30.0
rinextype: obs
fast_processing: 0
time_system: GPS
filename: 3019148c.23o
position: [-3941956.1858, 3368145.4203, 3702217.4723]
Xarray 読み解く
解釈が合っているか分かりませんが、DataFrame の index は一種類の値を参照するキーとしますが、このCoordinate は同じようなモノを複数の要素で参照することができるようになるようです。時刻でも参照したいし、衛星番号でも参照したい、という要望に応えられるようです。
時刻は３０秒周期で１時間のものなので、120個の要素があります。datetime で保持してくれているのはありがたいですね。
In [5]: rnxobs.time
Out[5]:
<xarray.DataArray 'time' (time: 119)>
array(['2023-05-28T02:00:00.000000000', '2023-05-28T02:00:30.000000000',
'2023-05-28T02:01:00.000000000', '2023-05-28T02:01:30.000000000',
...
'2023-05-28T02:59:00.000000000'], dtype='datetime64[ns]')
Coordinates:
* time (time) datetime64[ns] 2023-05-28T02:00:00 ... 2023-05-28T02:59:00
In [6]: len( rnxobs.time)
Out[6]: 119
あれ、120個のはずが119個とはおかしいな？と思い、下記で調べてみました。実際に、59分30秒のデータは入ってませんでした。地理院サーバ側のバグかな。
$ cat 3019148c.23o | grep \> | tail -3
> 2023 05 28 02 58 0.0000000 0 24
> 2023 05 28 02 58 30.0000000 0 24
> 2023 05 28 02 59 0.0000000 0 24
衛星は、この時刻では32機観測されたようです。
In [9]: rnxobs.sv
Out[9]:
<xarray.DataArray 'sv' (sv: 32)>
array(['E10', 'E12', 'E14', 'E21', 'E24', 'E26', 'E31', 'E33', 'G01', 'G02',
'G03', 'G07', 'G08', 'G14', 'G17', 'G21', 'G27', 'G30', 'J02', 'J03',
'J04', 'J07', 'R02', 'R03', 'R04', 'R09', 'R10', 'R15', 'R16', 'R17',
'R18', 'R19'], dtype='<U3')
Coordinates:
* sv (sv) <U3 'E10' 'E12' 'E14' 'E21' 'E24' ... 'R16' 'R17' 'R18' 'R19'
観測値は xarray の Data Variables で24種類あります。
In [12]: rnxobs.data_vars
Out[12]:
Data variables:
C1C (time, sv) float64 nan nan nan nan ... nan 2.171e+07 2.016e+07
L1C (time, sv) float64 nan nan nan nan ... nan 1.159e+08 1.078e+08
D1C (time, sv) float64 nan nan nan nan ... nan -2.064e+03 -365.7
S1C (time, sv) float64 nan nan nan nan nan ... nan 46.2 nan 45.9 40.5
C1W (time, sv) float64 nan nan nan nan nan nan ... nan nan nan nan nan
信号の種類、周波数の種類などで３文字のコードになっていますが、これはRINEX ファイルのヘッダで定義されている対応関係をそのまま値として格納しているようです。
$ head -40 3019148c.23o
3.02 OBSERVATION DATA M (MIXED) RINEX VERSION / TYPE
BINEX2RINEX 2.10 GSI, JAPAN 20230528 03:43:06UTCPGM / RUN BY / DATE
3019 MARKER NAME
GEODETIC MARKER TYPE
GSI, JAPAN GEOSPATIAL INFORMATION AUTHORITY OF JAPAOBSERVER / AGENCY
00000 TPS NETG5 5.2.6,08/May/2020 REC # / TYPE / VERS
TRM159900.00 GSI ANT # / TYPE
-3941956.1858 3368145.4203 3702217.4723 APPROX POSITION XYZ
0.0000 0.0000 0.0000 ANTENNA: DELTA H/E/N
G 20 C1C L1C D1C S1C C1W L1W D1W S1W C2X L2X D2X S2X C2W SYS / # / OBS TYPES
L2W D2W S2W C5X L5X D5X S5X SYS / # / OBS TYPES
R 16 C1C L1C D1C S1C C1P L1P D1P S1P C2C L2C D2C S2C C2P SYS / # / OBS TYPES
L2P D2P S2P SYS / # / OBS TYPES
E 16 C1X L1X D1X S1X C5X L5X D5X S5X C7X L7X D7X S7X C8X SYS / # / OBS TYPES
L8X D8X S8X SYS / # / OBS TYPES
J 12 C1C L1C D1C S1C C2X L2X D2X S2X C5X L5X D5X S5X SYS / # / OBS TYPES
2023 5 28 2 0 0.0000000 GPS TIME OF FIRST OBS
G L1C 0.00000 10 G01 G02 G03 G07 G08 G14 G17 G21 G27 G30 SYS / PHASE SHIFT
G L1W -0.25000 10 G01 G02 G03 G07 G08 G14 G17 G21 G27 G30 SYS / PHASE SHIFT
G L2X 0.25000 8 G01 G03 G07 G08 G14 G17 G27 G30 SYS / PHASE SHIFT
G L2W 0.00000 10 G01 G02 G03 G07 G08 G14 G17 G21 G27 G30 SYS / PHASE SHIFT
G L5X 0.00000 6 G01 G03 G08 G14 G27 G30 SYS / PHASE SHIFT
R L1C 0.00000 10 R02 R03 R04 R09 R10 R15 R16 R17 R18 R19 SYS / PHASE SHIFT
R L1P -0.25000 10 R02 R03 R04 R09 R10 R15 R16 R17 R18 R19 SYS / PHASE SHIFT
R L2C 0.00000 9 R02 R03 R04 R09 R15 R16 R17 R18 R19 SYS / PHASE SHIFT
R L2P -0.25000 9 R02 R03 R04 R09 R15 R16 R17 R18 R19 SYS / PHASE SHIFT
J L1C 0.00000 4 J02 J03 J04 J07 SYS / PHASE SHIFT
J L2X 0.25000 4 J02 J03 J04 J07 SYS / PHASE SHIFT
J L5X 0.00000 4 J02 J03 J04 J07 SYS / PHASE SHIFT
E L1X 0.00000 8 E10 E12 E14 E21 E24 E26 E31 E33 SYS / PHASE SHIFT
E L5X 0.00000 8 E10 E12 E14 E21 E24 E26 E31 E33 SYS / PHASE SHIFT
E L7X 0.00000 8 E10 E12 E14 E21 E24 E26 E31 E33 SYS / PHASE SHIFT
E L8X 0.00000 8 E10 E12 E14 E21 E24 E26 E31 E33 SYS / PHASE SHIFT
10 R02 -4 R03 5 R04 6 R09 -2 R10 -7 R15 0 R16 -1 R17 4 GLONASS SLOT / FRQ #
R18 -3 R19 3 GLONASS SLOT / FRQ #
GLONASS COD/PHS/BIS
30.000 INTERVAL
18 LEAP SECONDS
smtt off - ms jumps in time tag (smooth phase and range) COMMENT
END OF HEADER
Xarray の Data Variables の名前を取得すには、下記でできました。
In [18]: print([i for i in rnxobs.data_vars])
['C1C', 'L1C', 'D1C', 'S1C', 'C1W', 'L1W', 'D1W', 'S1W', 'C2X', 'L2X', 'D2X', 'S2X', 'C2W', 'L2W', 'D2W', 'S2W', 'C5X', 'L5X', 'D5X', 'S5X', 'C1P', 'L1P', 'D1P', 'S1P', 'C2C', 'L2C', 'D2C', 'S2C', 'C2P', 'L2P', 'D2P', 'S2P', 'C1X', 'L1X', 'D1X', 'S1X', 'C7X', 'L7X', 'D7X', 'S7X', 'C8X', 'L8X', 'D8X', 'S8X']
観測値のプロット
ではまず疑似距離をプロットしてみます。
satname = 'G17'
fig, axes = plt.subplots(4,1, figsize=(12,8), sharex=True)
axes[0].plot(rnxobs.time, rnxobs['C1C'].sel(sv=satname), label=satname)
axes[1].plot(rnxobs.time, rnxobs['L1C'].sel(sv=satname), label=satname)
axes[2].plot(rnxobs.time, rnxobs['D1C'].sel(sv=satname), label=satname)
axes[3].plot(rnxobs.time, rnxobs['S1C'].sel(sv=satname), label=satname)
axes[0].set_title('Pseudo range[m]')
axes[0].set_ylabel(r'$\rho$ [m]')
axes[1].set_title('carrier phase')
axes[1].set_ylabel(r'$\phi$ [cycle]')
axes[2].set_title('Doppler frequency')
axes[2].set_ylabel('D[Hz]')
axes[3].set_title('C/N0')
axes[3].set_ylabel('S [dB]')
for ax in axes:
ax.grid(True)
axes[3].set_xlabel('GPST')
plt.tight_layout()
とすると、こんな感じになります。
観測値の足し算引き算
では次に、観測値を足したり引いたりしてみましょう。そのためにL1信号の波長を求めておきます。これは光の早さと周波数（GPSの仕様）から求められます。L1信号の波長、19cm になった。
In [34]: CLIGHT = 299792458.0
...: L1_FREQ=1.57542e9
...: wlen = CLIGHT / L1_FREQ
...: print(wlen)
0.19029367279836487
で、こんな風に足し算引き算してみます。
疑似距離と搬送波位相の差を波数（位相、サイクル数）の単位で計算してみました。実際には伝搬遅延の影響が違ったり、受信機ノイズの影響があったりするので、まぁ、ざっくりとして見積もりです。
cp_l1 = rnxobs['L1C'].sel(sv=satname)
pr_l1 = rnxobs['C1C'].sel(sv=satname)
cp_pr_l1 = cp_l1 - pr_l1 / wlen
確認すると、diff を取ったときの最初のデータの時刻は 2:00:30
In [69]: cp_l1.diff('time') * (1/30.0)
Out[69]:
<xarray.DataArray 'L1C' (time: 118)>
array([-2875.11366667, -2865.25266667, -2855.2751 , -2845.1182 ,
...
-1188.323 , -1169.7848 ])
Coordinates:
* time (time) datetime64[ns] 2023-05-28T02:00:30 ... 2023-05-28T02:59:00
sv <U3 'G14'
一方で、doppler の方はデータは 2:00:00 からある。
In [70]: dp_l1
Out[70]:
<xarray.DataArray 'D1C' (time: 119)>
array([2880.188, 2870.113, 2860.43 , 2850.223, 2839.914, 2829.832,
...
1234.434, 1216.18 , 1197.684, 1179.18 , 1160.547])
Coordinates:
* time (time) datetime64[ns] 2023-05-28T02:00:00 ... 2023-05-28T02:59:00
sv <U3 'G14'
これら、データ数が少ないのですが、引き算できました。
In [71]: dp_l1 + cp_l1.diff('time') * (1/30.0)
Out[71]:
<xarray.DataArray (time: 118)>
array([-5.00066667, -4.82266667, -5.0521 , -5.2042 , -5.0474 ,
...
-9.0376 , -9.143 , -9.2378 ])
Coordinates:
* time (time) datetime64[ns] 2023-05-28T02:00:30 ... 2023-05-28T02:59:00
sv <U3 'G14'
値を見てみると、2:00:30 の値を選択して引き算を実行してくれているのが確認できました。素晴らしい。
In [76]: 2870.113 -2875.11366667
Out[76]: -5.000666670000101
こんな感じにプロットできます。案外、搬送波位相(ADR, accumurated delta range)の初期バイアスに近い値が求められているのかな。いい感じだ。^^)g
まとめ
RINEXの観測値を読むことができました。また、それをプロットすることもできました。来週も生き延びられるかな。。。デコーダ作らずに済むので少しだけ長生きできそうだ。
このgeorinex は便利そうなので、もう少し見てみます。使うにはxarrayをお勉強しなくては。
とまぁ、どこかでGPSの知識を生かして研究開発、ビジネスできないかなぁ。ご連絡お待ちしています。なんて。
(2023/05/28)
追記：
- コードはここに置く予定。https://github.com/torupati/gnssraw.git