MacでRTKLIBを使ってGNSS-PPK

# はじめに

[RTKLIB](http://www.rtklib.com)を使ってGNSSの観測結果の後処理（PPK）をすることで高精度な測量ができるそうで，とても興味があるのでやってみました．
ただ，RTKLIBはWindowsだとGUIアプリとして簡単に使えるのですが，MacやLinuxの場合だとCLIのみで，さらに自分でコンパイルしないといけないようなので，記録として残しておきます．

# 事前準備

当方Macに明るくないという致命的な問題がありますが，恐らく事前準備としてgccのインストールが必要です．というのも，デフォルト状態のMacだとgccはインストールされていないor`gcc`コマンドがあっても，単なるclangのエイリアスになっているようです．
必要に応じて
```
brew install gcc
```
などして，インストールしておくと良いと思います．　　
当方環境（M1 Mac Monterey）では，正規のgccは`gcc-11`というコマンドになっていました．

# コンパイル

まずはソースコードをダウンロードして展開します．今回は，`~/Desktop`を作業ディレクトリとします．

```
wget https://github.com/tomojitakasu/RTKLIB/archive/refs/tags/v2.4.3-b34.zip
unzip v2.4.3-b34.zip
```

あとは，PPKに必要なものだけコンパイルします．

まずは，何に使うのかよくわかりませんが，`iers`なるものをコンパイルします．

```
cd RTKLIB-2.4.3-b34/lib/iers/gcc
```

ここで，`makefile`をテキストエディタで開き，`F77 = gfortran`と書かれた部分を`F77 = gfortran-11`と書き換えます．  
あとは，コンパイルです．

```
make
```

次にコンパイルするのは，PPKを行うCLIプログラム本体のrnx2rtkpです．

```
cd ../../../app/consapp/rnx2rtkp/gcc
```

ここで`makefile`を下記の通りに書き換えます．

1. ２行目あたりに`CC = gcc-11`を追加
2.  `LDLIBS = ../../../../lib/iers/gcc/iers.a -lgfortran -lm -lrt`と書かれた行から`-lrt`を消す．
3. インストールディレクトリを`BINDIR = `の部分に書く．今回は，`BINDIR = /Users/hogehoge/Desktop/RTKLIB-2.4.3-b34/bin`としました．

`-lrt`を消してよいのかよくわかりませんが，stackoverflowにとりま消してみたらうまくいったという事例があったので，そのようにしてみます．
[https://stackoverflow.com/questions/1505402/library-not-found-for-lrt-with-qtcreator-mac-os](https://stackoverflow.com/questions/1505402/library-not-found-for-lrt-with-qtcreator-mac-os)

あとは，コンパイルしてインストールです．

```
make
make install
```

これにて，インストールディレクトリに`rnx2rtkp`というバイナリができていればおっけーです．適当に，インストールディレクトリをパスに追加しておくと便利です．

```
export PATH="/Users/hogehoge/Desktop/RTKLIB-2.4.3-b34/bin:$PATH"
```

# PPKやってみる

さて，本題のPPKをやってみます．今回は国土地理院電子基準点を基準局としてPPKをやってみます．  
ということで，移動局のGNSS測量でげっとしてきたRINEXファイルと同時刻の電子基準点データを用意します．電子基準点の場合はPCV補正データが使えるので，そちらもダウンロードしておきます．

ここでは，

* `rover.o`：移動体のRINEXファイル
* `base.22o`：電子基準点のRINEXファイル
* `base.nav.22*`：電子基準点のナビゲーションファイル（衛星軌道情報）
* `GSI_PCV.TXT`：PCV補正データ

とします．電子基準点データのダウンロードについては，ググるなりしてうまくやってください．

あとは，防災科研の内山大先生の教えに従ってPPKのパラメータと上記のファイルを`rnx2rtkp`に食わせってあげます．
[https://hdtopography.github.io/learning/book/GNSS/RTKLIB_uchiyama_20180729v3s.pdf](https://hdtopography.github.io/learning/book/GNSS/RTKLIB_uchiyama_20180729v3s.pdf)

PPKの設定は，愚直にコマンドラインに書いてもよいのですが，設定ファイルを準備したほうがわかりやすいと思うので，下記の設定ファイル`conf.txt`を用意します．

```conf.txt
# RTKNAVI options (2013/03/01 10:41:04, v.2.4.2)

pos1-posmode       =kinematic     # (0:single,1:dgps,2:kinematic,3:static,4:movingbase,5:fixed,6:ppp-kine,7:ppp-static)
pos1-frequency     =2      # (1:l1,2:l1+l2,3:l1+l2+l5)
pos1-soltype       =forward    # (0:forward,1:backward,2:combined)
pos1-elmask        =10         # (deg)
pos1-snrmask_r     =off        # (0:off,1:on)
pos1-snrmask_b     =off        # (0:off,1:on)
pos1-snrmask_L1    =35,35,35,35,35,35,35,35,35
pos1-snrmask_L2    =35,35,35,35,35,35,35,35,35
pos1-snrmask_L5    =0,0,0,0,0,0,0,0,0
pos1-dynamics      =off        # (0:off,1:on)
pos1-tidecorr      =Solid/OTL        # (0:off,1:on)
pos1-ionoopt       =brdc       # (0:off,1:brdc,2:sbas,3:dual-freq,4:est-stec,5:ionex-tec,6:qzs-brdc,7:qzs-lex,8:vtec_sf,9:vtec_ef,10:gtec)
pos1-tropopt       =saas       # (0:off,1:saas,2:sbas,3:est-ztd,4:est-ztdgrad)
pos1-sateph        =brdc       # (0:brdc,1:precise,2:brdc+sbas,3:brdc+ssrapc,4:brdc+ssrcom)
pos1-posopt1       =on         # (0:off,1:on)
pos1-posopt2       =on         # (0:off,1:on)
pos1-posopt3       =on         # (0:off,1:on)
pos1-posopt4       =on         # (0:off,1:on)
pos1-posopt5       =off        # (0:off,1:on)
pos1-exclsats      =           # (prn ...)
pos1-navsys        =63         # (1:gps+2:sbas+4:glo+8:gal+16:qzs+32:comp)
pos2-armode        =fix-and-hold # (0:off,1:continuous,2:instantaneous,3:fix-and-hold)
pos2-gloarmode     =2        # (0:off,1:on,2:autocal)
pos2-arthres       =3
pos2-arlockcnt     =0
pos2-arelmask      =20         # (deg)
pos2-arminfix      =10
pos2-elmaskhold    =0          # (deg)
pos2-aroutcnt      =5
pos2-maxage        =30         # (s)
pos2-slipthres     =0.05       # (m)
pos2-rejionno      =30         # (m)
pos2-rejgdop       =30
pos2-niter         =1
pos2-baselen       =0          # (m)
pos2-basesig       =0          # (m)
out-solformat      =llh        # (0:llh,1:xyz,2:enu,3:nmea)
out-outhead        =on        # (0:off,1:on)
out-outopt         =on        # (0:off,1:on)
out-timesys        =gpst       # (0:gpst,1:utc,2:jst)
out-timeform       =hms        # (0:tow,1:hms)
out-timendec       =3
out-degform        =deg        # (0:deg,1:dms)
out-fieldsep       =
out-height         =ellipsoidal   # (0:ellipsoidal,1:geodetic)
out-geoid          =internal   # (0:internal,1:egm96,2:egm08_2.5,3:egm08_1,4:gsi2000)
out-solstatic      =all        # (0:all,1:single)
out-nmeaintv1      =1          # (s)
out-nmeaintv2      =1          # (s)
out-outstat        =off        # (0:off,1:state,2:residual)
stats-eratio1      =300
stats-eratio2      =300
stats-errphase     =0.003      # (m)
stats-errphaseel   =0.003      # (m)
stats-errphasebl   =0          # (m/10km)
stats-errdoppler   =1          # (Hz)
stats-stdbias      =30         # (m)
stats-stdiono      =0.03       # (m)
stats-stdtrop      =0.3        # (m)
stats-prnaccelh    =10         # (m/s^2)
stats-prnaccelv    =10         # (m/s^2)
stats-prnbias      =0.0001     # (m)
stats-prniono      =0.001      # (m)
stats-prntrop      =0.0001     # (m)
stats-clkstab      =5e-12      # (s/s)
ant1-postype       =llh        # (0:llh,1:xyz,2:single,3:posfile,4:rinexhead,5:rtcm)
ant1-pos1          =90         # (deg|m)
ant1-pos2          =0          # (deg|m)
ant1-pos3          =-6335367.6285 # (m|m)
ant1-anttype       =*
ant1-antdele       =0          # (m)
ant1-antdeln       =0          # (m)
ant1-antdelu       =0          # (m)
ant2-postype       =llh       # (0:llh,1:xyz,2:single,3:posfile,4:rinexhead,5:rtcm)
ant2-pos1          =36.26146875          # (deg|m)
ant2-pos2          =136.9042704          # (deg|m)
ant2-pos3          =573.24          # (m|m)
ant2-anttype       =*
ant2-antdele       =0          # (m)
ant2-antdeln       =0          # (m)
ant2-antdelu       =0          # (m)
misc-timeinterp    =off        # (0:off,1:on)
misc-sbasatsel     =0          # (0:all)
misc-rnxopt1       =
misc-rnxopt2       =
file-satantfile    =
file-rcvantfile    =GSI_PCV.TXT
file-staposfile    =
file-geoidfile     =
file-ionofile      =
file-dcbfile       =
file-eopfile       =
file-blqfile       =
file-tempdir       =
file-geexefile     =
file-solstatfile   =
file-tracefile     =

inpstr1-type       =ntripcli   # (0:off,1:serial,2:file,3:tcpsvr,4:tcpcli,7:ntripcli,8:ftp,9:http)
inpstr2-type       =off        # (0:off,1:serial,2:file,3:tcpsvr,4:tcpcli,7:ntripcli,8:ftp,9:http)
inpstr3-type       =off        # (0:off,1:serial,2:file,3:tcpsvr,4:tcpcli,7:ntripcli,8:ftp,9:http)
inpstr1-path       =kaiyodai:tuomsat00@mgex.igs-ip.net:2101/CUT07:
inpstr2-path       =
inpstr3-path       =
inpstr1-format     =rtcm3      # (0:rtcm2,1:rtcm3,2:oem4,3:oem3,4:ubx,5:ss2,6:hemis,7:skytraq,8:gw10,9:javad,15:sp3)
inpstr2-format     =rtcm3      # (0:rtcm2,1:rtcm3,2:oem4,3:oem3,4:ubx,5:ss2,6:hemis,7:skytraq,8:gw10,9:javad,15:sp3)
inpstr3-format     =rtcm3      # (0:rtcm2,1:rtcm3,2:oem4,3:oem3,4:ubx,5:ss2,6:hemis,7:skytraq,8:gw10,9:javad,15:sp3)
inpstr2-nmeareq    =off        # (0:off,1:latlon,2:single)
inpstr2-nmealat    =26.37293571 # (deg)
inpstr2-nmealon    =127.143649075 # (deg)
outstr1-type       =off        # (0:off,1:serial,2:file,3:tcpsvr,4:tcpcli,6:ntripsvr)
outstr2-type       =off        # (0:off,1:serial,2:file,3:tcpsvr,4:tcpcli,6:ntripsvr)
outstr1-path       =
outstr2-path       =
outstr1-format     =llh        # (0:llh,1:xyz,2:enu,3:nmea)
outstr2-format     =nmea       # (0:llh,1:xyz,2:enu,3:nmea)
logstr1-type       =off        # (0:off,1:serial,2:file,3:tcpsvr,4:tcpcli,6:ntripsvr)
logstr2-type       =off        # (0:off,1:serial,2:file,3:tcpsvr,4:tcpcli,6:ntripsvr)
logstr3-type       =off        # (0:off,1:serial,2:file,3:tcpsvr,4:tcpcli,6:ntripsvr)
logstr1-path       =
logstr2-path       =
logstr3-path       =
misc-svrcycle      =10         # (ms)
misc-timeout       =30000      # (ms)
misc-reconnect     =10000      # (ms)
misc-nmeacycle     =5000       # (ms)
misc-buffsize      =32768      # (bytes)
misc-navmsgsel     =all        # (0:all,1:rover,2:base,3:corr)
misc-proxyaddr     =
misc-fswapmargin   =30         # (s)
```

細かいオプションの意味は，マニュアルを参照すればなんとなくわかると思いますが，よく確かめた方が良いオプションは，今回の場合は

* pos1-posmode
* pos1-frequency
* pos1-soltype
* pos1-elmask
* pos1-snrmask_L{1,2,5}
* out-solformat
* out-timesys
* out-height
* ant2-postype
* ant2-pos{1,2,3}
* file-rcvantfile

あたりでしょうか．マニュアルと内山大先生の解説によると移動局と基準局が同じアンテナ・機材の場合だとまた異なる設定になるようです．

設定ファイルが準備できたら，あとはPPKの実行です．

```
rnx2rtkp -k conf.txt -o out.pos rover.o base.22o base.nav.22g base.nav.22l base.nav.22n
```

PPKの結果は，GUIアプリがあればすぐに見れるのですが，今回はCLIのためすぐに見れません．取り急ぎ，クイックルックのためのPythonプログラムで確認してみます．

```quick_look.py
import pandas
import matplotlib.pyplot as plt

df = pandas.read_csv(
    "out.pos",
    comment="%", 
    delim_whitespace=True, 
    names=("date", "time", "lat", "lon", "height", "Q", "ns", "sdn", "sde", "sdu", "sdne", "sdeu", "sdun", "age", "ratio")
)
plt.scatter(df.lon[df.Q==1], df.lat[df.Q==1], c="green", s=5)
plt.scatter(df.lon[df.Q==2], df.lat[df.Q==2], c="orange", s=5)
plt.scatter(df.lon[df.Q==5], df.lat[df.Q==5], c="red", s=5)
plt.show()
```

おおよそFix解（Q=1）になっていればオッケーかと思います．

# RINEXファイルへの変換

いきなりRINEXファイルが保存される機材であればよいですが，例えば保存できるのは.ubxファイルだったりすることもあります．この場合は，はじめにRINEXに変換する必要があり，そのためには，convbinというプログラムもコンパイルが必要です．以下はそのTIPSです．

```
cd ~/Desktop/RTKLIB-2.4.3-b34/app/consapp/convbin/gcc
```

makefileを次のように編集

* `CC = gcc-11`を追加
* `BINDIR = `にインストールディレクトリを指定
* `LDLIBS = -lm -lrt`から`-lrt`を削除

```
make
make install
```