s-taka.org

高橋 賢

Lang:

シェルスクリプトを用いたRTK基準局の座標決定

category: gnss
tags: rtk

日本国内でのRTK基準局の座標決定には、国土地理院の電子基準点を利用することが多いと思います。しかし、いつも利用する基準点は決まっているのに、毎回、基準点を検索してダウンロードするのは億劫です。しかし、国土地理院では、ftp (file transfer protocol)でのデータダウンロードも可能にしてくださっています。

また、自ら観測データを集めて、RTKLIBのRTKPOSTで対話的に処理すると、RTKへの理解が深まります。しかし、アンテナ変更のたびに計算を行い、エクセルなどでfixした座標を抽出し、さらに平均を求める作業は億劫です。しかし、RTKLIBにはrnx2rtkpという非対話的に実行できるアプリケーションも同梱されています。

そこで、国土地理院電子基準点の観測データのダウンロードからRTK基準局の座標決定までの一連の処理を行うシェルスクリプトを作成してみました。

  1. 国土地理院の電子基準点ftpアカウント（ID、パスワード）を取得していること。
  2. 電子基準点の番号、緯度、経度、楕円体高があらかじめ検索済みであること。
  3. curl、convbin、rnx2rtkp、awkがパスの通っているディレクトリにあること。

を仮定しています。衛星配置は12時間周期、または、24時間周期ですので、RTK基準局の座標決定には24時間以上の観測が必要です。一方、基準局のアーカイブにはRINEX形式が用いられることが多いですが、RINEX形式では複数ファイルをそのままでは結合できません。複数観測ファイルを結合するために、ここでは生データであるu-blox形式ファイルを用います。生データはストリーム形式ですので、そのファイルをそのまま結合できます。

下記の日付YEAR、 MONTH、 DAY、電子基準点局番号POS、座標LAT、LON、HEIGHT、国土地理院アカウントGSIID、GSIPW、および、RTK基準局観測データのURIRTKURIは皆様の環境に合わせて適宜、書き換えて下さい。

#!/bin/bash

YEAR=2021　     # date
MONTH=03
DAY=26
POS=1157        # GSI reference station number
LAT=34.39704732 # position
LON=132.33981921
HEIGHT=217.55
GSIID= 国土地理院アカウント # 書き換えてください
GSIPW= 国土地理院パスワード # 書き換えてください

T1=`date +%s -d "$YEAR$MONTH$DAY"`
T2=`date +%s -d "${YEAR}0101"`
DAYNO=`echo $T1 $T2 | awk '{printf "%03d", ($1-$2)/60/60/24+1}'`
GSIURI=ftp://terras.gsi.go.jp/data/GRJE_3.02
echo Getting GSI observation data on $YEAR/$MONTH/$DAY at $POS
curl -Ss -u $GSIID:$GSIPW -O $GSIURI/$YEAR/$DAYNO/$POS${DAYNO}0.${YEAR:2:2}o.gz

echo -n "Getting RTK observation data: "   # GNSS生データを得る一例です
RTKURI=https://phys.info.hiroshima-cu.ac.jp/gnss/f9p/
for hourcode in {a..x}; do
    echo -n "$YEAR$MONTH$DAY$hourcode "
    curl -Ss -O $RTKURI/$YEAR$MONTH/$YEAR$MONTH$DAY$hourcode.ubx
done
echo
cat $YEAR$MONTH${DAY}?.ubx > $YEAR$MONTH$DAY.ubx && rm $YEAR$MONTH${DAY}?.ubx

echo Converting RTK raw to RINEX
convbin $YEAR$MONTH$DAY.ubx -v 3.02 -ti 1 && rm $YEAR$MONTH$DAY.ubx

echo Executing RTK
rnx2rtkp -p 2 -h -t -u \
    $YEAR$MONTH$DAY.obs $YEAR$MONTH$DAY.nav $POS${DAYNO}0.${YEAR:2:2}o.gz \
    -l $LAT $LON $HEIGHT \
    > $YEAR$MONTH$DAY.pos

echo Averaged position:
awk '$6==1{n+=1; lat+=$3; lon+=$4; height+=$5} END {printf "%.7f %.7f %.3f\n", lat/n, lon/n, height/n}' $YEAR$MONTH$DAY.pos

# EOF

RTK基準局の観測データが取れている日では、別の日で処理しても、経度の小数点以下7桁目がほとんど変わりませんでしたので、センチメートルオーダーで座標特定ができていると思います。私のところでは、生データを5 Hzで記録していますが、1 Hzでも充分かも知れません。