User:Nahainec/Triangulation/Program

はじめに

マップエディタ上に補助線をひいて交点を求めるためのプログラムです。

使い方

1. 観測地点と目標物の方位をファイルに記述
2. プログラム実行
3. できたGPXファイルをアップロード

入力ファイル

P <メモ> <緯度> <経度>    原点座標 緯度・経度の単位は度(分秒は使わない)
M <方位> [<長さ>]         磁方位 単位は度 長さの単位はm 省略時は1000
T <方位> [<長さ>]         真方位 単位は度 長さの単位はm 省略時は1000

M,Tは直前の最後のPを基準にするので 一つの観測点から複数の目標を観測した場合は 一つのPで観測点を示したあと複数のMまたはTで目標を示すことができます。

プログラム実行

上で書いたファイルを標準入力から入れます。 結果のGPXファイルは標準出力から出力されます

アップロード

できたGPXファイルはトラックログと同様にアップロードできます。

制限

磁方位

磁方位を真方位へ変換する偏角の計算が日本近傍を前提としたものなので 他の地域では誤差が大きくなります。

距離

原点からある方位にある距離だけ離れた地点の計算がとても大雑把です。

調べたのですが式が理解できなかったので、理解できる範囲で書いています

肉眼で見える範囲なら平面と仮定できるかなと。
どなたか

• これでどの程度の誤差があるのか計算
• より正確な実装

してもらえませんか?

入力

入力ファイルのチェックはしていないので注意してください。

プログラム

#! /usr/bin/awk -f

# 座標,方位から作図用の GPX ファイルを生成する

# 入力ファイル
#
# P <メモ> <緯度> <経度>    原点座標 緯度・経度の単位は度(分秒は使わない)
# M <方位> [<長さ>]         磁方位 単位は度 長さの単位はm 省略時は1000
# T <方位> [<長さ>]         真方位 単位は度 長さの単位はm 省略時は1000

BEGIN {
    pt = 0;

    printf( "<?xml version=\"1.0\" ?>\n" );
    printf( "<gpx>\n" );
}

# 原点座標
( $1 == "P" ) {
    baselat = $3;
    baselon = $4;
    basename = $2;

    wpt( baselat, baselon, basename );

    # 偏角 (日本近傍限定)
    dlat = baselat - 37.0;
    dlon = baselon - 138.0;
    baseD2000 = 7.619033 + 0.360367 * dlat - 0.127867 * dlon + 0.007033 * dlat * dlat - 0.005333 * dlat * dlon - 0.011250 * dlon * dlon;
}

# 磁方位
( $1 == "M" ) {
    azconv( $2 - baseD2000, $3, $2 );
}

# 真方位
( $1 == "T" ) {
    azconv( $2, $3, $2 );
}

END {
    printf( "  <trk>\n" );
    printf( "    <name>Bounding Box</name>\n" );
    printf( "    <trkseg>\n" );

    for( i = 0; i < pt; i++ )
    {
        trkpt( alat[i], alon[i] );
    }

    printf( "    </trkseg>\n" );
    printf( "  </trk>\n" );
    printf( "</gpx>\n" );
}

function trkpt( lat, lon )
{
    printf( "      <trkpt lat=\"%f\" lon=\"%f\">\n", lat, lon );
    printf( "        <time>%s</time>\n", strftime( "%Y-%m-%dT%H:%M:%SZ", systime() - 9 * 60 * 60 ) );
    printf( "      </trkpt>\n" );
}

function wpt( lat, lon, name )
{
    printf( "  <wpt lat=\"%f\" lon=\"%f\">\n", lat, lon );
    printf( "    <name>%s</name>\n", name );
    printf( "  </wpt>\n" );
}

function azconv( az, len, daz )
{
    len = len / 1000.0;
    if( len < 0.001 )
    {
        len = 1.0;
    }

    s = sprintf( "%s %d", basename, daz );

    # 以下は大雑把な計算  lenが大きくなると誤差が増大する
    vlat = len * cos( az * 3.141593 / 180 ) / 111.6;
    vlon = len * sin( az * 3.141593 / 180 ) / cos( baselat * 3.141593 / 180 ) / 111.6;

    wpt( baselat + vlat, baselon + vlon, s );

    add_array( baselat, baselon );
    add_array( baselat + vlat, baselon + vlon );
    add_array( baselat, baselon );
}

function add_array( lat, lon )
{
    alat[pt] = lat;
    alon[pt] = lon;
    pt++;
}