はじめに

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