目次：INDEX

ライブラリ・モジュール

ツール

スクリプトの実行方法など

Pythonを用いてGNSS関係の計算を行う上で必要となるライブラリです。 様々な状況に応じて動的にエフェメリスを取り出す・・・なんて処理をC/C++で書きたくなかったのでライブラリ化しました。 GPSだけではなく、他の測位衛星のデータも処理できるようにモジュール構造とクラス構造を工夫したつもりです。 個人的にはこれを使ってしばらく試行錯誤してみます。

基本的使い方

適当なフォルダに解凍し、環境変数PYTHONPATHを通してください。 その後、どこからでもimportできます。 各モジュールの利用方法はモジュールの最下部に記述しているセルフテスト部をご覧ください。

準拠楕円体の定義クラス及び重力加速度クラスを整備した、測地学クラスです。 準拠楕円体にはWGS84とGRS80を定義済みです。 重力加速度クラスにはデフォルトでGRS80用およびWGS84用を用意しています。 正規重力や楕円体高によって変わる重力加速度を計算することができます。 重力加速度を提供するGravityBaseクラスは準拠楕円体定義を受け取って動作するジェネリッククラスです。 従って、準拠楕円体さえ定義すれば、重力加速度の計算が簡単にできるようにしています。

コードはC#で開発しており、dllをインポートすればVisual Basicでの利用も可能です。

基本的使い方

緯度31.2°における正規重力値を取得するコード例を以下に示します。 GravityクラスはデフォルトではGRS80を利用していますので、この例ではGRS80楕円体での計算値になります。 現時点では近似式を用いたメソッドのみ実装していますが、実用上問題はないと思います。

＊コードのHTML変換にはC#toHTMLを利用させて頂きました。

本クラスライブラリは、座標を表すBLH（緯度/経度/楕円体高）・XYZ座標系を処理する構造体，緯度・経度方向の長さを定義する構造体，領域を表現する構造体，NMEAフォーマットログデータを一括処理できるクラスを含んだGNSS向けのライブラリです。 NMEAのログを使用して何らかの計算を実行するプログラムに活用可能です。 距離の計算やBLHとXYZの相互変換に対応しています。 また、座標構造体は足し算引き算に対応しています。

コードはC#で開発しており、dllをインポートすればVisual Basicでの利用も可能です。

以下に本ソースコードが含む、クラス・構造体を示します。 割と力を入れたのはXYZ座標系とBLH座標系及びNMEAパーサークラスです。 座標系に関しては相互変換や距離を求めるメソッドなどを実装しており、手前味噌ながら上手くできていると思います。 NMEAパーサーは取得済みのログファイルを処理するには余り不自由しない程度には実装しています。 一行ずつ処理する機能は未実装ですが、現時点でも大した行数もなくNMEAログを扱えるようになるのでメリットは大きいと思います。

[注意]以下のクラスダイアグラムは2012/5/19リリースのクラス群に比べて古いものです。 若干の変更や機能拡張が行われていることに注意してください。

含んでいるクラス・構造体の一覧　及び　座標系を定義する構造体

WGS-84の定数定義クラス　と　長さや領域を定義する構造体・クラス

NMEパーサークラス　と　NMEAクラスが包含している観測データ格納・処理用クラス

[注意] プロジェクトファイルをVisual Studioで開くと、初期状態では依存しているライブラリの参照設定のエラーが出ているかと思います。 その場合は、binフォルダ内にある該当dllファイルを参照するように設定しなおしてください。

リリース日	プログラムパッケージ	更新内容
2012/10/8	GNSS_20121008.zip	以下の更新を行いました。 長方形領域定義RectangleFieldクラスに、中心座標と半径で定義するコンストラクタを追加 上記の追加したコンストラクタは、指定した半径の円が収まるような長方形を定義するコンストラクタです。
2012/8/24	GNSS_20120824.zip	以下の更新を行いました。 長方形領域定義RectangleFieldクラスに、領域比較用メソッドIsContain(RectangleField)を追加 NMEAログを処理するNmeaReaderクラスに実装していたOpenWithDialog()を、ファイルを開けたかどうかの成否を返すように仕様を変更 NmeaReaderクラスに、任意の文字列がNMEAであるかをチェックするIsNmea(string)を実装 測位情報を扱う、PositioningInfoおよびPositioningInfosクラスの名前空間をGNSS直下へ移す PositioningInfoクラスへGetData(int index)を実装 Length構造体へToString()を実装 XYZ構造体をEcefへ改名 SDFという、マイナーなフォーマットのGPSログを処理するクラスを新設 RINEX・NMEA・SDF　名前空間を一つの名前空間TextDataへ移した GPS・ガリレオ等の衛星の名前空間はSatellite名前空間へ引っ越し GEONET名前空間はObservationNetwork空間へ引っ越し
2012/7/23	GNSS_20120723.zip	以下の更新を行いました。 GRS80への対応 測地系の定義を測地学プロジェクトライブラリへ移築 長方形領域定義RectangleFieldクラスに、領域を拡張するためのExtend(double pad)メソッドを追加
2012/7/6	CSLibTest_20120706.zip	GNSSクラスライブラリはC#で作成していますが、これをVisual C++から利用したサンプルプロジェクトファイルです。 参照は図に示すようにソリューションエクスプローラーより追加してください。 プログラムを実行すると下の図のようになります。 VC++からdllを利用する方法についてはこちらのページも参考にして下さい。
2012/7/3	GnssLib_20120703.zip	最新のライブラリとテストプロジェクトを同梱しました。 テストプロジェクトはテスト用静的クラスを利用して2点間の距離を計算・表示するコンソールアプリで、コマンドプロンプトから実行可能です。 実行すると下図の様に表示されます。 スペースの関係上縮小させています。 画像のみを別タブで開くと見やすいと思います。 以下更新内容です。 XYZ構造体へ文字列で初期化するコンストラクタを追加 XYZ構造体に実装していた地心距離rプロパティをNormへ改名 BLH構造体へ文字列で初期化するコンストラクタを追加 WGS84パラメータを表す静的クラスを、測地系間の変換に備えるために普通のクラスにして継承ができるように変更 テキストファイル処理時の利便性向上のために、BLHとXYZ構造体のコンストラクタ生成時に文字列を受け取れるように変更しました。 また、測地系間の座標変換が簡単になるように、WGS84クラスからstatic修飾子を取り去って普通のクラスとしました。 WGS84クラスはをWgs84に改名しています。 このWgs84クラスが返すパラメータはG873をデフォルトにしています。 普通のクラスにしたことによって継承ができるようになったため、デフォルトをG1150やG730へ変更することも容易になりました。 インターフェイスの実装も行ったため、将来的に他の測地系を実装するのも簡単になっています。 ただし、普通のクラスになってしまったためパラメータを利用する際はインスタンスの生成が必要となり、コード中では1行余計に食うようになっています。 増えた手間よりも利便性の向上が勝っていますのでご容赦ください。
2012/6/9	GNSS_20120609.zip GNSS_help_20120609.chm	最新のライブラリとその最新ヘルプファイルです。 以下更新内容です。 NMEAクラスのOpen()メソッドへエラー対策を追加 NMEAクラスにファイル名を保持させるように変更 RectangleFieldクラスにコピーコンストラクタを追加 BLH構造体に距離を求めるメソッドを3手法（実質は2つ）を追加 XYZ構造体の+と-の演算子の定義が間違っていたのを訂正 GnssTestクラスに距離を求めるアルゴリズムの比較用メソッドTwoBlhToDistanceTest()を追加 今回のアップデートはXYZ構造体の演算子定義ミスという致命的バグの修正を含んでいます。 たったの3要素しかないベクトルの足し算と引き算を間違えるとは…。 これにより距離の計算精度が出るようになりました。 また、短距離での距離計算ではベクトルの引き算でも精度十分でかつ高速動作可能ということが分かりました。 メモ：XYZ構造体に実装している地心距離rのプロパティですが、これは将来Normへ変更する予定です。
2012/6/1	GNSS_help_20120601.chm	ライブラリの最新ヘルプファイルです。
2012/6/1	GNSS_20120601.zip	GGAセンテンスの処理プログラムにバグがあったのを修正しました。 また、テスト用の静的クラスGnssTestクラスを追加しました。 なお、プロジェクト言語はC#ですがVisual Basic(VB)でもdllを参照して機能を呼び出すことができることを確認しました。 コードの変換は不要の様です。 [いつになるかもわからない更新予定] 1) 測地系評価メソッドの追加 2) C++で作成した衛星軌道演算・測位演算のラッパ―を実装 3) 測地系変換処理
2012/5/19	GNSSclass20120519.zip	作りためたGPS関係のクラスをGNSS名前空間上に移してプロジェクトとしてまとめました。 せっかくなので拡張に備えてGNSS名前空間内にQZSSなどの名前空間も用意しています（今んとこは空ですが）。 利用するには、Visual Studioのソリューションに「既存のプロジェクトの追加」で追加して参照設定するか、 もしくはbinフォルダ内にあるGNSS.dllとGNSS.XMLファイルを参照設定に追加して下さい。 クラス内のメソッドや依存関係が以前リリースしたコードと若干異なっていることに注意が必要ですが、名前空間の整理もしたつもりなので便利になったと思います。 今後は、測地系のパラメータ比較用のメソッドを実装する予定です。
2011/12/1	GPSandGISclass.chm	GPSクラスの説明書です。 Sandcastle Help File Builderを使用して、ヘルプファイルを出力してみました。 クラスのメンバを確認することが可能です。 html版はこちらです。 なお、3か月以内に公開予定のGISクラスについても同梱しています。 何かコメントがあればよろしくお願い致します。
2011/10/20	GPScs_v2.zip	NMEAパーサを使ってu-bloxのバイナリを含んだログデータを処理しようとしたのですが、できませんでした。 そこで行頭がNMEAセンテンスである行のみを抽出したファイルを処理しようとしたのですが、これもデータ中に含まれていた不正文字によってエラーを起こしました。 今回の更新ではその辺のバグ対策（応急処置）を施しています。 なお、ライセンス文書も一緒にしておきました。
2011/9/22	GPS.cs	初公開です。

基本的使い方

本クラスを用いたサンプルコードを以下に示します。 この例ではNMEAのログファイルを読み出して緯度・経度・時刻を読み出し、それを基にDEMの標高データ系列を取得＆保存しています。

本クラスを用いた処理の一例

本ライブラリはGISに関するライブラリです。 数値標高モデル（DEM）を始めとするラスタ地図に関するクラスや、Google Earth用のKMLに関するクラスを内蔵しています。 DEMを画像にするなんて世に数多あるソフトウェアでできますが、自作のプログラムに組み込める便利なライブラリは見当たりませんでした。 という訳で、作ってみた次第です。

ラスタ地図クラスは任意の型の構造体やクラスや標準型（int, long, float, double）を格納できるマップクラスです。 任意の型用および比較が可能な型用さらに標準型用の3種類を用意しています。 例えば、深さ方向を持たせた構造体を使えば、3次元データを表現できます。 また、C#では構造体に関数を実装可能ですので日変化する地磁気等の表現が可能です。 データの取得方法は、緯度経度もしくはマップ内の整数アドレスで指定できます。

ラスタ地図クラスを利用するDEMクラス群も定義しています。 本クラス群には、DEMを管理するマネジメントクラスや読み込み用クラスが含まれます。 DEMのマネジメントクラスを用いれば、多数のDEMを結合して一枚の大きな地図を作ることが可能です。 以下に示している図はそうやって作った図です。 本クラスには、GPSで計測した緯度経度に沿った標高値を取得するとか、測位環境を地形から評価するとか、河川の支流領域を抽出するとか使い方は色々あると思います。 DEMの読み込みクラスが対応しているフォーマットは、国土地理院が発行しているGMLファイルのみとなっています（2012/6時点）。 ただし、同じフォーマットであればデータポイント数が違っていても、地域が日本でなくても処理可能です。 読み込みクラスは、簡単な構造ですので任意のフォーマットへ対応した新しいクラスを定義することは容易だと思います。 DEMの読み込み可能なファイル数はメモリ容量の制約を受けますが、少なくともファイル情報なら3400個以上，5 mメッシュデータフルで700個程度は読み込めることを確認済みです。 マルチコアマシンで動かせば素早い処理が可能です。 なお、本クラスは任意の構造体を格納したマップを作成する際のサンプルでもあります。 ご参考にされて下さい。

KMLに関するクラスには、点を表すPointクラスと面を表すPolygonクラスを同梱しています。 これらのクラスを利用してGoogle Earth上へデータを書きだすことが可能です。 例えば、GPSのログデータを点で表現するなどの応用が考えられます。

図　DEMを結合して作った画像
（国土地理院の10 mメッシュデータを用いています。）


図　保持しているDEMが一目でわかる
DEMファイル保持状況をKMLファイルへ出力する機能も備えている。


図　傾斜角マップと傾斜角に閾値を設定してマスクしたマップ

基本的使い方

プロジェクトに含まれるサンプルコード（GisTest.cs）より、フォルダを指定して走査されたDEMファイルの一覧をKMLファイルへ出力するサンプルコードを抜粋し、以下に示します。 このサンプルでは、読み込んだファイル数が0でなければファイルへ出力させるという処理を行っています。

無論、DEMファイル読み込み後はBLH構造体を用いて位置座標を指定することでその地点の標高値を取得することができます。

＊コードのHTML変換にはC#toHTMLを利用させて頂きました。

実験データを整理するために整備した時刻関係の演算メソッドを集めたモジュールです。
GPS週番号や通算日の計算に使用可能です。 リクエストがあればGPS時刻やGPS週秒の演算処理も追加したいと思います。

基本的使い方

インポートしてから使ってください。。

図　Pythonによる時刻演算スクリプトのセルフテストの様子

実験データを整理するために整備した時刻関係の演算メソッドを集めたプログラムです。
最新バージョンにおいて、以下のメソッドを整備しています。

基本的使い方

クラスにはしていませんので、普通にrequire 'time_culcv6.rb'と宣言してソースヘッダでインクルードしてから使用して下さい。
関数の返値は基本的には数値となります。

元々はマイコンでの時刻管理のために整備したプログラムですが、これにrubyで作った便利関数を一部移植したものです。

基本的使い方

時刻ライブラリのC++版です。 C言語版に比べ、より使いやすくよりスタイリッシュにより高機能になっています（笑）。 DOW（day of year:年間の通算日）とGPS week（GPS週番号）が追加されています。

リリース最新版の搭載機能

#ifndef    HDR_time_personal
#define    HDR_time_personal                    //このファイルをインクルードしているかはこれが定義されているかでわかる。

namespace GPS
{
    class Time
    {
    // メンバ変数
    private:
        // static member
        static const int               _SUNDAY                  = 0;            // 便宜上の曜日定義。GPS時刻に合わせて、日曜日を0とする。
        static const int               _MONDAY                  = 1;
        static const int               _TUESDAY                 = 2;
        static const int               _WEDNESDAY               = 3;
        static const int               _THURSDAY                = 4;
        static const int               _FRIDAY                  = 5;
        static const int               _SATURDAY                = 6;
        static const int               _WEEK_ERROR              = 7;
        static const int               _DAY_A_WEEK              = 7;            // 1週間の日数
        static const long int          _TIME_COUNT_A_DAY        = 86400l;       // 一日の時間の長さ[s] 3600s * 24h = 86400s
        static const long int          _TIME_COUNT_A_WEEK       = 604800l;      // 1週間の時間[s], 24h * 7day (1week) = 604800 sec
        static const double            GPS_START_EPOCH_JD       = 2444244.5;    // ユリウス日表現によるGPSの開始エポック
        static const long int          GPS_START_EPOCH_MJD      = 44244l;       // 修正ユリウス日表現によるGPSの開始エポック
        static const unsigned long int GPS_START_EPOCH_DATE     = 19800106ul;   // GPS時刻の開始エポック（1980/1/6 sunday）
        // other
        long int            _inc;           // 一時的に記憶する、加算量
        unsigned long int   _clock;         // 時刻の数値化したもの。例：“020409”-> [20409]（2時4分9秒）
        long int            count;          // 時刻の一日の通算値。例：CLOCK == 1235（0時12分35秒）-> [12 * 60 + 35]. 一日の時刻を秒で与える。
        unsigned long int   date;           // 日付の数値化したもの。例：2009年4月2日 -> 20090402
        long int            GPS_time;       // GPS時刻。
        long int            _GpsTow;        // Time of Week:GPS週秒。日曜日の午前0時を基準に、秒カウント。
        int                 GPS_week;       // GPS週番号（うるう秒分だけ、UTCの週変わりよりも早くカウントアップされます）
        int                 DOY;            // day of year:通日
        int                 week;           // 曜日。日曜日を0とし、土曜日を6とする。エラーがある場合は7となる。
        int                 year;           // 年
        int                 month;          // 月
        int                 day;            // 日
        bool                error;          // 演算エラー。1:error, 0:normal
        int                 timeSystem;     // 時刻系（デフォルトでは、GPST or UTC）
        bool                _countChangeFlag; // _clockの更新用に用いる
        bool                _clockUpdateReq;  // _clockの更新用に用いる
    public:
        static const int    GPST                                = 16;           // UTCを基準にしたうるう秒
        static const int    UTC                                 = 0;
        static const int    JST                                 = 32400;
    // アクセッサ
    private:
        void     setClock(unsigned long int a){ _clock = a; return; }
        void     setCount(long int a){ count = a;           return; }
        void     setDate(unsigned long int a){ date = a;    return; }
        void     setGpsTime(long int a){ GPS_time = a;      return; }
        void     setGpsTow(long int a){_GpsTow = a;         return; }
        void     setGpsWeek(int a){GPS_week = a;            return; }
        void     setDOY(int a){     DOY = a;                return; }
        void     setWeek(int a){    week = a;               return; }
        void     setYear(int a){    year = a;               return; }
        void     setMonth(int a){   month = a;              return; }
        void     setDay(int a){     day = a;                return; }
    public:
        long int getClock()
        {
            if(this->_countChangeFlag) (void)this->Update(); // 時刻全般を更新（フラグも処理される）
            if(this->_clockUpdateReq)                        // countの更新がなされていないといけないので、この処理は後
            {
                this->_clock = Count_to_Clock(this->count);  // hhmmssを更新
                this->_clockUpdateReq = false;
            }
            return this->_clock;
        }
        long int getCount(){            if(this->_countChangeFlag) (void)this->Update(); return this->count; }
        unsigned long int getDate(){    if(this->_countChangeFlag) (void)this->Update(); return this->date; }
        long int getGpsTime(){          if(this->_countChangeFlag) (void)this->Update(); return this->GPS_time; }
        long int getGpsTow(){           if(this->_countChangeFlag) (void)this->Update(); return this->_GpsTow; }
        int      getGpsWeek(){          if(this->_countChangeFlag) (void)this->Update(); return this->GPS_week; }
        int      getDOY(){              if(this->_countChangeFlag) (void)this->Update(); return this->DOY; }
        int      getWeek(){             if(this->_countChangeFlag) (void)this->Update(); return this->week; }
        int      getYear(){             if(this->_countChangeFlag) (void)this->Update(); return this->year; }
        int      getMonth(){            if(this->_countChangeFlag) (void)this->Update(); return this->month; }
        int      getDay(){              if(this->_countChangeFlag) (void)this->Update(); return this->day; }
        void     setError(bool a){  error = a;             return; }
        bool     getError(){                               return error; }
    // メソッド（メンバ関数とも言う）
    public:
        Time();                                                         // コンストラクタ, 時刻系はUTCとなる
        Time(unsigned long int date_p);                                 // コンストラクタ　オーバーロード1, 時刻系はUTCとなる
        Time(unsigned long int date_p, unsigned long int GpsTow);       // コンストラクタ　オーバーロード2, 時刻系はGPSTとなる
        Time(unsigned long int date_p, long int clock_p);               // コンストラクタ　オーバーロード3, 時刻系はUTCとなる
        Time(const Time& obj);                                          // コピーコンストラクタ
        ~Time();                                                        // デスコンストラクタ
        void    Increase_Clock(long int);                               // 時刻を加減算する
        static unsigned long int Count_to_Clock(long int);              // 通算時間から時刻へ変換した数値を生成
        static unsigned long int Clock_to_Count(unsigned long int);     // 時刻から通算時間へ変換した数値を生成
        void    Set_yyyy_mm_dd(unsigned long int date_p);               // 時刻の初期化メソッド
        void    Set_yyyy_mm_dd(unsigned long int date_p, unsigned long int TOW_p);  // オーバーロードを利用している
        void    Set_yyyy_mm_dd(unsigned long int date_p, long int clock_p);         // オーバーロードを利用している
        void    Day_after(void);                                        // 一日日付を進める
        void    Day_before(void);                                       // 一日日付を戻す
        const char*     GetWeekName(void);                              // 曜日名を返す
        void            SetTimeSystem(int);                             // 時刻系をセットする
        const char*     GetTimeSystemName(void);                        // 時刻系名を返す
    private:
        char            Week_culc(Time *);                              // 曜日を計算する
        long int        yyyymmdd_to_ModifiedJulius(int y, int m, int d);// 修正ユリウス日を求める
        int             date2ddd(Time *);                               // 1年の内の通算日を求める
        int             Date2GPSweek(Time *);                           // GPS週番号を求める
        static long int Compare(const Time& l,const Time& r);           // 比較メソッド
        void            UpdateGpsTime(long int inc);                    // GPS時刻の更新用関数
        void            Update(void);                                   // 時刻の更新用関数
    // 演算子
    public:
        Time&    operator=(const Time& obj);
        long int operator-(const Time& r) {return Compare(*this, r);}   // -演算子オーバーロード
        bool     operator<=(const Time& r) const{return Compare(*this, r) <= 0;}
        bool     operator>=(const Time& r) const{return Compare(*this, r) >= 0;}
        bool     operator<(const Time& r) const{return Compare(*this, r)  <  0;}
        bool     operator>(const Time& r) const{return Compare(*this, r)  >  0;}
        bool     operator==(const Time& r) const{return Compare(*this, r) == 0;} 
        bool     operator!=(const Time& r) const{return Compare(*this, r) != 0;}
    };
}

#endif

基本的使い方

NMEAフォーマットで記録したGPS受信機のログかを解析するパーサーのRuby版です。 ログから緯度・経度を抜きだしたり、GSVセンテンスを処理することが可能です。

使い方は非常に簡単で、例えば以下の様にして呼び出します。 ここで、lineにはファイルから読みだしたログが1行分入っているものとします。
nmea = GetNMEAstruct(line)
nmeaには解析結果が入ったオブジェクトが格納されます。

天文学に関する.NETライブラリです。 これからちょっとずつ拡張していくつもりです。 日の出時刻の荒っぽい計算値が欲しい方はどうぞ。

放送暦を基に衛星軌道を計算するPythonスクリプトです。 GNSS汎用ライブラリ for Pythonのサンプルを兼ねています。

基本的使い方

コマンドの実行時に、スクリプトファイル名を指定する様にして下さい。 なお、コマンド冒頭の"python"は省略しないと動作しないかもしれません。

k bitサンプリングのIFデータが隙間なく格納されているバイナリデータをk bitずつ切り出して0詰めした8 bitデータを保存するPythonスクリプトです。 ソフトウェアGPS用のデータ変換にお使いください。

基本的使い方

コマンドの実行時に、スクリプト名に続いてバイナリファイルのパスを指定して下さい。 なお、コマンド冒頭の"python"は省略しないと動作しないかもしれません。

コマンドを実行した様子

ハンディGPS受信機で測位した結果はカシミール3Dを使えば簡単に地図上にプロットできます。 カシミールに取り込んだGPSのログは各種フォーマットでエクスポートできますが、本ソフトはその中のsdfフォーマットファイルを処理して数値標高マップから標高を拾い出して結合させてテキストファイルに書き出します。 用途としては、GPSの高度方向の測位結果評価などでしょうか？ 個人的には気圧高度推定と組み合わせた経路推定アルゴリズムの検証用に使用しています。

基本的使い方

まずは国土地理院からダウンロードできる基盤地図情報（標高・GML形式）を予めダウンロードして解凍して一つのフォルダにまとめておいて下さい。 起動後は下記の図に従って起動後にsdfファイルをドラッグアンドドロップするだけでOKです。

図1 プログラムの起動


図2 DEMフォルダの選択


図3 表示メッセージ


図4 起動したらドラッグアンドドロップしよう


図5 変換が完了すると出るメッセージ


図6 ファイルが作出されている

国土地理院からダウンロードできる基盤地図情報（標高）を、Excelで確認するのに使用したスクリプトです。 DEMデータはXML形式になっていますが、これを簡易な方法で中身を確認する方法としてExcelで処理できるようにcsvファイルへ変換します。 地図データを取り込めるソフトには、 カシミール3D や ArcGIS や MANDARA がありますので、よほどのもの好きでないと必要ないかもしれません。 自作プログラム作成時のデータの確認等にご利用下さい。

=begin
--------------------------------------------------------------
国土地理院の標高地形データをExcel用にCSVにするプログラム

Excelで等高線を描くと行列を入れ替えても上下が入れ替わったままだけど、何とか表示できる
左上が座標0,0とされていることを確認した。

    開発者：K.Morishita Kumamoto-University
    開発履歴：2011/6/20  テストに成功
--------------------------------------------------------------
tab size: 4
--------------------------------------------------------------
=end

fname = ARGV[0]
txt = File.readlines(fname)
line = 0
line += 1 while txt[line].index("<gml:high>") == nil      # データサイズが書かれた行を探す
str = txt[line]                                           # データサイズの記入された文字列を取得
p str
str.sub!("<gml:high>","")                                 # 必要のない文字列を削除
str.sub!("</gml:high>","")
size = str.split(" ")
x = size[0].to_i + 1                                      # 格納されているのは最後のアドレスなので1を足すとサイズになる
y = size[1].to_i + 1
puts "sizeof(x,y) == #{x},#{y}"                           # 取得したサイズを出力

line += 1 while txt[line].index("<gml:tupleList>") == nil # 標高データの開始タグ位置を取得
line += 1                                                 # こっからデータが始まる
File.open("DEM2CSV_converted.csv", 'w') do |f|
    _x = 0
    _y = 0
    while txt.length > line do
        f.write  ",#{txt[line].chomp.split(",")[1].to_f}" # 標高データを取得&書き込み
        _x += 1
        if _x == x
            f.write "\n"
            _y += 1
            _x = 0
            break if _y == y
        end
        line += 1
    end
end

基本的使い方

国土地理院からダウンロードしたファイル名を指定して実行するだけです。 国土地理院の基盤地図情報（デジタル地図）ダウンロードはこちらからどうぞ。

本スクリプトは、平成23年東北地方太平洋沖地震の後に停止したGEONETの電子基準点を抽出するために作成したスクリプトです。 GEONETのサイトからFTPダウンロードして得られるファイルの数は膨大です。 例えば、2011/3/11の観測データをダウンロードすると2,474個のファイルが含まれていました。 観測データはこの半分ですが、これらを手作業でチェックするのは現実的ではありません。

基本的使い方

先ず、ブラウザかFFFTPの様なソフトを用いて、国土地理院から電子基準点の観測データをダウンロードします。 全観測点のダウンロードは、データ量が半端ないので回線が細いと40分以上かかります。 ダウンロード後は適当なフォルダに解凍して下さい。 その後、サンプルと同じファイル構造へと変更して下図に有るようにスクリプトを実行すればOKです。

GEONETの基準点コード一覧表

観測データと処理スクリプト

観測データのチェック実行

抽出された観測データの基準点住所を確認する

GEONETの基準点コード一覧表と突き合わせた結果

スクリプトを実行すると得られるファイルの中身
先ず、右下のファイルを生成して、次に真ん中のファイルを作ります。 Excelで開くなら、時刻の表示設定を“時刻”へ設定して下さい。

国土地理院が発行している電子基準点のRINEXファイルは観測日から3日経たなければ1日分の結合したファイルがリリースされません。 短時間の処理ならばらばらにリリースされるRINEXファイルを用いても良いのですが、そうも言っていられないケースもあります。 しかし手動で結合というのはバカらしいですしプログラムに任せましょう、というコンセプトで本スクリプトを作りました。 3日も待てない方はどうぞお使いください。

ところで、“RINEX 結合”でググるとTEQCというフリーソフトがかかるのですが、私には使用法が十分には分かりませんでした。 それなら、プログラムが何をやっているのか管理できる自作ソフトの方がましかと思って作っています。 国土地理院のWEBサイトではTEQCを薦めてあるので、将来的には使い方をマスターしたい所ですが現時点では別に自作ソフトで構わないし。。。

基本的使い方

まずは結合させたいRINEXファイルを一つのフォルダに入れ、そこにここでリリースしているスクリプトをコピーして下さい。 その後、そのフォルダを作業フォルダにしたコマンドプロンプトを起動させてスクリプトファイル名を入力することで実行します。 プログラムが走り始めると、カレントディレクトリ中のファイルを全走査してRINEXファイルをとにかく結合します。 なお、実行にはpyscripterを利用しても構いません。

スクリプトのあるフォルダをカレントディレクトリにしてコマンドプロンプトを立ち上げ（cdで移動してもOK）

スクリプトの実行画面

実行後にはcombined.txtというファイルが作出されています

時刻と時間幅を記載したファイルに基づき、指定された時間のRINEX観測データを抜き出します。 RINEXヘッダを保持した状態で（もちろん、時刻情報は更新します）保存します。 現時点での用途としては、静止中の観測データのみを切り出せるかなぁといったとことです。 ただし、この場合は時刻のリストを他の手段で作る必要があります。 私の場合だと、RINEXと同時にNMEAも保存しておいて、NMEAのVTGセンテンスから速度情報を得て時刻のリストを自動的に作っています。

基本的使い方

まずは分割させたいRINEXファイルと時刻のリストを一つのフォルダに入れ、コマンドプロンプトからスクリプトを実行してください。 引数の指定方法はreadme.txt内に記述しています。 実行後、元の観測ファイル名に「picked_rinex_03d」（03dは3桁０詰めであることを表しています）を付加したファイルが生成されます。

コマンドプロンプトでの実行の様子

コマンドプ実行の結果出力されるRINEXファイル
ハイライトされているファイルが出力されたRINEXファイルで、csvファイルは時刻のリストファイルです。

時刻のリストをExcelで開いたところ

時刻のリストはテキストエディタでも編集可能です

u-blox受信機が出力したログファイルをRTKLIBで処理するとRINEX形式の放送暦が得られます。 このプログラムは、その放送暦のパラメータを整理して読みやすいCSV形式に変換するRubyスクリプトです。 実際に放送された放送暦と、IGSの提供する放送暦の比較を行いたい場合に有用かもしれません。
「IGSの放送暦を整理してCSVファイルに出力する」を元に開発しています。

処理後のCSVファイルは、統計処理ソフトRでの処理がやり易いように、PRN番号を使用した層別ラベルを含んでいます。

基本的使い方

コマンドの実行時に、ファイル名を指定する様にして下さい。

コマンドを実行した様子

多数の観測データがある場合、GUIを操作しながら測位計算理するのは面倒です。 とある実験で移動と停止を繰り返しつつ観測データを収集したのですが、これを処理するのがメンドウだったのでスクリプトを書いてみた次第です。 測位演算はRTKLIBのrnx2rtkp.exeを繰り返し呼び出すことで実行しています。 一つのフォルダに観測データと航法データを集めておけば一度に多数のRINEXファイルを処理できます。 ググっても同系統のプログラムは簡単には見つかりませんでしたので有用だと思います。

基本的使い方

コマンドの実行時に、ファイル名を指定する様にして下さい。

解凍したフォルダは、RTKLIBの実行ファイルの中に配置する

上記の部分を実行の度に確認し、必要なら変更する

コマンドを実行した様子（実行中）

実行すると、posファイルが作成されます

操作の様子

RTKLIBが出力した測位演算結果ファイル（*.posファイル）を処理して、測位座標に関して統計処理を行うPythonスクリプトです。 多数の観測データがある場合、GUIを操作しながら統計処理するのは面倒です。 とある実験で移動と停止を繰り返しつつ観測データを収集したのですが、これを処理するのがメンドウだったのでスクリプトを書いてみた次第です。

用途別に、「RtklibPosProcessing.py」と「RtklibPosProcessingWithTimeList.py」という、2つのスクリプトを同梱しています。 「RtklibPosProcessing.py」はスクリプトのあるフォルダ内全ての測位結果（*.pos）が対象となり、 「RtklibPosProcessingWithTimeList.py」は一つの観測結果を、時刻を指定する形で統計処理します。 個人的には、「RtklibPosProcessing.py」をstatic解の解析に、 「RtklibPosProcessingWithTimeList.py」をKinematic解の解析に使用しています。

基本的使い方

コマンドプロンプトを立ち上げて、スクリプトファイル名を入力後にエンターキーで実行されます。 詳しくは配布ファイル内のreadmeや、下にある動画をご覧ください。

ファイル構成（これら以外に、ライブラリファイルが必要です）

コマンドの実行状況

スクリプト実行の結果、2つのファイルが出力される

RtklibPosProcessingWithTimeListの実行の様子（こちらは引数が必要です）

RtklibPosProcessingWithTimeListの実行の結果出力されたlog.csv

RtklibPosProcessingWithTimeListの実行の結果出力されたlog lite.csv

操作の様子

u-blox受信機が出力したバイナリデータを含むログファイルからNMEAセンテンスを抽出するスクリプトです。 正規表現を使ったパターン合致によって抽出していますので途中にバイナリデータが入ったりしません。 NMEAそのものを抽出しても何に使えるか分かりませんが、整備したパーサには他の機能もあるので色々試してみてください。

基本的使い方

コマンドの実行時に、ファイル名を指定する様にして下さい。

コマンドを実行した様子

GPS衛星の軌道計算をする際にいちいちブラウザを使ってファイルをダウンロードしなくても済むようにスクリプトを作りました。 期間の指定により大量の放送暦を一度にダウンロードすることが出来ます。 熊本大学の回線では、1日分の放送暦当たり5秒程度でダウンロードできました。
なお、コードの開発そのものはWin7 x64で行っていますが、スクリプトですので他のマシンでも動作します。

IGSの放送暦をダウンロードした様子

IGSがリリースしている放送暦のパラメータを整理して読みやすいCSV形式に変換するRubyスクリプトです。 コメントをかなり入れてますので、更に処理するパラメータを増やすことも容易だと思います。 ダウンロードファイルには、IGSの放送暦・Rubyによる処理プログラム・処理後CSVファイル・グラフ化しているxlsxファイルを同梱しました。 xlsxファイルを見て貰うと分かりますが、昇交点赤経のパラメータが意外な動きをしているのが分かります。 って、合ってますよね?

[2011/1/10 追記]ファイル容量制限に引っ掛かって定期的に削除されていたので同梱する放送暦・画像ファイルを減らしました。

基本的使い方

コマンドを実行した様子

左：元のRINEXフォーマットファイル　右：処理後のCSVファイル

以下の様に経度，緯度，標高の順に並んだデータを読み込んで、これをGoogle Earthで読み込めるようにkmlファイル化します。 ダウンロードファイルには、あるGPS衛星一日分のサンプルデータを添付しています。 MATLABを使わなくても軌道を確認出来て便利です。 ECEF座標系のx,y,zのデータであっても、割と簡単なソースコードの改造で同様の処理が行えるはずです。

**データフォーマット例**
173.398474103809,-49.5641186950898,20085089.7646485
-179.113828242437,-52.1760665266914,20117614.7985721
-170.594272392161,-53.6038773264183,20150831.0478571
-161.65496191649,-53.7147177897012,20184165.1599795
-153.071294391038,-52.506460875898,20217045.3987616
-145.481787881562,-50.1098409607123,20248911.3256447

以下の様にx，y，zの順に並んだデータを読み込んで、これを緯度・経度・楕円体高に変換してログファイルに出力します。 このプログラムは、PRN2GPS衛星の3ヵ月分のサンプルデータを同梱しています。 計算で求めた衛星の位置をKMLファイル化するために使えます。

計算には地表付近では十分な精度を持つ近似式を利用していますので地表付近のデータなら大丈夫なのですが、衛星の位置になると誤差が無視できないと思います。 あくまでGoogle Earthで表示して衛星の動きを確認するという大雑把な把握に使う目的でお使用ください。

**データフォーマット例**
-17161519.32,14031432.89,14255041.22
-21751946.48,14187762.58,4204428.294
-22034958.15,12592697.41,-7009818.374
-19338810.66,7637797.372,-16290262.97
-16064923.38,-550937.301,-21115181.62
-14181466.87,-10019779.71,-20238205.18
-14074306.28,-17978999.7,-13967987.89

ECEF座標系の座標値から緯度・経度・楕円体高を計算するGUIソフトです。 計算アルゴリズムは坂井さんの「GPSのための実用プログラミング」p.28の近似式を使っています。 電子基準点をその公開されているECEF座標とハンディGPS受信機を使って探し出すのに使ったことがあります。 書籍には近似精度は実用上十分とありましたが、確かに10cm以内の誤差しかないような気がします。 あ、ちなみに地球中心近くの座標を入力しても正常には計算してくれません。 あしからず。

NMEAフォーマットで記録したGPS受信機のログから、UTC時刻・緯度・経度だけを抽出するサンプルプログラムです。 緯度と経度はデフォルトだとddmm.mmmmのフォーマットになっていますが、これをdd.ddddddと度に単位を変換して出力します。 位置情報の取得にはGGAセンテンスを使用しました。

サンプルプログラムにはRuby初心者でも分かる様に出来る限りコメントを付けました。 また、機能拡張が簡単に出来るように書いていますので、NMEAセンテンスをいじるのにはかなり便利なツールだと思います。

サンプルデータにはへミスフェア（Hemisphere）社製のA100と云う受信機で記録したログとその処理結果を添付しています。
以下は処理例です。

073727.00,130.729382166667,32.813639
073728.00,130.7293825,32.8136386666667
073729.00,130.729383,32.8136386666667
073730.00,130.729383333333,32.8136385
073731.00,130.729383666667,32.8136383333333
073732.00,130.729384333333,32.8136383333333
073733.00,130.729384666667,32.8136381666667

[2011/1/10 追記]ファイル容量制限に引っ掛かって定期的に削除されていたので同梱する画像ファイルの品質を下げました。

NMEAフォーマットで記録したGPS受信機のログを指定数分の一に間引くことでスリムにするプログラムです。 GGAセンテンスをトリガーにして動作するように書いていますが、コードをテキストエディタで変更すれば他のセンテンスでも上手く行きます。 デフォルトでは、トリガーがGGA，付随してZDAがあればそれも出力するようにしています。

GPS受信機から得られたNMEAフォーマットのログデータから未測位データを削除します。 カシミール3DにはNMEAログデータを読み込む機能がありますが（プラグインで対応）、未測位データがあると下図の様に律義に緯度経度0,0まで線を引いてくれます。 このうざい機能を避けるために、このプログラムは未測位データを削除します。 本ソフトはWindowsマシンで実行できるGUIソフトです。

rubyというプログラミング言語を使った、年月日からGPS週番号と曜日を求めるプログラムです。 rubyを既にインストールしている環境であれば直ぐにでも使えます。 下図のスクリーンショットに有る様にコマンドを打って貰うと直ぐに計算結果が出力されます。

アルゴリズム的には、入力された日付から一度ユリウス日（JD）を求めて、GPS週番号の起点であるJD2,444,244.5を引いて、さらにそれを7で割るという単純なものです。 JDはある日を起点とした連続通算日です。 ちなみに、GPSの標準エポックは1980/1/6 Sunday 00:00:00.0です。 GPS週番号はGPS時刻系に基づいて計算する必要があります。 従って、うるう秒を考慮して日曜日の日付変更前後の日時は入れない方が無難でしょう。