目次
はじめに
障害物検知用センサの計測精度を評価したりする場合、その基準とする位置座標を得るためにGPS測量を行うことがあります。RTK-GPSであれば数センチ程度の精度でその位置座標が取得できるので、ミリ波レーダやスキャンレーザの計測テストを行う際は、まず先に基準とするポイントを決めて、その位置座標をGPSで測るという作業を全ポイントでやっておきます。(これが中々大変な作業。。)
この時に問題になるのが、取得された緯度、経度情報をX-Yの平面直交座標に座標変換する処理です。通常なら、その座標変換を行うための原点やスケールなどのパラメータを自分で設定しなければならないのですが、Pythonで使えるパッケージの一つであるpyprojを使えば、この座標変換処理を数行のコードを書くだけで簡単に実行してくれて非常に便利です。
今回の記事では、pyprojの使い方と、それによるGPSのNMEAフォーマットの一つである$GPGGAデータの解析を行ってみた際の結果を紹介します。
参考資料
pyprojの使い方については以下の記事を参考にしました。
また、解析対象であるGPSのNMEAデータのフォーマットについては以下の記事に細かくまとめられていました。
こんな書籍もあるようで、面白そうですね。
サンプルデータ
こんなNMEAデータのテキストログをネットで見つけたので、今回はこれをサンプルデータとして使用しました。２秒間隔で更新されており、データフォーマットの種類もいろいろ記録されているので、NMEAデータの解析を覚えるのには丁度いいと思います。
http://www.oiccam.com/reno/gps/tahoe/nmea/2005-09-25-north-tahoe_nmea.txt
$GPGGAフォーマット
今回扱うのはNMEAデータの一つである$GPGGAというフォーマットで、以下のようになっています。分析したいデータによって抽出すべきフォーマットを選ぶのですが、緯度、経度と一緒に海抜高さや使用衛星数、水平精度低下率を分析したいときはこのフォーマットを抽出します。
NMEAテキストログからの$GPGGAデータ抽出
テキストログの中身は以下のようになっています。気持ち悪いですね(笑)
画像はメモ帳で開いた様子なのでだらだら書き込まれていますが、実際は$から始まって次の$の手前までが一つのデータであり、１行についき１種類のデータで改行されながら書き込まれていきます。
これをまずPythonで開き読み込みます。そして$GPGGAから始まる行だけを取り出していき、それをカンマで区切ることで$GPGGAが持つ各データを得ることが出来ます。最終的には、カンマで区切った各データに対してカラムを割り当てて、時刻をインデックスとしたデータフレームにします。
この処理を行うコードは以下のようになっています。
with open(sample_file_name, 'r') as f: str_nmea_log = f.readlines() file_name_split = sample_file_name.split('-') # read NMEA log and extract $GPGGA for log_line in tqdm(str_nmea_log): if log_line.find('GPGGA') != -1: data_counter += 1 log_line_split = log_line.split(',') if data_counter == 1: date_time_str = file_name_split[0]+'-'+file_name_split[1]+'-'+file_name_split[2]+\ '-'+log_line_split[1] init_date_time = datetime.datetime.strptime(date_time_str, '%Y-%m-%d-%H%M%S') date_time_array = np.append(date_time_array, np.array([init_date_time]), axis=0) gpgga_array = np.append(gpgga_array, np.array([log_line_split]), axis=0) prev_date_time = init_date_time else: date_time = prev_date_time + datetime.timedelta(seconds=2) date_time_array = np.append(date_time_array, np.array([date_time]), axis=0) gpgga_array = np.append(gpgga_array, np.array([log_line_split]), axis=0) prev_date_time = date_time # create DataFrame of GPGGA gpgga_format = ['Data Type','UTC Time','Latitude','North/South','Longitude', 'East/West','Quality','Satellites Num','HDOP','Altitude','Alit Unit', 'Geoid Height','Geo Unit','Comm Time','Check Sum'] gpgga_data_frame = pd.DataFrame(data=gpgga_array, index=date_time_array, columns=gpgga_format)
緯度経度からX-Y平面への座標変換
データの抽出が完了したら、次はいよいよ緯度経度からX-Y平面座標への変換をします。
NMEAの$GPGGAフォーマットに記録されている緯度経度は、～度～分の単位で表されているのですが、X-Y座標への変換を行うには、まずこれらを～度の単位に変換する前処理が必要になります。上記の$GPGGAフォーマット表にある通り、例えば緯度が3541.1439と記録されていればそれは35度41.1439分であり、dddmm.mmmmとなっています。なので、これを～度に変換するにはdddとmm.mmmmに分割して、mm.mmmmを60で割ってあげれば～度になり、あとは分割した際のdddと足せば完了です。
以下のようなコードで処理を実行します。
def transform_lonlat_to_xy(gpgga_data_frame): # convert lat_deg_min to lat_deg lat_dm_array_str = gpgga_data_frame['Latitude'].values lat_dm_array_float = [float(lat_str) for lat_str in lat_dm_array_str] # str -> float lat_d_array_int = [int(lat_dm_float/100) for lat_dm_float in lat_dm_array_float] lat_m_array_float = [(lat_dm_float/100 - lat_d_array_int[idx_lat]) for idx_lat, lat_dm_float in enumerate(lat_dm_array_float)] lat_d_array_float = [lat_m_float*100/60 for lat_m_float in lat_m_array_float] lat_d_array = [(lat_d_array_int[idx_lat]+lat_d_float) for idx_lat, lat_d_float in enumerate(lat_d_array_float)] # convert lon_deg_min to lon_deg lon_dm_array_str = gpgga_data_frame['Longitude'].values lon_dm_array_float = [float(lon_str) for lon_str in lon_dm_array_str] # str -> float lon_d_array_int = [int(lon_dm_float/100) for lon_dm_float in lon_dm_array_float] lon_m_array_float = [(lon_dm_float/100 - lon_d_array_int[idx_lon]) for idx_lon, lon_dm_float in enumerate(lon_dm_array_float)] lon_d_array_float = [lon_m_float*100/60 for lon_m_float in lon_m_array_float] lon_d_array = [(lon_d_array_int[idx_lon]+lon_d_float) for idx_lon, lon_d_float in enumerate(lon_d_array_float)]
緯度経度を～度の単位に変換できたら、今度はpyprojを使ってX-Yの２次元平面直角座標に変換します。
pyprojをまずはインポートします。インストールされていなければ、pipコマンドでインストールしておいてください。
import pyproj
インポートできたら、次は座標系、測地系の設定をします。座標系、測地系の詳しい説明は下記の記事を参照してください。
座標系とは｜座標系の種類・平面直角座標系について解説 - 空間情報クラブ｜株式会社インフォマティクス
測地系とは｜日本測地系・世界測地系 - 空間情報クラブ｜株式会社インフォマティクス
ここでは、世界で最も広く使われているGRS80楕円体を測地系に設定し、それをX-Y平面直角座標系に変換します。
測地系の設定から座標変換、最後に計算したX, Y座標をデータフレームに追加するところまでの処理を以下のコードで実行します。
# transform from grs80 to x-y grs80 = pyproj.Proj(init='EPSG:6668') rect6 = pyproj.Proj(init='EPSG:6680') x_tmp, y_tmp = pyproj.transform(grs80, rect6, lon_d_array, lat_d_array) # North:+, South:- north_south = gpgga_data_frame['North/South'].values y = [-y_tmp[idx_n_e] if n_or_e == 'S' else y_tmp[idx_n_e] for idx_n_e, n_or_e in enumerate(north_south)] # East:+, West:- east_west = gpgga_data_frame['East/West'].values x = [-x_tmp[idx_e_w] if e_or_w == 'W' else x_tmp[idx_e_w] for idx_e_w, e_or_w in enumerate(east_west)] # add x and y to data frame gpgga_data_frame['x'] = x gpgga_data_frame['y'] = y
あまり難しくない実装で行えるので便利ですね。ちなみにこれをpyprojの助けを借りずに自力で座標変換しようとするなら、下記の記事にあるような計算を実装しなければならないのでとても面倒です。
X-Y平面座標の可視化
データの前処理が済んだところで、最後にデータを少し分析してみます。
計算されたX-Y座標の軌跡を描き、それを位置特定品質で色分けしたのが下記の散布図です。
マゼンタ色が単独測位、黄色がディファレンシャルGPSを意味しており、単独測位だと10～15m程度の精度、ディファレンシャルなら2m以内の精度と言われています。場所によって両者が切り替わっていますが、これらは衛星の配置や数によってどんどん変わっていきます。ディファレンシャルGPSについては以下の記事に詳しい説明があります。
https://www.rex-rental.jp/knowledge/gps/gps_003.html
今度は測位に使用された衛星数で色分けしてみましょう。
色が黄色に近づくほど使用衛星数が多いことを表しますが、位置特定品質の色が変わるポイントを比較すると、衛星数が多くなる場所ほどディファレンシャルGPS測位になり精度が向上していそうに見えますね。ただし、軌跡の左半分を見るとあまり関係なさそうにも見えます。今度は、測位状態の良し悪しを判断する指標の一つであるHDOP(水平精度低下率)を見てみます。
この数値が大きくなるほど測位精度が悪くなっていると判断できますが、こう軌跡全体で見るとこれといって精度が悪化しているようには見えないですね。もっと細かい範囲で解析すればもしかしたら変化が起きているかもしれないですね。 ここまで可視化してきた特定品質、使用衛星数、HDOPに加えて海抜高さやジオイド高さなど、測位の精度に関わってきそうなデータ同士で相関を調べてみたら、何か面白い特徴が見つかるかもしれないので、今度見てみようと思います。
