EurekaMoments

新米エンジニアが一人前を目指す修行の日々を記していくブログです。

GPS NMEA Analysis with Python package pyproj

Introduction

障害物検知用センサの計測精度を評価したりする場合、その基準とする位置座標を得るためにGPS測量を行うことがあります。RTK-GPSであれば数センチ程度の精度でその位置座標が取得できるので、ミリ波レーダやスキャンレーザの計測テストを行う際は、まず先に基準とするポイントを決めて、その位置座標をGPSで測るという作業を全ポイントでやっておきます。(これが中々大変な作業。。)
この時に問題になるのが、取得された緯度、経度情報をX-Yの平面直交座標に座標変換する処理です。通常なら、その座標変換を行うための原点やスケールなどのパラメータを自分で設定しなければならないのですが、Pythonで使えるパッケージの一つであるpyprojを使えば、この座標変換処理を数行のコードを書くだけで簡単に実行してくれて非常に便利です。
今回の記事では、pyprojの使い方と、それによるGPSのNMEAフォーマットの一つである$GPGGAデータの解析を行ってみた際の結果を紹介します。

Reference

pyprojの使い方については以下の記事を参考にしました。

pyproj [いかたこのたこつぼ]

また、解析対象であるGPSのNMEAデータのフォーマットについては以下の記事に細かくまとめられていました。

GPSのNMEAフォーマット

こんな書籍もあるようで、面白そうですね。

Raspberry Pi Gps Using Python 2.7 or 3.4: For Raspbian Jessie Linux Using Gpsd Gps3

Raspberry Pi Gps Using Python 2.7 or 3.4: For Raspbian Jessie Linux Using Gpsd Gps3

Sample Data

こんなNMEAデータのテキストログをネットで見つけたので、今回はこれをサンプルデータとして使用しました。2秒間隔で更新されており、データフォーマットの種類もいろいろ記録されているので、NMEAデータの解析を覚えるのには丁度いいと思います。
http://www.oiccam.com/reno/gps/tahoe/nmea/2005-09-25-north-tahoe_nmea.txt

About $GPGGA format

今回扱うのはNMEAデータの一つである$GPGGAというフォーマットで、以下のようになっています。分析したいデータによって抽出すべきフォーマットを選ぶのですが、緯度、経度と一緒に海抜高さや使用衛星数、水平精度低下率を分析したいときはこのフォーマットを抽出します。
f:id:sy4310:20180901152859p:plain

Extraction $GPGGA data from NMEA text log

テキストログの中身は以下のようになっています。気持ち悪いですね(笑) f:id:sy4310:20180901154529p:plain 画像はメモ帳で開いた様子なのでだらだら書き込まれていますが、実際は$から始まって次の$の手前までが一つのデータであり、1行についき1種類のデータで改行されながら書き込まれていきます。
f:id:sy4310:20180901165517p:plain
これをまずPythonで開き読み込みます。そして$GPGGAから始まる行だけを取り出していき、それをカンマで区切ることで$GPGGAが持つ各データを得ることが出来ます。最終的には、カンマで区切った各データに対してカラムを割り当てて、時刻をインデックスとしたデータフレームにします。
f:id:sy4310:20180901165713p:plain
この処理を行うコードは以下のようになっています。

    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)

Transformation Longitude and Latitude to X and Y position

データの抽出が完了したら、次はいよいよ緯度経度から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の2次元平面直角座標に変換します。
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の助けを借りずに自力で座標変換しようとするなら、下記の記事にあるような計算を実装しなければならないのでとても面倒です。

平面直角座標への換算 計算式

Visualize X-Y Position and Analysis

データの前処理が済んだところで、最後にデータを少し分析してみます。
計算されたX-Y座標の軌跡を描き、それを位置特定品質で色分けしたのが下記の散布図です。
f:id:sy4310:20180901181256p:plain
マゼンタ色が単独測位、黄色がディファレンシャルGPSを意味しており、単独測位だと10~15m程度の精度、ディファレンシャルなら2m以内の精度と言われています。場所によって両者が切り替わっていますが、これらは衛星の配置や数によってどんどん変わっていきます。ディファレンシャルGPSについては以下の記事に詳しい説明があります。
https://www.rex-rental.jp/knowledge/gps/gps_003.html

今度は測位に使用された衛星数で色分けしてみましょう。
f:id:sy4310:20180901182416p:plain
色が黄色に近づくほど使用衛星数が多いことを表しますが、位置特定品質の色が変わるポイントを比較すると、衛星数が多くなる場所ほどディファレンシャルGPS測位になり精度が向上していそうに見えますね。ただし、軌跡の左半分を見るとあまり関係なさそうにも見えます。今度は、測位状態の良し悪しを判断する指標の一つであるHDOP(水平精度低下率)を見てみます。

f:id:sy4310:20180901183317p:plain
この数値が大きくなるほど測位精度が悪くなっていると判断できますが、こう軌跡全体で見るとこれといって精度が悪化しているようには見えないですね。もっと細かい範囲で解析すればもしかしたら変化が起きているかもしれないですね。 ここまで可視化してきた特定品質、使用衛星数、HDOPに加えて海抜高さやジオイド高さなど、測位の精度に関わってきそうなデータ同士で相関を調べてみたら、何か面白い特徴が見つかるかもしれないので、今度見てみようと思います。