GPSレシーバーからのデータを加工

今回は前回に引き続き、USB接続のGPSレシーバーから取得したデータをPythonにてGeoJSON形式に変換し、QGISに表示するまでを説明します。

前回までの内容

前回は、GPSレシーバーの構造や仕組みや出力電文のNMEA0183のフォーマットを確認し、実際にTeraTermで電文を受信しました。

[blogcard url=”https://takamoto.biz/other/gps1/”]

今回の実現イメージ

今回は、USB接続のGPSレシーバーから取得したデータをPythonにてGeoJSON形式に変換し、QGISに表示します。
今回の実現イメージを以下に示します。
GPSレシーバーから出力されるNMEAデータ（テキスト）は、TeraTermのログファイルとして取得します。

各種仕様

使用するNMEAデータの仕様

使用するGPSからのデータ（以降、NMEAデータ）は、以下の緑網掛けの項目です。

NMEA0183フォーマットを解説しているサイトは多くありますが、GISupplyさんのサイトがシンプルでわかりやすいです。また、古野電気さんのGPS製品のサイト（PDF）も参考になります。

GeoJSONの仕様

NMEAデータそのままではQGISに表示できないため、今回はGeoJSON形式に変換します。
GeoJSONは、地物をJSON形式で表現するフォーマットです。
以下、Wikipediaからの引用です。

GeoJSONはJavaScript Object Notation (JSON) を用いて空間データをエンコードし非空間属性を関連付けるファイルフォーマットである。 属性にはポイント（住所や座標）、ライン（各種道路や境界線）、 ポリゴン（国や地域）などが含まれる。 他のGISファイル形式との違いとして、Open Geospatial Consortiumではなく世界各地の開発者達が開発し管理している点で異なる。

引用元：https://ja.wikipedia.org/wiki/GeoJSON

以下は、2つの属性（Prop1、Prop2）を持つ点データの定義例です。

（2つの属性を持つ点のGeoJSONデータ定義例）
{“type”:”FeatureCollection”,
　　”features”: [
　　　　{“type”:”Feature”,
　　　　　　”geometry”:{“type”:”Point”,”coordinates”:[139.1,35.3]},
　　　　　　”properties”:{“Prop1″:”111″,”Prop2″:”aaa”}
　　　　}
　　]
}

トップの”FeatureCollection”が地物の集まりを表し、その下に”Features”配列が存在します。
”Features”配列のメンバーとして”Feature”が存在し、”Feature”には形状を表す”Geometry”と属性を表す”properties”が存在します。図で表すと以下のようになります。

Pythonプログラム

プログラムは、メインとGeoJSONを書き込むクラスの2つに分けました。
メインプログラムは以下となります。

import csv
import datetime
from GeoJSONWriter import GeoJSONWriter

# TeraTermログ（NMEA）ファイル
NMEA_FILE = ".\\teraterm.log"

# ファイルオープン
with open( NMEA_FILE )  as fileNMEA:
    # CSVオブジェクトの生成
    csvNMEA = csv.reader(fileNMEA)
    # 最初のRMC受信フラグとディクショナリの宣言
    enableRmc=False
    gjWriter = GeoJSONWriter('.\\gps.geojson')
    dicGpsdata = []
    for rowNMEA in csvNMEA:
        # 行Loop
        if rowNMEA[0].startswith('$GPRMC') and len(rowNMEA)==12:
            # GPS/Transit Data
            # 3.状態で有効データを判定
            if( rowNMEA[2]=='V' ):
                continue
            enableRmc=True
            # 2.時刻、10.日付（UTC→JSTに変換）
            receiveDatetime = datetime.datetime.strptime(rowNMEA[9]+rowNMEA[1],"%d%m%y%H%M%S.%f") + datetime.timedelta(hours=9)
            # 4.緯度を十進数に変換
            latitude = float(rowNMEA[3][0:2]) + float(rowNMEA[3][2:]) / 60
            # 6.経度を十進数に変換
            longitude = float(rowNMEA[5][0:3]) + float(rowNMEA[5][3:]) / 60
            # 8.移動速度をノットからキロメートルに変換
            speed = -1.0 if rowNMEA[7]=='' else float( rowNMEA[7] ) * 1.852
            # 9.方位
            direction = -1.0 if rowNMEA[8]=='' else float( rowNMEA[8] )
        elif rowNMEA[0].startswith('$GPGGA') and len(rowNMEA)==15:
            # Fix Data
            # 8.衛星使用数
            satelites = 0 if rowNMEA[7]=='' else int( rowNMEA[7] )
        elif rowNMEA[0].startswith('$GPGSA') and len(rowNMEA)==18:
            # active satellites
            # 3.タイプ（測位次元）
            dimension = 0 if rowNMEA[2]=='' else int( rowNMEA[2] )
            # 16.位置精度低下率(PDOP)
            pdop = 100.0 if rowNMEA[15]=='' else float( rowNMEA[15] )
            if( enableRmc ):
                dicGpsdata.append(str(receiveDatetime)+','+str(round(latitude,5))+','+str(round(longitude,5))+','+
                                str(round(speed,2))+','+str(direction)+','+str(satelites)+','+str(dimension)+','+str(round(pdop,2)))
                gjWriter.setGeometry(round(longitude,5),round(latitude,5))
                gjWriter.setPropertie('receiveDatetime' ,receiveDatetime)
                gjWriter.setPropertie('speed' ,round(speed,2))
                gjWriter.setPropertie('direction' ,direction)
                gjWriter.setPropertie('satelites' ,satelites)
                gjWriter.setPropertie('dimension' ,dimension)
                gjWriter.setPropertie('pdop' ,round(pdop,2))
                gjWriter.Write()
                # RMC待ち状態とする
                enableRmc=False
        elif rowNMEA[0].startswith('$GPGSV') and len(rowNMEA)==20:
            # satellites in view (Skip)
            pass
        else:
            # 不整レコード
            print("不正レコードです：[" + str(csvNMEA.line_num) + "]" +str(rowNMEA))

以下に、メインプログラムを簡単に説明します。

[table id=11 /]

次に、GGeoJSON書き込み用のクラスは、以下となります。

class GeoJSONWriter:
    # コンストラクタ
    def __init__( self ,inOutFile ):
        # 出力ファイルのオープン（とりあえず失敗・例外を考慮しない）
        self._jsonfd = open(inOutFile ,'w' ,encoding='utf-8')
        # ジオメトリとプロパティをクリア
        self.__Geometry = None
        self.__Properties = None

    # デストラクタ
    def __del__(self):
        # フッターを書き込み
        self._jsonfd.write( "\n" + "]}")
        # ファイルのクローズ
        self._jsonfd.close()

    # ジオメトリの設定
    def setGeometry(self ,inLon ,inLat):
        # FIGTYPE=4 : 文字
        self.__Geometry = "\t" + "{" + "\"" + "type"+ "\"" + ":" +"\"" + "Feature" +"\"" + "," + "\n"
        self.__Geometry+= "\t" + "\"" + "geometry"+ "\"" + ":{"
        self.__Geometry+= "\"" + "type"+ "\"" + ":" + "\"" + "Point" + "\"" + ","
        self.__Geometry+= "\"" + "coordinates"+ "\"" + ":"
        self.__Geometry+= "[" + str(inLon)+","+str(inLat) + "]"

    # プロパティの設定
    def setPropertie(self ,inName ,inValue):
        # 2つ目以降か？
        if self.__Properties==None:
            # 1つ目
            self.__Properties = "\t" + "\"" + "properties"+ "\"" + ":{"
        else:
            # 2つ目以降
            self.__Properties+= ","
        if type(inValue) is list:
            # 文字列リスト
            self.__Properties+= "\"" + inName + "\""+":" + "\"" + "".join(inValue) + "\""
        else:
            # 文字列項目
            self.__Properties+= "\"" + inName + "\""+":" + "\"" + str(inValue) + "\""

    # ファイルの書き込み
    def Write(self):
        if self._jsonfd.tell()==0 :
            # 最初はヘッダーを書き込み
            self._jsonfd.write("{" + "\"" + "type" + "\"" + ":" + "\"" + "FeatureCollection" + "\"" + "," + "\"" + "features" + "\""  + ": [" + "\n")
        else:
            # 最初の書き込み以外
            self._jsonfd.write(",\n")
        self._jsonfd.write(self.__Geometry + "}," + "\n")
        self._jsonfd.write(self.__Properties + "}"+"}")
        # ジオメトリとプロパティをクリア
        self.__Geometry = None
        self.__Properties = None

以下に、GeoJSON出力クラスの内容を簡単に説明します。

[table id=12 /]

結果として出力されるGeoJSONのイメージは、以下となります。

{“type”:”FeatureCollection”,”features”: [
{“type”:”Feature”,
“geometry”:{“type”:”Point”,”coordinates”:[139.747,35.63974]},
“properties”:{“receiveDatetime”:”2020-04-23 12:42:56″,”speed”:”-1.0″,”direction”:”-1.0″,”satelites”:”4″,”dimension”:”3″,”pdop”:”6.4″}},
{“type”:”Feature”,
“geometry”:{“type”:”Point”,”coordinates”:[139.74684,35.63914]},
“properties”:{“receiveDatetime”:”2020-04-23 12:42:57″,”speed”:”0.0″,”direction”:”-1.0″,”satelites”:”4″,”dimension”:”3″,”pdop”:”6.4″}},

QGISでの表示

GeoJSONデータも作成できましたので、QGISを使って地図に表示してみましょう。
QGISのデータソースマネージャからGeoJSONファイルを指定します。

今回取得したデータとOSM地図を重ね合わせると、以下のような表示となりました。
私が実際に歩いたルートは赤破線ですので、誤差はさほど大きくないことがわかります。

誤差の原因を分析

誤差の原因を分析してみましょう。
今回出力したGeoJSONには属性として”測位日時”、”測位次元数”、”使用衛星数”、”精度低下率”（PDOP）を保持しています。
この情報を使って誤差の原因を分析してみます。

以下の地図は、10秒ごとに以下の情報をラベル表示しています。
　分：秒（測位次元数，使用衛星数，精度低下率）

①のエリアで誤差が大きくなっています。
この時のデータは以下のようになっており、捨てるべきデータだったことがわかります。

• 測位次元数が2次元
• 使用衛星数が0
• PDOPが50.0（数字が大きいほど誤差が大きい）

②のエリアでも誤差が大きくなって水没しています。
以下の条件が1つでも満たされると誤差が大きくなっているように見えます。（この時は、水際ぎりぎりを歩いていたので、水没とはいえ大きい誤差ではないと思っています。）

• 測位次元が2次元
• 使用衛星数が6未満
• PDOPが3以上

③のエリアでは、ほぼ正確な位置が取得できています。
使用衛星数も多く、PDOPも低い値となっています。

④のエリアでは、使用衛星数も多くPDOPも低いのですが、誤差が大きくなっています。
この原因は、ビルによるGPS衛星の電波の反射だと考えられます。
10階建て程度であまり高くはないのですが、ビルの谷間を歩いていました。その結果、GPS的には正しく測位していますが、ビルに反射した電波を使用しているため誤差が大きくなったと考えられます。

まとめ

ここまで、以下を確認しました。

• GeoJSON形式の仕様確認
• PythonプログラムによるNMEAデータのGeoJSON形式への変換
• QGISを使った測位結果の表示と誤差分析

今回使用したGPSモジュールは、旧式のせいか測位開始まで40秒程度の時間が必要です。
測位にはコツがあって、開始までの40秒間は動かずに待ち続け測位が開始（GPSモジュールのLEDランプが点滅状態となります）されてから動くようにすると、精度が高くなります。
最初に測位したときはそのことを知らず、いきなり歩いてしまったので、全般的に誤差が大きくなり、再実施を余儀なくされました。（下図）

ただ、なぜここまで誤差が大きくなったのかは不明です。