SS1の日記: アホの子が考えた気圧の計算法 2
アホの子が考えた天気図の気圧のテキトーな計算方法
注:以下の話は間違っています。用法,用量を考えて参考にしてください。
以前にやった仕事で,車上で記録されたGPSデータから現地の気圧を推定する必要がでてきて,そのとき考えた話。
求めるデータは,1キロ走行毎の現地気圧。基礎にしたデータはこんな感じ。
1.テレメータログ,GPS緯度・経度,時刻,気温・・・(で,気圧データが無かったので推定する必要に迫られたわけです)。
2.気象台データ 気温,大気圧(高度補正済み) ネットで公開されてます。
3.地図データ(該当地点の標高) これもあったよね。確か。
4.大気圧と標高の関係式(知りたいのは現地気圧だったので。補正式は理科年表に載ってる)
私の考えた,てきとーなアルゴリズムは次の手順
1.該当地点近傍3点を抽出
各気象台の座標と現地座標から絶対距離を求め,その長さでソート。近傍3点をリストアップする。
2.上記3点の気象台データから該当時刻の気圧データを出し,3点に接する平面を求める
地球上で平面というのはアレなんだけど「メルカトル地形図上の」平面ということで。
3.平面の値を求めたら,それを基に現地点の気圧を平面状の大気圧から求めて,標高(気柱高度)の換算をして求める。
こんな感じ。単純な3点からの絶対距離による加重平均にしなかったのは,該当地点が選んだ3点を結ぶ3角形の内側にある保障が無かったから。目的の地点が3角形の内側で無いのに加重平均すると,ちょっとおかしな計算結果になってしまうため。んで,気象台は大抵県庁所在地というか各県の真ん中辺にあるので,海岸部だと3点目がその外側に無くて,そのとき,ヘンな計算にならないような計算方法にしたかったわけです。
たとえば,木更津の気圧を推定するのに,近傍3点に千葉,東京,横浜を選ぶと木更津は3角形の中に納まらない。だからといって東京の代わりに八丈島とかハワイのデータを使えるのか,という問題になる。
で,上記の近傍3点に接する平面上の値とすれば3角形の内側なら加重平均と変わらない結果になるし,外側でも,それの延長上の値に出来る。
・・・という考え。で,ごりごりと二次方程式を解く。
1.近傍三点の抽出
緯度経度(x,y)から絶対距離(d)を求め,近い順のリストを作り近くの3点をリストアップ。
計算式:
d=√(x^2+y^2)
(・・・いや。極座標系にそれはないんだけど。もちろん)
2.気圧平面の計算
平面の式は,こんな感じ?
式1: ax+by+c=z
近傍3点(P, Q, R)の気象台の緯度経度と気圧データ(時間補正は単純なので先にやっとくか前後で計算して最後に補完する)
{
P: x1, y1, z1
Q: x2, y2, z2
R: x3, y3, z3
}
3.上記の連立方程式を解く
んでまあ,頭の良い人ならベクトル演算で解いちゃうわけですが,頭悪いので,シコシコと解いてみた。
4.上記式1を因数分解(って呼び方でいいんだっけ・・・)
式1に各データを代入する
{
fp: ax1+by1+c=z1,
fq: ax2+by2+c=z2,
fr: ax3+by3+c=z3,
}
ごりごり因数分解する。
cを求める式に変換。
{
fpc: c=z1-ax1-by1 --式C
fqc: c=z2-ax2-by2
frc: c=z3-ax3-by3
}
c=fpc=fqc=frcなので,fpc-fqc=0で,cを消す。同様にfpc-frc。
{
fpq: 0 = ( z1-ax1-by1 ) - ( z2-ax2-by2 )
fpr: 0 = ( z1-ax1-by1 ) - ( z3-ax3-by3 )
}
で,a を求める式に変換。
0 = z1 -ax1 -by1 -z2 +ax2 +by2
0 = -ax1 +ax2 -by1 +by2 +z1 -z2
aを左辺にもってって,係数をまとめる。。。
ax1 -ax2 = -by1 +by2 +z1 -z2
a(x1-x2) = b(-y1+y2) + (z1 -z2)
a = {b(-y1+y2) + (z1 -z2)}/(x1-x2)
fprも同様に・・・
a = {b(-y1+y3) + (z1 -z3)}/(x1-x3)
で,上記2式から,bを求める。手順はさっきと一緒。
{
a = {b(-y1+y2) + (z1 -z2)}/(x1-x2) --式A
a = {b(-y1+y3) + (z1 -z3)}/(x1-x3)
}
0 = {b(-y1+y2) + (z1 -z2)}/(x1-x2) - {b(-y1+y3) + (z1 -z3)}/(x1-x3)
= b(-y1+y2)/(x1-x2) + (z1 -z2)/(x1-x2) + b(-y1+y3)/(x1-x3) - (z1-z3)/(x1-x3)
= b(-y1+y2)/(x1-x2) + b(-y1+y3)/(x1-x3) + (z1-z2)/(x1-x2) - (z1-z3)/(x1-x3)
= b{(-y1+y2)/(x1-x2)+(-y1+y3)/(x1-x3)} + {(z1-z2)/(x1-x2)-(z1-z3)/(x1-x3)}
bを移して。。。
b{(-y1+y2)/(x1-x2)+(-y1+y3)/(x1-x3)} = -{(z1-z2)/(x1-x2)-(z1-z3)/(x1-x3)}
= (z1-z2)/(x1-x2)+(z1-z3)/(x1-x3)
bの係数を割って。。。
b = {(z1-z2)/(x1-x2)+(z1-z3)/(x1-x3)} / {(-y1+y2)/(x1-x2)+(-y1+y3)/(x1-x3)} --式B
で,上記の式Bを式C,式Aに代入して変数a,cを求める。
a = {b(-y1+y2) + (z1 -z2)}/(x1-x2) --式A
c = z1-ax1-by1 --式C
5.検算
検算1:
式1(ax+by+c=z)に任意の値abcを入れてダミーの3点を計算。そこから上記式ABCで計算して,任意の値abcが求められるか確かめる。
検算2:
任意の気象台を測定点とし周辺3点の気象台データで近似してみる。
#HP電卓がダンボールの底にあるので検算してません。
6.バグとか・・・
よく考えると,各気象台の風向風速データ(その気象台における気圧勾配などの傾向がわかる)も考慮しないとまともな気圧の推定は出来ないんですが,当時はそこまで思い至らなかったっす。それから緯度の1度と経度の1度がほぼ等距離になるのは赤道上の話なので,これも考慮する必要があります。なお,緯度,経度のおよその距離は理科年表に載ってるのでそれを参照すればいいそうです(昔,NIFTYの某フォーラムで教わった)。とはいえ,これも思い至らなかったので無視した・・・ ダメじゃん(笑)。報告書に「ちゃんと知りたいなら気圧データも取ってね」とか言い訳を書きまくった記憶があります(汗)。
・・・などというような,理論に問題がある上に頭の悪い算数問題は二度とやりたくないので,気象台データから天気図を描くアルゴリズムって誰か知りませんか?
つーか,これで近似したことにするには気象台の気圧データってメッシュが粗すぎるのだよねー。アメダスって気圧測定してないか公開してないかだし。いや,まじめに地球シミュレータで扱うようなグリッドの細かい気象モデルだと現在のメッシュでは荒すぎて,アメダスデータと対応するような,まともな気象モデルつくれないんじゃないかと。つーか,メッシュ細かくすると現在の気柱モデルに準じた補正大気圧って誤差が大きすぎで,それで公開できないらしい。こうなると卵が先なのか鶏が先なのか・・・
もっとアホな子が考えた場合 (スコア:0)
バージョンアップ、偏微分方程式を解けるようにしてみました(ぼうよみ
もっともっとアホな子が考えた場合 (スコア:0)
謎企業の中の人「なんで?」
あほ「解らないのか?だって、考えてみろよ
”なんだか熱くなってきたわ、アタシぬいじゃおっかなー”とか
”きゃー風さんのエッチ!(スカートがふわり)”
という気候の悪戯、偶然を予測して現場に急行できるンだぜ?」
謎企業の中の人「まじでwww簡単な計算で正確な気圧計算する方法考えてみるわ
ちょっとまちなー」
数ヵ月後・・・
Google天気図Betaサービス開始!