トップ 差分 一覧 ソース 検索 ヘルプ PDF RSS ログイン

独り言日記(2009/12)

独り言日記

地球の自転が1日/回転でない理由と地軸の傾きの役割(2009/12/31)

おまけではありますが、地球の自転が1日/1回転でない理由として同一時間で昼と夜が逆転しない、というのがあります。後、地軸が傾いているために日本では四季が生まれるというのがありますね。

このあたりは偶然の運動だと思いますが、よくできています。惑星によっては回転の周期と地軸の傾きは差があるらしく、地球はたまたま好条件がそろっていたので生物が住めたのでは、というのは感慨深いところです。

しかし宇宙関連は知識として面白いですねぇ。たしか中学の理科でこれくらいのことなら習ったはずだったんですが、当時は自分の理解力が薄かったです(^_^;;。星空も再現してみたいけど、どんどん意図からそれてる気がするので そろそろ軌道修正せねば。そういえば、昔PC-98で星空を再現するソフトがあったなぁ。フリーも製品版もあったと記憶してます。結構出来がよかったものでした。3DCGソフトは再現するにはうってつけの環境ですので、興味と腕に覚えのある方はプラネタリウムを実装してもらいたいですね〜。

やはりこのあたりにご興味のある方もいらっしゃるようで、各種ソフトがありますね。

http://www.vector.co.jp/vpack/filearea/win/edu/science/space/

雲の写真(2009/12/31)

元旦は冷え込むそうなので、本日は絶好の雲日和でした。

携帯カメラにて撮影。12/31 約14〜15時の時間帯で撮っています。縦線が入っているのは仕様です。天使降臨みたいな光が差し込む表現も撮れました。

これを3DCGで再現するためのリファレンスとしたいのですが、やっぱり半球的なマッピング(観測点を中心とした半球じゃなくて、地球の中心からの半球)がよい感じになりそうですね。

露出調整はカメラで自動的に行われるため 空撮影では色合いがバラバラになりますねぇ。

続・太陽の入射方向計算(2009/12/31)

太陽の入射角方向の計算の続きです。

入力情報として、

日付を年(int m_year)、月(int m_month)、日(int m_day)、時間を時間(int m_hour)、分(int m_minutes)、秒(int m_second)、場所を緯度(double m_lat)、経度(double m_lng)とします。日時指定はUTCです。

また、defineとして以下を定義しています。

///< 地球が自転を1回行うときの経過日数
#define EARTH_DAY_ROT     0.9973

///< 太陽から見たときの地軸の傾き(度数)
#define EARTH_AXIS_ANGLE  23.4

///< 円周率
#define S_PI             3.1415926535

手順としては以下の流れ。

  1. 春分の日(03/20) 12:00:00の経過秒数を計算
  2. 指定日時の経過秒数を計算
  3. 経度と経過秒数を合わせた水平方向の回転角度を計算
  4. 緯度(θ)と水平方向の回転角度(φ)の極座標からデカルト座標に変換
  5. 地球上の対象位置における3軸の直交ベクトルを計算(Y-upで上向きが空方向)
  6. 地軸の傾きにて直交ベクトルを回転させる
  7. 太陽の公転を考慮した太陽が入射する水平方向の向きを計算
  8. [7]の太陽の入射向きを、[5]の直交ベクトルで表される座標系に変換

あくまでも入射角度を求めるだけなので、地球の大きさや太陽-地球の距離などは不要です。では、プログラムの流れに沿って解説します。

春分の日(03/20) 12:00:00の経過秒数を計算

double dx, dy, dz, sinT, cosT, sinF, cosF, r, cosV, sinV;
double ex, ey, ez;
int baseSec, iVal, allYearSec, nowSec, oneDaySec, iSec;
double eRotSec, dRot;
vec3 axisX, axisY, axisZ, sunDir, v;
mat4 rMat, eMat;
vec4 v4;

// 1日の経過秒数
oneDaySec = 24 * 60 * 60;

// 1年における経過秒数
allYearSec = oneDaySec * 365;

// 地球の自転の基準は「春分の日(3/20)にて、UTC 12:00に緯度0/経度0」
// このときの1月1日 00:00からの経過秒数を計算する
iVal = m_CalcDays(m_year, 3, 20);
baseSec = iVal * oneDaySec + m_CalcSec(12, 0, 0);

緯度0/経度0の位置にて、春分の日(03/20)の12:00:00にちょうど真上に太陽がきますので、ここを地球の自転および太陽回りの公転の開始とします。

このときの01/01 00:00:00から数えたときの経過秒数を計算しています。

指定日時の経過秒数を計算

// 指定日時での経過秒数
nowSec = m_CalcDays(m_year, m_month, m_day) * oneDaySec +
     m_CalcSec(m_hour, m_minutes, m_second);

指定の日時での01/01 00:00:00から数えたときの経過秒数を計算します。

そして、以下の計算が春分の日から数えた経過秒数になります。

// 春分の日の12:00を基準とした経過秒数
iSec = (nowSec - baseSec + allYearSec) % allYearSec;

経度と経過秒数を合わせた水平方向の回転角度を計算

地球の地軸を中心とした日時による回転と、経度指定の回転は同一方向での回転になりますのでまとめます。

先に、地球が一回転するときの経過秒数を求めます。精度のため小数値もほしいのでdouble計算。「EARTH_DAY_ROT = 0.9973」です。

eRotSec = EARTH_DAY_ROT * (double)oneDaySec;

eRotSecは、0.9973 * (24 * 60 * 60) = 86166.72となります。

春分の日からの経過秒数であるiSecをeRotSecで割った余りにて、地球の自転による極座標での位置(回転角度)が求まります。

// 地球の経度方向の回転角度を計算
dRot = fmod((double)iSec, eRotSec) * 360.0 / (double)eRotSec;

dRot = 0の場合は春分の日の12:00:00の位置になります。真上から太陽系を見たときに+Z方向が春分の日の開始位置であるので-90度引きます。これに経度であるm_lngを加えると水平回転の回転角度となります。

dRot = m_lng + (dRot - 90.0);
while(dRot < 0.0) dRot += 360.0;
while(dRot >= 360.0) dRot -= 360.0;

緯度(θ)と水平方向の回転角度(φ)の極座標からデカルト座標に変換

これは極座標→デカルト座標に変換する処理になります。

緯度方向の回転角度は「θ = 90 - m_lat」、水平方向の回転角度は「φ = dRot」になります。

// 緯度経度より、極座標からデカルト座標に変換 (Y-up)
sinT = sin((90.0 - m_lat) * S_PI / 180.0);
cosT = cos((90.0 - m_lat) * S_PI / 180.0);
sinF = sin(dRot * S_PI / 180.0);
cosF = cos(dRot * S_PI / 180.0);

r = 1.0;
dx = r * sinT * cosF;
dz = -(r * sinT * sinF);
dy = r * cosT;

右手系(Y-up、Zは手前方向)の極座標からデカルト座標(普通のXYZ直交軸の座標系)変換を行っている点に注意してください。

半径rの球とした場合は、ここで計算した(dx, dy, dz)が緯度経度と日時による回転を加えた地球上の位置、となります。ただし、この段階ではまだ地軸の傾きは考慮していません。

r=1.0としているので(dx, dy, dz)は単位ベクトルに相当します。これは、地球上のこの位置における上向きのベクトルに相当します。

地球上の対象位置における3軸の直交ベクトルを計算(Y-upで上向きが空方向)

この位置におけるXYZ直交座標を計算します。この座標位置の+X方向が東、-X方向が西、-Z方向が北、+Z方向が南になるようにします。Y方向は垂直上向きになるので(dx, dy, dz)がそのままY軸方向に相当します。

北方向を-Z方向に当てはめたいので、傾く前の地軸方向である(0, 1, 0)と(dx, dy, dz)の外積をとります。この計算結果がX方向の向きになります。残りは、YとX軸の外積で求まりますね。ということで、以下のような計算になります。

// 対象となる位置でのXYZ直交軸を計算
axisZ = vec3(0, 1, 0);
axisY = vec3((float)dx, (float)dy, (float)dz);
axisX = product(axisZ, axisY);
axisZ = product(axisX, axisY);
axisX = normalize(axisX);
axisY = normalize(axisY);
axisZ = normalize(axisZ);

productが外積計算の関数、normalizeがベクトルの正規化関数です。

地軸の傾きにて直交ベクトルを回転させる

// XY平面にて-23.4度の地軸の傾きを加える
rMat = mat4::rotate(vec3(0, 0, 1), (float)(-EARTH_AXIS_ANGLE * S_PI / 180.0));

v4 = vec4(axisX, 0.0f) * rMat;
axisX.x = v4.x;
axisX.y = v4.y;
axisX.z = v4.z;

v4 = vec4(axisY, 0.0f) * rMat;
axisY.x = v4.x;
axisY.y = v4.y;
axisY.z = v4.z;

v4 = vec4(axisZ, 0.0f) * rMat;
axisZ.x = v4.x;
axisZ.y = v4.y;
axisZ.z = v4.z;

ここで計算したrMatは以下のようなZ軸を中心とした-23.4度の回転行列になります。

0.917755 -0.397148 0.000000 0.000000
0.397148  0.917755 0.000000 0.000000
0.000000  0.000000 1.000000 0.000000
0.000000  0.000000 0.000000 1.000000

回転は春分の日の12:00に、太陽から真正面に地球を見据えたときのZ軸中心の回転となります(+Z方向が地球から太陽に向かう方向とする)。それぞれの直交軸をrMatの回転で変換してあげると、地軸の傾きを考慮した地球の指定位置上の座標系となります。

eMat = mat4::identity;
eMat[0][0] = axisX.x; eMat[0][1] = axisX.y; eMat[0][2] = axisX.z;
eMat[1][0] = axisY.x; eMat[1][1] = axisY.y; eMat[1][2] = axisY.z;
eMat[2][0] = axisZ.x; eMat[2][1] = axisZ.y; eMat[2][2] = axisZ.z;

太陽の公転を考慮した太陽が入射する水平方向の向きを計算

これは冬至(12/21)からの経過日数より太陽周りの地球の位置(向き)を計算しています。

// 冬至を基準とした1年の経過日数
iVal = (m_CalcDays(m_year, m_month, m_day) - m_CalcDays(m_year, 12, 21) + 365) % 365;

// 太陽の入射の向きを計算
dRot = ((double)iVal / 365.0) * 360.0;
sinF = sin(dRot * S_PI / 180.0);
cosF = cos(dRot * S_PI / 180.0);

sunDir.x =  (float)cosF;
sunDir.z = -(float)sinF;
sunDir.y =  (float)0.0;

sunDirは太陽から地球への入射方向になります。太陽からの光の入射は、地球に対して公転面に平行に入ることになります(上記の計算のようにXZ平面での向き計算でOK)。これは、太陽と地球の距離がかなり離れていて、かつ太陽が地球に比べて圧倒的に大きいからです。

[7]の太陽の入射向きを、[5]の直交ベクトルで表される座標系に変換

残すは、太陽光の向きsunDirを地球表面上のeMatの座標系に変換します。

v4 = vec4(sunDir, 0.0f) * inv(eMat);
sunDir = vec3(v4.x, v4.y, v4.z);

invは4x4行列の逆行列を求める関数です。eMatは直交ベクトルで構成されますので、あらかじめ転置してinv関数を使わないようにしても問題ありません。

ここで計算されたsunDirが地球の指定位置に入る太陽からの入射ベクトルになります。

と、これで一通りの太陽の入射方向計算は完了となります。ちょっと極座標が入りますが 基本はベクトル行列計算で、じっくり考えていくと特に難しいというものではないというのが分かるかと思います。

これを応用すると、人工衛星が地球を周期的に旋回するときの地球の自転を相殺した場合の軌道の計算や、GPSで行っている計算方法などがだいたい想像できそうですね。もともとはGPSで使用する目的のためにこの部分を掘っていたのですが、太陽光の向き計算に応用できた、という流れで解説してみました。

さて、本日は今年最後の日です。駄文で長文な日記を読まれた方、お付き合いありがとうございます。来年もできる限り使えるかどうかは微妙な(^_^;;)ネタを書いていって、読み物としても面白い内容にしていければと思っております。

来年はじめに(できれば正月期間にて集中して)、Shadeの太陽光計算プラグインをバージョンアップしますね。太陽の向きに加えて、日記で書いていた内容ですが背景色反映/太陽光の色/雲、を入れる予定としています。

ちょっとしたシミュレーションですので、たとえば地球の自転が1日/回転だとどういう現象が起きそうか、地軸が傾いていることでおきる現象は?などを大晦日の晩に記載して締めくくろうと思います。

ということで、残すところわずかですが よいお年を。

太陽の入射方向計算(2009/12/30)

本日は移動日でちょっとだけ時間があるので記載開始です。おそらく今年最後の日記ネタになります。数回に分けて、緯度経度と日時が与えられた場合の太陽の入射方向を計算する流れを書いていきます。

実際のソースは、以前公開したShadeの太陽光計算プラグインのEarthUtil.cpp/hが関連する部分になっています。

関連知識は主にWikipediaで仕入れています。関連ワードとして地球、自転、公転、地軸、などでググれば情報は大量に出てきますので、それからあってるかどうか試行錯誤していただければ幸いです。プログラマ向けの手順解説が少ないと感じましたので私のほうで書いてみようかと。

月の日数

1〜12月の日にちについて。私はまったく覚えていないのですが、たしか語呂合わせの覚え方があったはず。

131
228
331
430
531
630
731
831
930
1031
1130
1231

合計すると365日になります。

うるう年は、西暦が「4の倍数年」でかつ「100で割り切れない or 400で割り切れる」で2月が29日になります。今回はうるう年は無視してます。

入力データ

  • 緯度経度(緯度は-90〜0〜+90、経度は0〜360とする)
  • 日時(UTC(協定世界時)で指定のこと)

緯度経度については、以前書いた日記の内容を参考に。地球上の絶対位置指定になりますね。経度は-180〜0〜+180でも同義です。

日時に関してはUTC指定です。「月/日 時:分:秒」で指定。これは標準時間となり日本の場合は+9時間進んでいることになります。この日時での経度0を基準として計算します。ですので、日本時間が与えられている場合に標準時間を求めるには-9時間引いてあげる必要があります。

以前公開した太陽計算プラグインでは、時間を-9して

hour2 = hour - 9;
hour2 = (hour2 + 24) % 24;

としていたのですがこれは間違いでした。-9引いてマイナスになる場合は日を繰り下げる、日が1より小さくなる場合は月を繰り下げる、としないといけませんね。

日本時間の月、日、時間が指定されているとすると以下のような感じでUTCに変換します。

dHour  = 時間(0-24);
dMonth = 月(1-12);
dDay   = 日(1-31);
standardMedian = 時差(日本の場合は9)

dHour = dHour - standardMedian;     // UTCに変更
if(dHour < 0) {
    dHour += 24;
    dDay--;
    if(dDay < 1) {
        dMonth--;
        if(dMonth < 1) dMonth = 12;
        switch(dMonth) {
            case 1:
                dDay = 31;
                break;

            case 2:
                dDay = 28;
                break;

            case 3:
                dDay = 31;
                break;

            case 4:
                dDay = 30;
                break;

            case 5:
                dDay = 31;
                break;

            case 6:
                dDay = 30;
                break;

            case 7:
                dDay = 31;
                break;

            case 8:
                dDay = 31;
                break;

            case 9:
                dDay = 30;
                break;

            case 10:
                dDay = 31;
                break;

            case 11:
                dDay = 30;
                break;

            case 12:
                dDay = 31;
                break;
        }
    }
}

月日指定を1/1からの累積日数に変換

365日で太陽の周りを一周しますので、月日を1/1からの累積日数に変換します。以下は年指定してうるう年も計算する部分を入れてますが ここは未使用なのでコメント化してます。

/**
 * 指定の年月日より、1月1日からの累計日数を取得。
 * @param[in]  year    年(西暦)
 * @param[in]  month   月
 * @param[in]  day     日
 * @return 累計日数
 */
int CEarthUtil::m_CalcDays(const int year, const int month, const int day)
{
    int rDays, uruu;

    // うるう年の計算
    uruu = 0;
//  if(!(year & 3) && (((year % 100) != 0) || ((year % 400) == 0))) uruu = 1;

    rDays = 0;
    if(month > 1) rDays += 31;
    if(month > 2) rDays += 28 + uruu;
    if(month > 3) rDays += 31;
    if(month > 4) rDays += 30;
    if(month > 5) rDays += 31;
    if(month > 6) rDays += 30;
    if(month > 7) rDays += 31;
    if(month > 8) rDays += 31;
    if(month > 9) rDays += 30;
    if(month > 10) rDays += 31;
    if(month > 11) rDays += 30;

    rDays += day;

    return rDays;
}

春分の日の3/20は「m_CalcDays(year, 3, 20)」で1/1からの日数を計算できます(yearは未使用なので適当に入れる)。

指定時刻を00:00:00からの経過秒数に変換

これは説明不要ですね。1時間は60分、1分は60秒、を単位変換してます。

/**
 * 指定の時刻より、00:00:00からの累計秒数を取得
 * @param[in]  hour    時間
 * @param[in]  minutes 分
 * @param[in]  second  秒
 * @return 累計秒数
 */
int CEarthUtil::m_CalcSec(const int hour, const int minutes, const int second)
{
    return (hour * 60 * 60) + (minutes * 60) + second;
}

1日の経過秒数は「24*60*60」、1年の経過秒数は「365*(24*60*60)」。

これに加えて、地軸の傾き23.4、地球が自転にて1回転するときの経過日数0.9973、を使用します。

ということで、続きは次回へ。

続・雲表現のテスト(2009/12/29)

雲表現で自己影をつけてみました(先日は散乱の近似を入れていなかったため、うまくできてなかったのでした)。空を見上げる形だとそれっぽくなりますね。

上記の動画。

Gems5や論文では、空は平面じゃなくて半球の一部になっているのですが後はそれを対応すればとりあえずは、太陽光・背景・雲、でプラグインとして公開できそうです。

とりあえずこれは少し置いておいて、太陽の入射角度計算の解説を入れて今年は締めくくることができればと思ってます。

雲表現のテスト(2009/12/28)

そういえば、GameProgrammingGems5にて雲表現の説明があったなと思ってちょっと掘ってみました。まだボリュームっぽさがないですが、、、。ボリュームレンダリングの要領なのですぐにいけると思ったら、どこかで躓いたのか影がうまくつかなかったので雲の影なしです。

Shadeの背景シェーダーを使った雲を描画するテスト。

上記を動画にしました。雲が流れます。東京の緯度経度にて12/20 09:00〜18:05をシミュレート。雲は物理的に正確なものじゃなくて、ウソっぱちです。

夕方あたりでいきなり光源色が変化してうそっぽいですが(^_^;;。ここはバグってるっぽい。

背景描画で使用している「sunsky」で計算されたRGBは、輝度(露出)を2^(-14)にするとそれらしい色合いになりました。Shadeで露出がかなり大きいものをレンダリングして試していたのですが、白飛びした結果から画像処理として露出調整しようと思ったら手段がないですね。このへんは本体機能もしくはプラグインでもいいのでほしいところ(ちなみにプラグインだと作れますよ)。で、GIレンダリングすると、背景での光量にかかわらずレンダリング結果は後で露出調整でいじくりまわせると記憶していたのですが、Shadeの場合は影が吹っ飛んでしまいますね。露出を下げたら影はしっかり見えていると思うのですが、どうだっただろうか?後で、別レンダラで確認してみようかな。

また、太陽が背景のIBLとして焼きついていれば平行光源(Shadeでは無限遠光源)はいらないと考えているのですが、これも実験がいりそうです。今のところ、Shadeでは光源としてIBL(背景のhdr)に極端に強い要素を入れてレンダリングした場合に ハードエッジな強い影が出ないようです。

露出調整は、レンダリング結果をhdr形式で出力してHDRViewで調整するのが今のところ楽な感じがします。

http://www.debevec.org/FiatLux/hdrview/

話を戻しまして、雲はPerlinNoiseを使ったもので古典的なやつです。今はGPUでリアルタイムにもっとリアルな雲を表現してますね。

実は、ニコニコ動画でたまたまベヨネッタの動画を見ていて、背景すごいなぁと思って感化されたんですよ。このゲーム、大神を作ったチームが開発してるんですね。音楽がすごいいいなぁ。そういえば、大神伝という大神の続編っぽいものがNDSで出るらしい、というのをつい最近知りました。これは自分の中では購入決定かもしれません。

太陽を考慮した天空での背景を計算する(2009/12/27)

太陽を考慮した背景色計算についてナイスな情報を教えていただきました。

藤田さんのサイトより。

http://lucille.atso-net.jp/blog/?p=115

リンク先の「sunsky」にソースがあります。

http://www.cs.cmu.edu/afs/cs/user/ajw/www/software/

藤田さんのサイトにて、オリジナルをさらに噛み砕いたソースがありましたので、それを使ってさっそくShadeの太陽光計算プラグインにて背景色をつける実装をしてみました。

東京駅の位置にて12/20 16:00。

12/20の08:00〜17:05までのアニメーション。Shadeの背景シェーダー機能を使っています。

背景だけでなく、太陽光の向き計算も上記ソースで可能だったので私の計算と照合してみると、だいたい同じでした(微妙に差があるのですがこれは春分の日の基準とかが異なるからかも)。

ですが、Shade11ではシェーダー部分で不都合があるのか、レンダリングにてマルチスレッド使用時でかつマルチコアの環境だと、ブロック状のノイズが発生してました。1スレッドにすると問題なし。Shadeの場合はアニメーション(シェーダー対応)というよりも、静止画としてhdr吐き出しモードを設けたほうが賢明かな。そうすると、他のレンダラに持っていけますし。

それと、sunskyのほうですが太陽が沈む瞬間にて夕方の時間で太陽色が突然RGB(0, 0, 0)になるのですが、夕方から夜の移行はどうするのだろう?なんとなく、トーンマップのレンジを調整してあげる必要がありそうな気もしてます。というか、人間の目はたしか自動でトーンマップして太陽からの明るさを抑えたりしていたと思いますが、それの実装が必要なのかな。

内容を理解した上で実装したいので、いずれにせよ論文そのままから一歩先に進んで スペクトルあたりや散乱関連含めて掘ってみよう。

ところで、いじっていると太陽光を考慮した青空や夕焼けだけじゃなんとなく物足りなく感じてしまいました。欲を言えば雲がほしいかもしれない(^_^;;。「ミー散乱」も考慮した計算もやりたいですね。

地球の自転と公転(2009/12/26)

地球の自転

地球の自転は、「地軸」を中心に約1日で1回転する運動になります。「地軸」は北極点と南極点を結んだ直線で、地球を上から下にずばっとだんごのように突き刺さったものと考えてください。

ここで2点条件があります。1つめは「地軸」は傾いているということ。地球の緯度0・経度0の地点を中心として太陽から真正面に地球を見たとしましょう。

このとき、時計回りに23.4度傾いています。ここで基準となるスタート地点となるのは緯度経度が0になる点です。つまりは上記画像の見え方をしているのが基本になります。

2つめの条件は約1日で1回転です。「約1日」なんです。地球が一回転するときの経過日数は1日でなく、約0.9973日です。この微妙なずれがものすごい大事です。「地球は1日に1回転するんだよ」という回答は間違いですのでご注意を。

これを秒に直すと

86166.72 = (24 * 60 * 60) * 0.9973

です。1日は24時間で1時間は60分、1分は60秒、ですので上記計算になります。

自転のキーワード数値は、23.4度の傾きと、0.9973日/回転、の2つ。後々、これを使うときにいかに重要か説明しますね。

で、自転というのは経度の回転と同じだなと気付くかと思います(経度も反時計まわりでプラス方向の回転ですので)。

太陽周りの公転

では次は公転です。太陽の周りを地球は回っていますがこれは1年(365日)で一周とします。上で説明した自転の時の図からずっとカメラをさげていくと以下のようになります。

「春分」の位置が基準となるスタート地点とします。

このとき、地軸の傾きはどの位置でも変わりません。傾いたまま太陽周りを一周すると考えていただければ。なお、実際の地球と太陽の距離はもっと離れており、太陽は地球に比べてかなり大きいです。ですので、上記は説明のための便宜上の図と考えてください。

上から見ると以下の感じ。

春分から開始して90度回転ごとに「夏至」「秋分」「冬至」となります。

そして日にち的には以下のようになります。
---
春分03/20
夏至06/21
秋分09/21
冬至12/21
実際は、年によって1日ずれる場合もあるようです。まぁ、だいたいで。

地球は地軸を中心に自転しながら、太陽を中心に公転します。両方とも反時計回りの回転です。

で、プログラムの説明ですので基準がないとどう回転すればよいか分からない、ということで「春分の日(03/20)のUTC 12:00:00の緯度0・経度0のポイントでは、頭上に太陽が来る」という条件を利用します。スタート地点は、03/20 12:00:00からです。計算で単位が複数あると面倒ですので秒で考えていきましょう。この日時からの経過秒数より、まずは緯度経度を指定している場合の地球の自転による位置はどう変化したか、を計算することになります。

ちょうど頭上に太陽が来る(このときは影がちょうど真下に出るため、地面に垂直に刺した棒の影は見えないことになります)というのは日本ではないですが、UTCでの時間指定で 03/20 12:00:00の緯度0・経度0の地点に もし船をこいでいけたとすると、ここが世界一暑い場所ってことになりそうですね。先日作った太陽光計算プラグインで計算させると、ほぼ真上から太陽が入ることになると思います。気温も計算できると面白いだろうなぁ、はたして何度なんだろうか、、。私、昔グアム行ったときもたまらんかったので(暑さは得意なのですが、日焼けがすごかった)それよりも暑いってことだよなぁ。

とりあえずの必要な情報は以上かな。次回は、これらを踏まえた上での計算に入ります。

追伸

上記とは直接的には関係ないですが、空色の計算は難しい、ですがなんとなく見えてきました。キーワードは「レイリー散乱」と「ミー散乱」です。波長ごとに分解された色は大気の厚みにより吸収・散乱するのですが、「空はなぜ青いのか」「夕焼けは赤いのか」などはこれで説明が付くみたいです。

緯度と経度(2009/12/26)

まずは基本的な緯度と経度について。地図としても正しくなるようにしたいので、DirectXのサンプルテクスチャを使用した図を作ってみました。

緯度経度は極座標で考えます。地球の中心が原点です。緯度は赤道周りを0度とみなしたときの北極へ向かう垂直方向の回転を行い、これを北緯としています。北緯90度だと北極になりますね。赤道から南極に向けての回転は南緯となります。ここでは-90度を南極点とします。つまり、緯度は-90度〜+90度になります。上図では青色の矢印です。

経度は水平方向の回転で、ロンドンのグリニッジ天文台が基準になるというのはご存知のとおり。ここを基準に東方向に回ると東経、西方向に回ると西経になります。東経方向へは0度から180度、西経方向は0度から-180度とあらわせます。ようはぐるっと一周で、0度〜360度(東経 0〜180、西経 -180〜0、で一周)になりますね。上図では赤色の矢印です。

そして、緯度経度がともに0になる位置は、アフリカ大陸の西の海上にあることにあります。地球の自転ではこの基準をスタート地点として行います。地球の自転は北極を見下ろす上から見て反時計周りに回転します(上図でいえば、東経向きとなる赤色の矢印方向)。

水平方向の反時計周りの回転をφ・垂直方向の反時計周りの回転をθ、とすると、緯度θ・経度φ、で地球上の位置があらわされることになります。

これ、まんま極座標ですね。

次回は地球の自転と太陽周りの公転について、を記載します。年末までは、指定の日時・場所における太陽入射角度の計算を紐解いていきましょう。

これらの計算は「基準」を定めることが大事で、地球の自転・太陽を中心とした公転・緯度経度による回転、はそれぞれ絡み合っている関係になります。これに、加えて日時も考慮する必要がありますので複雑怪奇に思えますが順番を追っていけばそんな難しくはありません。

ということで、今年は残すところ一週間ですね。

太陽光計算 for Shade10/11(2009/12/25)

日時と場所(緯度・経度)を指定することで太陽光の向きを計算し、Shadeの無限遠光源に割り当てることができるプラグインを「太陽光計算」にてアップしました。Shade 10/11のWin32/64ビット、Macの32/64ビットのstandard/professionalにて動作します。

いつもどおりソースコード付きです。初回バージョンですのでたいしたことはしておらず、照射方向を当てはめるだけです。

本当は、背景をシェーダー的に青空や夕焼けに変化させる処理を入れたかったのですが、波長計算なども考慮する必要があり 思ったよりも難しかったので見送りました(残骸がソースに残ってますが今は触れないでください(汗))。

この日記で解説する教材になるもので「地球の自転と太陽周りの公転はどう計算するの?」の説明で使う副産物です。今は、日時や緯度経度を入れて遊んでみてください。太陽は東から昇って西に沈む、夏と冬の日の沈む時間の違い、などがシミュレートできますので実際の世界と見比べてみるのも面白いかと。

Shadeでの背景シェーダーですが、実はプラグインSDKではすでに実装されており(Shade9あたりだったかなぁ)、GIレンダラ(パストレ)でもここがIBLで反映されるようになってます。これを利用しているShade本体の機能や外部プラグインは皆無だと思いますが、実はSDKレベルでは結構なことができるんですよ。これにタッチして背景を青空にしたり夕方は夕焼けにしたいところではありますね。

ですが、GPSでの衛星軌道計算を先にしたいのでバージョンアップは来年になるかなぁ。

mqoImporter/Exporter Ver.1.3.0 for Shade10/11(2009/12/23)

メタセコイア形式のmqoファイルをShadeとやりとりするmqoImporter/ExporterプラグインをVer.1.3.0に更新しました。Shade 11にて動作するように修正しています。Shade 10でも動作することを確認しました。

http://ft-lab.ne.jp/cgi-bin/wiki.cgi?page=mqoImporter_shade

なお、Mac版はXCode 3.2にてビルドしています。ソース類も同梱していますので、プラグイン作成の勉強用にご利用くださいませ。

続・mqoImporter/Exporter Shade11での動作について(2009/12/22)

Shade11にてmqoエクスポート時にテクスチャが出力されないという不都合について。今までエクスポート時にファイルを出力する場合に、たとえば画像の場合は「image_interface::save("test.png");」のようにするとShade10.5以前はエクスポートファイルと同じところがカレントディレクトリになり、そのディレクトリ内にファイル出力できていたのですが、Shade11ではカレントが変化してしまっている、というのが原因のようでした。

なので、この部分をフルパスにして出力するとうまくいってます。ほかに不都合が出ないか確認してからShade11対応として公開しますね。

同等のエクスポートプラグインを作っている方は注意が必要です。1ファイルのみ出力ならいいのですが、同時にテクスチャを出すという動きの場合は同じ現象が起こるかもしれません。出力時に絶対パスで出していればよいですが、インポート・エクスポート時の絶対パスの取得方法はShade9〜10でたしか変わってしまってたので過去バージョンにさかのぼってのチェックもいるかも。

Shade11ではdo_exportの引数より「plugin_exporter_interface::get_file_path」で出力ファイルのフルパスを取得できますが、「plugin_exporter_interface::get_text_stream_interface」で取得したtext_stream_interfaceでの「get_name」ではパスは取得できてません(これはカラ文字が返される)。stream_interfaceからはインポータと同等ならパスが取得できたと記憶していますが、未確認です。

mqoImporter/Exporter Shade11での動作について(2009/12/21)

mqoExporter for ShadeのShade11動作ですが、1点問題を見つけました。

Windows32ビット環境にてmqoエクスポート時にテクスチャが出力されていない(もしくは別ディレクトリに出力されている?)ようです。一度デスクトップにテクスチャが展開されたことがあるので、もしかしたら出力されていないのではなくてどこか別ディレクトリに吐き出されているのかもしれません。Shade10ではもちろん動作していてShade11に持ってきて起こっていた部分ですので、回避できそうかチェックするようにします。

やっぱしバージョンアップでの挙動が変わるのは避けられないですねぇ。簡単にしか動作確認しませんでしたが、もっときっちりやらないと怖いね。

続・太陽光をシミュレートする(2009/12/19)

下の日記に書いた太陽光計算のムービーを作ってみました。

位置は東京駅を想定 : 12月19日 08:00 - 17:15。15分刻みに1フレームを進めています。

位置は東京駅を想定 : 8月19日 08:00 - 18:30。15分刻みに1フレームを進めています。

で、あくまでも計算が正しいと仮定した場合のシミュレーション結果ですが、12/19の日付では17時過ぎには日が落ちる、08/19の日付では18:30ごとに日が落ちる、で明らかに差がありますね。

これが沖縄だと冬でももう少し日が落ちる時間が遅いです。

まだテスト的なプラグインですので、もうちょっとまともに使えるようになったら勉強用にソースごと公開予定です。

太陽光をシミュレートする(2009/12/19)

地球の自転と太陽周りの公転を考慮すると、「指定の日時」「指定の緯度経度」の情報を与えると、どの方向から太陽光が入ってくるのか、というのが判断できます。

Shade11(SDKはShade10のもの)で太陽光をシミュレーションするプラグインを作って実験してみました。

以下は、1月1日9時/東京駅での太陽の当たり具合をレンダリングしたもの。

以下は、1月1日12時/東京駅での太陽の当たり具合をレンダリングしたもの。

以下は、8月1日9時/東京駅での太陽の当たり具合をレンダリングしたもの。

で、時間によって日時によって太陽の当たりっぷりが違うことが分かります。計算があってるかの確認のため、GoogleSketchUpの影計算と照らし合わせてみました。スケッチアップはデフォルトでの位置は東京駅かな(位置の変更はできないのでこのあたりは不明)。まったく同じ感じになりましたのでこれでOKかと。

この計算を日記ネタにしてみようと。最終目的は衛星軌道を求めることで、地球の自転は考慮しなければならないためこれは副産物となります。

これらの照射計算に関しては「地球の自転は24時間で1回転ではない」「地軸の傾きが23.4度」というのを元に紐解くとそんなに手間ではないです。

後は、夕焼けはなぜ赤いのか、などの情報を入れたいけどまだ調べてないです。Shadeではエフェクタプラグインで背景となる色を変更できるので、それ使えばレンダラでも反映することができるかもしれませんね。これも確認してみようかと。

iPhone Developer Program(2009/12/17)

米Appleからメールが来ました。で、作業を進めていくと日本のAppleStoreでiPhoneの開発ライセンスを購入できました。なんだ、入り口だけ英語でショッピングは日本語でOKなんだ、安心しました。

ということで、おそらく年末最後の散財としてiPhoneの開発ができるように。後で分かったのですが、iPhoneエミュレータと書いてたのですがシミュレータの立場なんですね。

それと、料金を払って使えるのは開発したプログラム(ソフト)をiPhone(ハード)に持っていくためのライセンス料みたいな感じです。シミュレータで遊ぶ分には無料で使えるみたい。ただ、GPSやカメラなどをいじりたいのでいずれにせよ私の場合はiPhone自身を購入しないと先に進めないです(<あくまでも遊びで)。

その前に、先に手持ちのGPSで衛星軌道を求めてみたいので今年の実験はこれがらみで締めくくりかなぁ。

Shade 11(2009/12/17)

発売日より一足早く、Shade11が届きました。注文したのが4日前だけど(^_^;;。

で、私が作ってたプラグインが動作するか確認。mqoImporter/Exporterは問題なく動いているようです。

発売記念でまたプラグイン作ろうと思ってたのですが、明日には間に合わないな、、、。予定していたのがあるので、おそらく今年中に公開すると思います。当分はShade10プラグインSDKで問題なさそうですので、しばらくはこれで。

誰か作ってるかもしれないけど、日付時間と緯度経度を入れての太陽光シミュレートがほしいよね。今までの日記の話題でうすうす気付いている方もいたかと思いますが、副産物でこれができるなぁと。

Google Mapsのルート検索(2009/12/10)

またまたGoogleが新しいサービスを公開しましたね。徒歩のルート検索です。

http://maps.google.co.jp/

ためしに、「品川駅」から「吉野家 品川駅港南口店」でやってみると

な感じでルートが示され、「改札口に向かう」「屋内通路を進む」「歩道橋を渡る」「下り エスカレーターに乗る」「右折する」「横断歩道を渡る(目的地は前方左側です)」なんて案内も出てきてます。

これはすげ〜〜。散歩好き(かつ、方向音痴なのに迷うのが好き)な自分としてはありがたいなぁと。iPhoneとかの携帯端末でも対応しているとのこと。

また、最寄駅から自宅までを徒歩ルートすれば(結構距離があります)、結構それなりの近道が出てきました。

ちなみにGoogle Maps APIでは、GDirectionsがルート検索のクラスのようです(徒歩も対応)。

しかし、既存のナビサービスはこれが出てしまうと商売上がったりになりそうではありますね。なんせ無料で使えるわけですから。

Googleはなんだか最近どんどん目新しいサービスを無償で発表しているのですが、何かあったのだろうか、、、。個人的には次に何を発表するのだろう?というのが楽しみではありますが、見方によってはクラッシャーっすね。ナビ系サービスはいろいろ淘汰されそう、、、。

iPhone SDK(2009/12/09)

iPhone 3GSのスペックを眺めていたらほしくなってしまい、つい、iPhone SDKのデベロッパプログラムに登録してしまいました。といっても、まだ手続き中なんですが。昔、XCodeの最新版がADC(AppleDevelopperConnection)の入会者しかダウンロードできないときがあり、そのときの情報でやりとりできました。ADCに入っていないとiPhone SDKはゲットできません。

開発にはお金がかかりまして、スタンダードが10800円です(でも、エミュレータで遊ぶ分には無料でできるっぽい)。IntelMacでしか開発は出来ません。入会手続きで英語サイトに飛んだのですが適当に入力。住所を英語で書くなんてわかんね〜よ、というか今のアパートの正しい英語綴りが分からん(^_^;;。たぶん後々クレジット入力のメール案内が来ると思うので、それだとなんとかなりそうだ。

最終確認後のショッピングカート画面が英語で、日本のAppleのカートは別なのかな?日本のAppleStoreのアカウントも持っていて、ログインして確認したのですがこちらには情報がないし。というか、日本のApple、このあたり日本語でやりとりできるようにしてくれ、と切にお願いしたいです(^_^;;。英語が苦手であきらめモードの自分としては、ものすごい敷居が高いんですけど、、、。

ということで、来年頭あたりにiPhoneは購入する予定で考えています。今はソフトバンク携帯なんですが、GPS機能がなかったものなのでいい機会かなと。ちょうど二年前に購入したものなんですが、思えば携帯は進化してますねぇ。iPhone 3GSでは、GPSはA-GPS対応、加速度センサー、電子コンパスまでついているみたいで、開発プラットフォームとしてはかなり魅力的だったので。後、ゲームも遊びたいしね。1月末まではキャンペーン中かで安いんですよね。

iPhone SDKはエミュレータがついているはずですので、とりあえずはこれで遊んでおこう。

、、、と思って、ADCでいろいろたどってたらSDKはなんの登録もなしでダウンロードできました。料金の支払いはまた別なのかな?

そして、なんの手続きもなくエミュレータで遊べますね。

日本語表示してみました(Cocoaでプログラム)。

Google Goggles(2009/12/08)

Googleの「Google Goggles」。

http://www.itmedia.co.jp/enterprise/articles/0912/08/news027.html

携帯で撮った写真画像を解析して、そこにタグ付けされた情報を出すという技術。「東のエデン」の第一歩がこれに近いんですが、思ったよりも早く世に出てきそうすね。

どうやればこれが実現できるか、は、実は以前掘っていたのでだいたいはつかめているのですが、先を越された〜〜。しかし、最近はGoogleは快進撃ですねぇ。なんか急にいろいろ発表しだしてるし。

GPSの位置判定(2009/12/07)

Wikipediaを見ると、衛星位置、衛星内の時計での時刻と測定位置での時刻との差分から求まる距離で、測定位置(緯度、経度、高度)を推定するという主旨の記載がありました。方向から求まるんじゃないのね、昨日の書き込みは嘘書いてました、すみません。

http://ja.wikipedia.org/wiki/%E3%82%B0%E3%83%AD%E3%83%BC%E3%83%90%E3%83%AB%E3%83%BB%E3%83%9D%E3%82%B8%E3%82%B7%E3%83%A7%E3%83%8B%E3%83%B3%E3%82%B0%E3%83%BB%E3%82%B7%E3%82%B9%E3%83%86%E3%83%A0

後、ここのサイトの解説もプログラムが書いていて分かりやすかったです。

http://www.enri.go.jp/~fks442/K_MUSEN/

NMEA-0183には衛星位置や距離情報はなかったのですが、ここはすでに調理された後のものが入ってる、ってことになるのかな。

となるとすると、取得できる観測位置から見た場合の衛星の方向情報からさらに推定して衛星位置を求めるというのは 大変かもしれない(誤差も含めて)。それよりも、距離(または衛星と測定位置での時間差分)に相当するパラメータがないと位置を推測する材料が足りないですね。測定位置から見た場合に、どの方向に衛星がいるのかはパラメータとして存在するのでこれをとっかかりになんとか3Dの地球儀の上に旋回させてみたいところですがさてはて。う〜む、届きそうで届かないか、まだ判断できてません。地球の自転もあわせて調べていたら結構こんがらがってきました。

根が深そうなので、また土日に掘ってみよう。

地球を回っている衛星の数は?(2009/12/06)

話がそれていってしまいますが、、、地球を回っている衛星の数はどれくらいあるのでしょう?

http://www.jaxa.jp/pr/inquiries/qa/satellite.html#01

より、3000個以上あるらしい。

NASAの以下にて、地球の周りを回っている衛星軌道を3D化したものをJavaで見ることができます。

http://science.nasa.gov/Realtime/jtrack/3d/JTrack3d.html

これは900個以上の衛星を視覚化してるとのことですので、実際はこれの3倍はあるということになりますね。よく衝突しないもんだ、、、。

白い点は星ではなくて全部人工衛星です。GPSの衛星番号08を発見。

地球の大気圏すぐのところに固まっているのと、はずれにあるものとが多いですね。外にある大きい輪はなぜか赤道に沿って平面的です。赤道に沿って周回する衛星が多いのか、これは面白い。

衛星の情報を取得する(2009/12/06)

GPSで取得できること(NMEA-0183プロトコルに記載されていること)として、衛星情報があります。「SkyTraq」というツールにて衛星データを見たものは以下のように。

番号のついた円が衛星をあらわします。番号は衛星番号で、全世界の打ち上げ衛星で通し番号がある感じかな。地球があって、その上空に衛星をマッピングできると。青・緑・赤の円の色の意味ですが、これは信号に対するノイズ量を対数であらわしたもので「SNR」と呼ばれるものです(単位はdB(デシベル))。この値が大きいほど高品質(ノイズが少ない)ということに。

SNR値は以下のように表示されていました。

このときのNMEA-0183での記述は以下のようになっています。

$GPGSA,A,3,25,10,02,04,05,07,08,29,15,,,,1.8,0.9,1.5*33
$GPGSV,3,1,11,02,80,148,25,10,64,341,26,05,55,314,33,04,40,136,18*74
$GPGSV,3,2,11,08,35,107,26,07,28,065,21,15,26,226,25,29,22,305,12*76
$GPGSV,3,3,11,24,20,314,,25,19,047,20,13,10,042,14*45

「$GPGSA」「$GPGSV」が衛星に関する情報が記述された行になります。

詳しくは以下のページなどを参考に。

http://bg66.soc.i.kyoto-u.ac.jp/forestgps/nmea.html

「$GPGSA」より、受信衛星数は「25,10,02,04,05,07,08,29,15,,,,」より9個分。この「9」の情報は「$GPGGA」での衛星数と同じですね。9個の衛星から今のGPSレシーバの位置を推定している、ということになりそうです。

そして、それぞれの番号が有効な衛星番号です。この場合は、「25/10/02/04/05/07/08/29/15」。

これらの衛星ごとの位置情報やNSR情報は複数列挙された「$GPGSV」より解析できます。

$GPGSV,3,1,11,〜 
$GPGSV,3,2,11,〜 
$GPGSV,3,3,11,〜 

の3つの場合は、それぞれに対して3つのメッセージ($GPGSVの出力)を持ち、2つめの連番で連続していることをあらわします。1つのメッセージに対して4つの衛星情報しか記載できないため、4+4+3=11の衛星情報を並べています。

「$GPGSV」が分かりやすくなるように改行します。

$GPGSV,3,1,11,
   02,80,148,25,  ←衛星番号02(仰角 80度 /  方位 148度 / SNR 25dB)
   10,64,341,26,  ←衛星番号10
   05,55,314,33,  ←衛星番号05
   04,40,136,18   ←衛星番号04
   *74
$GPGSV,3,2,11,
   08,35,107,26,  ←衛星番号08
   07,28,065,21,  ←衛星番号07
   15,26,226,25,  ←衛星番号15
   29,22,305,12   ←衛星番号29
   *76
$GPGSV,3,3,11,
   24,20,314,,   ←衛星番号24(NSRが空白なので無効)
   25,19,047,20, ←衛星番号25
   13,10,042,14  ←衛星番号13
   *45

これと上の画像の衛星番号/SNRが一致しているのが分かります。位置は「仰角/方位」でマッピングできそうですね。ですが、整数の度数指定なので精度は期待できずだいたいの位置になりそう。「$GPGSA」の列挙では衛星番号13はないため、これは欠番(SNR値が低いです、29番も低いですがなぜこれが有効なのだろう)。

ただ、おそらくですが上記SkyTraqのレーダーみたいな地球の画像は飾りで、
実際の仰角/方位から求められた絶対的な衛星位置とは
一致していないものだと思います。
ドラゴンボールのドラゴンレーダーみたいなもので、「仰角/方位」は今の位置を
中心として北を垂直上方向に見立てて「方位」で0〜360度回転、
「仰角」で0〜90度(90度は観測位置の真上)、としているかと。 
ようは、観測点から見たときに
半球上のどの方角から衛星の情報が入ってきているのか、なので
これに加えて、各衛星が何時何分何秒にどの位置にいるのかというのを
(これは決まっているかと思います、同じルートを周回してるので)
入れてあげないといけないのかなと。
ここで3D野郎は「半球!?」で何かを思い出しそうですね。
このあたり、かなり3DCGに近い世界です。
パストレーシング(もしくはフォトンマップ)か?と
思ったあなたはおそらく正しいです(^_^;;。
サンプリングで精度を上げているというのがGPSの肝っすね。
なぜ、GPSで位置が推定できるのか、という根本的な部分も、
この入射方向に由来するものですのでいずれ紐解いていこうと思います。

ということで、Google Mapsには自分の位置である緯度経度指定する(これは「$GPGGA」を使用する)、ということにはなりますが、地球規模で言えば衛星データもマッピング可能です。精度が悪い場合は、たいていは受信に使用した衛星数が少ないかNSRが低いという可能性がありますので、これに閾値を設けてマッピングにフィードバックさせることで多少はぶれを抑えられるかも。ロガーでハード的に保存される情報はこのあたりはチェックしてない感じがします。gpx出力ではこれを出してくれないため、ロガーの情報はちょっと情報不足ですね。

追記:

以下のページでPRNとあらわされるものが衛星番号になるようです。

http://www.navcen.uscg.gov/navinfo/Gps/ActiveNanu.aspx

番号は01〜32で増えることはないらしい(GPS衛星の場合)。これでGPS衛星は地球全部をカバーできているとのこと。

ただ、衛星は寿命があり設計寿命が7.5年、寿命が尽きたGPS衛星は破棄され(打ち落とすんじゃなく、軌道をずらして落下を待つ)、新しい衛星を打ち上げた際にPRN番号はそのままで世代交代してるっぽい。ゴミとなった衛星はやはり衝突などの問題になるようですね。打ち落とすとさらに被害が出る、というか過去にそれをやってゴミを撒き散らしたことがあるらしいので、大きな課題だそうで。

ひさびさにすれ違い通信(2009/12/06)

もう長いことドラクエ9は触ってなかったのですが、ひさびさにすれ違い通信がいまだに健在なのか調べてみたくなって、DSiをかばんに入れて近所を散歩してみました。

結果としてはすぐに3枠埋まりました。まだ健在のようです。でも、長時間プレイしている方ばかりで300〜400時間台ばかり。私は190時間ほど、です。休日か寝る前にコツコツ時間を使ってたんですが、さすがにもう飽きたな〜〜。

配信クエストもひさびさに取得しに行ったら5イベントくらい増えました。ということは一ヶ月ちょっと触ってなかったってことになりますかね。ですが、寝る前の数時間で消化。配信クエストがすべて終わるまではちょくちょくやっていこう。他のゲームを買う予定もないですし。後、個人的にはドラクエ9でのすれ違い通信がいつ廃れてしまうのか、も知りたいかも。

gpxフォーマット(2009/12/05)

GPSロガーで出力されるgpxフォーマットについて。

http://www.topografix.com/gpx.asp

にある

http://www.topografix.com/GPX/1/1/gpx.xsd

がスキーマの記述で何が書けるのかはここを見れば理解できます。

なお、GT-730F/Lに付属の「GPS Photo Tagger」でGPSレシーバから読み込んで出力したgpxファイルを元に解析してみます。緯度経度は実際にGPSに記録されたやつですが、個人情報(ようは自宅位置)が出ないように遠隔の地での計測値を入れています。

以下、gpxファイルの一部を取り出したもの。XML形式です。

<?xml version="1.0" encoding="UTF-8"?>
<gpx version="v1.1.5"
    creator="iTravel Tech Inc. - http://www.itravel-tech.com"
    xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
    xmlns="http://www.topografix.com/GPX/1/1"
    xsi:schemaLocation="http://www.topografix.com/GPX/1/1
        http://www.topografix.com/GPX/1/1/gpx.xsd
        http://www.topografix.com/GPX/gpx_overlay/0/3
        http://www.topografix.com/GPX/gpx_overlay/0/3/gpx_overlay.xsd
        http://www.topografix.com/GPX/gpx_modified/0/1
        http://www.topografix.com/GPX/gpx_modified/0/1/gpx_modified.xsd">
       
<trk>
    <name>Track2009/12/05_13:13</name>
    <trkseg>
        <trkpt lat="35.853715" lon="139.681340">
            <ele>39.213059</ele>
            <time>2009-12-05T05:20:42Z</time>
            <desc>lat.=35.853715, lon.=139.681340, Alt.=39.213059m. Speed=1.000000m/h.</desc>
            <speed>0.277778</speed>
        </trkpt>
        <trkpt lat="35.853739" lon="139.681313">
            <ele>35.598709</ele>
            <time>2009-12-05T05:20:47Z</time>
            <desc>lat.=35.853739, lon.=139.681313, Alt.=35.598709m. Speed=3.000000m/h.</desc>
            <speed>0.833333</speed>
        </trkpt>
    </trkseg>
</trk>

</gpx>

「<?xml version="1.0" encoding="UTF-8"?>」はおなじみのXMLタグです。「gpx」タグにてcreatorがgpxを出力した製作者情報、ほかは決めうちでよさそうです。「xsi:schemaLocation」は、行を分けてますが実際は1行につながってます。

Wikipediaによると

http://ja.wikipedia.org/wiki/GPX

以下のように緯度/経度/高度が最小限の情報でtrkpt(トラックポイント)が複数並ぶことで経路をあらわせます。

<trk>
    <trkseg>
        <trkpt lat="緯度" lon="経度">
            <ele>高度</ele>
        </trkpt>
        <trkpt lat="緯度" lon="経度">
            <ele>高度</ele>
        </trkpt>
    </trkseg>
</trk>

これ以外のタグは、出力するメーカーごとに異なる模様。ですが、計測時間がないと判断ができないため日付時間のtimeタグは普通に追加してそうな感じがします(他のgpx出力を試したわけじゃないので不明ですが)。

「GPS Photo Tagger」では、以下の解釈になりそうです。

<trk>
    <name>トラッキング時の名称(任意)</name>
    <trkseg>
        <trkpt lat="緯度" lon="経度">
            <ele>高度</ele>
            <time>2009-12-05T05:20:42Z</time>
            <desc>説明</desc>
            <speed>速度?</speed>
        </trkpt>
        <trkpt lat="緯度" lon="経度">
            <ele>高度</ele>
            <time>2009-12-05T05:20:47Z</time>
            <desc>説明</desc>
            <speed>速度?</speed>
        </trkpt>
    </trkseg>
</trk>

timeタグは「YYYY-MM-DDTHH:MM:SSZ」の形式で「T」「Z」は固定で、日付と時間指定をしています。時間はUTCで指定。時間に+9すると、日本での時間になります。

で、とりあえず、緯度/経度/高度/日付時間があれば使えそう。speedタグはXMLスキーマには宣言されていないですので独自かな。

これを後はGoogle Mapsに流し込むだけですので(緯度経度だけ)、後は簡単です。

また、GPSでのリアルタイム記録を指定の秒数間隔ごとに蓄えておき、記録後に独自ソフトでgpx形式に出力するとGPSロガーソフトウェアになりますね。ハードに蓄積されるログ情報の引き出し方は分かりませんでした。これもプロトコルが統一されていればよいのですが、、、。

今はGPSをネタに日記を書いていますが、最終的にはこれが3Dとつながる予定です。GPSから得られる情報では絶対位置に加えて高度があるわけですので、これをワールド座標を作って当てはめると、、、仮想空間と現実をリンクできる第一歩になるんじゃないかなぁと(座標の単位は何になるのだろう、、、mmやmじゃないよなぁ、緯度経度から変換できるのかな)。

1つ前の日記に書いたA-GPS/D-GPSでの補間ですが、これは複数のサンプリングの平均を取って精度を絞っていくのにもつながるのかなと思ってます。この部分、三次元的な補間に感じます(という意味合いも込めて、3D空間にGPSの結果を当てはめたい)。

基地局も仮に衛星的な役割をするんだったら、GPSのプロトコルでは(衛星データも記録されるので)どのように出力されているのだろう?これも調べてみるかな。ロガーからの情報ではこれを得られないので、実際にPCを持ち運んで調査しないといかんですが。

で、足りない情報がまだまだあるのですよ。GPS精度の問題もたしかにあるのですが、もっと大事なあるものが欠けています。これはGPSでも無線通信でも電子コンパス(向きの判断)でも補間できないものです。さて、なんでしょう(謎かけみたいですが)?これが今の個人研究の最大のテーマかもしれません。

ということで、これも追々日記に書いていきます。

続・GPSの記録実験(2009/12/05)

GPSのログを取りにリベンジしてきました。今度は、

  • GPSが受信できなくなるような室内には入らない(GPSの受信が途絶えるとその後に精度が下がるようでしたので)
  • ビルが少ない、見晴らしのよい道を選ぶ

に注意して散歩&マッピングしました。

赤い線がGPSでの記録をつないだものです。どうでしょうか、ほぼ道に沿えて記録できています。道路の歩道のどちら側にいたのかもだいたい把握でき、あっております。ですので、誤差は道路の幅以下になりそうです(5メートルくらいかな)。

「東京外環自動車道」の下をくぐって歩道橋をわたったのですが、そこもばっちり記録できています(中央緑色のところを横切ってる箇所)。

場所によって大きく差がある感じですね。見晴らしのいいところだと結構いい感じになってます。ビル街でもこれくらいの精度になるように補間できないだろうか。結局は、空が見えやすいというのが大事な条件ですかね。

A-GPSというので精度をあげているらしく(携帯電波を使用)、これはセンターと通信して補正してるみたいです。対してカーナビはD-GPSで(FM放送を使用)、補正していることは同じとのこと。もうちょっと詳しく調べてみるかな。

GPSの衛星情報に加えてA-GPSの携帯電波を使った補正で精度を上げている、ということになりますね。室内でうまくいくのは、実は衛星じゃなくてA-GPSの所業なのかもしれません。

GPSの記録実験(2009/12/01)

今年も後一月です。

自重せずに、GPSのログを取ってきました。本日昼食に散歩もかねて吉野屋に行ったときの道のりをマッピングしています(実際の道のりです)。

道中の途中からGPSをOnにしましたが、ウォーミングアップ時間が結構かかっているようです。旗のマークの部分がGPSの記録が開始されたときで、電源Onから20分くらいというところでしょうか。

水色の線が実際に歩いた道のりです。ピンクの線がGPSにてログとして記録された道のりです。

「GPSの電源On」からスタート。吉野屋までの道のりは「産業道路」といいまして一直線の道路になっています。

吉野屋の位置がGPSの記録位置とは道路はさんだ距離ずれていたのと、吉野屋内でGPSにつながらなかったようです。吉野屋を出た後の軌跡が南に寄っており、これは実際に歩いていないルートです。

ここから西川口の駅までは一直線なのですが、GPSではジグザグ。

西川口から蕨駅までは電車に乗りました。ですが、ログ上では電車内ではGPSは一切つながりませんでした(一直線になっているのは、始点と終点を結んでいるから。精度がいいからではないです)。

蕨駅近辺はビルが立ち並んでいるのですが、ルートとは違うマッピングをしているようです。

その後の道のりはほぼGPSと同じでした。この近辺はビルがなく、見晴らしがいいです。で、開始地点に戻った後に電源Off。

結果として、

  • ウォーミングアップの時間が結構かかっている
  • ビルなどが立ち並ぶところでは精度がかなり落ちる
  • 停止する(GPSが受信できなくなる)と、その後は挙動が怪しい
  • 電車ではGPSが使えなかった

というのが分かった気がします。場所によっては一区画以上ずれますね。

大雑把な道筋でいえば「だいたいあってる」なんですが、たとえばリアルタイムに誰かがGPSをつけて歩いているのをドンピシャで追跡する、という映画さながらなことは難しいかなという気が(そんなことをするつもりはないですが(^_^;;)。

車のカーナビの場合は、GPSだけでなく通信か何かの補正処理も入れているのでそれなりの精度になるのだそうですが、GPS単体ではまだまだ得手不得手が見えるかなぁと。

東京でマッピングしたらもっとぶれるかもしれないですね。最終はマッピングするのが目的ではないのですが、さてはて絶対位置を得る道具として使えるかどうか、いろいろつめていく必要がありそうです。

先月末に書いた「Google Maps Navigation」ですが、日本じゃビルだらけなのでナビとして使うのは難しいのかもしれないですね。今日本にあるGPSを使ったサービスは「近辺の店や最寄り駅を探す」っていう用途でしょうから、道が一区画ずれるとかならナビとしては厳しいでしょうし。

Future's Laboratory 技術格納庫 2004-2013 Yutaka Yoshisaka.