気象データを解く〜データをパソコンで解析する その2

  • 2012.08.14 Tuesday
  • 23:01
気象データを解析しようと思うが、
どこから入手すればいいのか、どうやって解析すればいいのかについて覚え書き

できるだけ多くの気象に興味のある人(特にこれから研究しようとしている学生さん)に、
気象の知識だけでなく、コンピュータで気象解析ができるようになってくれればと思っています。


その1では、ファイルを開くところや複数オクテットの数値演算を行ってきた。

その2では符号付きやデータ切りだしを考えていく。

なお、ソースはあくまで参考程度に必要最低限しか書いていないので、
バグを含んでいる可能性があるので、丸ごとコピーせずに考えながら処理していただきたい。



続きを読む >>

関連投稿
気象データを解く〜データを入手する
気象データを解く〜データを開く、前に用語とか整理
気象データを解く〜データを開く まずは基本的な第0節から
気象データを解く〜データを開く 第1節を読む
気象データを解く〜データを開く 第3節を読む
気象データを解く〜データを開く 第4節を読む
気象データを解く〜データを開く 第5節を読む
気象データを解く〜データを開く 第7節を読む
気象データを解く〜データを開く 続くデータを見ていく
気象データを解く〜データを開く surfファイルを読む

気象データを解く〜データをパソコンで解析する その1
気象データを解く〜データをパソコンで解析する その3
気象データを解く〜データをパソコンで解析する その4
気象データを解く〜データをパソコンで解析する その5


単位変換や気象計算覚え書き
単位変換覚え書き〜風向風速

 

符号が付いている場合(第4節24オクテット目や第5節など)は先頭が0か1かで判断して、
残りは符号なしの場合と同じように解くことができる。

符号か符号でないかの判別は理論演算ANDを使うとわかる。
(理論演算の詳しい説明は本などを読んでください)

つまり、あるオクテットの先頭を80(128)でANDした結果が80なら負数、0なら正数である。
ソースを書くなら、
 
If (rByte(0I) And &H80) = &H80 Then
  '負数である
End If

ということになる。
これを踏まえてByteToInt関数を改良してみると
 
Private Function ByteToSignedInt(ByRef Source() As Byte, ByVal Start As Int32, ByVal Length As Int32) As Int64
'Int32でもよいが、オーバーフローを考えるとInt64の方がいい
‌ 
        Dim rVal As Int64 = 0L
        Dim intB As Int32
        Dim intE As Int32 = Start + Length - 1I
        Dim bSng As Boolean = True
‌ 
        '先頭オクテットの確認
        If (Source(0I) And &H80) = &H80 Then '符号があるか
            bSng = True '符号有り
            Source(0I) = Source(0I) And &H7F '符号を取る
        End If
‌ 
        For intB = Start To intE
‌ 
            rVal *= CLng(&H100)
            rVal += CLng(Source(intB))
‌ 
        Next intB
‌ 
        If bSng = True Then '符号有り
            rVal *= -1'符号を付ける
        End If
‌ 
        Return rVal
‌ 
End Function

こんな感じだろうか。


これを応用してIEEE754を解くことができる。

IEEE754についてはGRIBの仕様書にも書いてあるので、詳しくは触れない(ソースコードを載せることもない)が、
ここでは値の取り出し方法を考えてみることにする。

以前解いた42 51 8B DEを使うと
(rByte(x) rByte(x+1) rByte(x+2) rByte(x+3)にデータがあるとする)

まず先頭の42を80でANDして符号か符号でないかを判断する
 
Dim bSng As Boolean = False
‌ 
If (rByte(x) And &H80) = &H80 Then '符号が付いているか
   bSng = True '符号あり
End If


次のE部が面倒である。なぜかというと2進数に直すと

0100 0010 0101 0001 1000 1011 1101 1110

緑色の箇所が先頭と次のオクテットの2つにまたがっているからである。

これをプログラムで抜き出すためにはどうしたらよいか。

まず、それぞれ必要な箇所を抜き出してみる。
つまり、先頭は7FでAND、次は80でANDするとそれぞれ42と0が返る。

続いてE部は8ビットなので、1000 0100(132)にしたいが、どうすればこうなるのか。

単純に42 + 0をしても42のままである。

どうにかして先頭の42を1ビット分桁を上げられればよいのだが・・・。

少し考えてみよう。たとえば、10進法で、1を1桁上げたいと考えるなら10をかければ10となる。
これと同様に2進数でも1桁上げたい場合は、何かをかけてあげればよいことになる。
では何をかけるか?
1を1桁あげると10、つまり10進でいうところの2になる。つまり、進数の数2をかければ1桁をあげられるわけだ。

早速計算してみると、

42(66) x 2 + 0 = 84(132)


となり84という値が得られる。2進で書くと 1000 0100 一つずつ桁があがっていることがよく分かる。

さらっと流したが、第0節を解読するときにF(16)や、前回のその1で100(256)をかけ算していたのは、16進や256進の桁上げを意味している。

蛇足してしまったが、ソースで書くなら、
 
Dim ex1 As Int32 = CInt(rByte(x) And &H7F)
Dim ex2 As Int32 = CInt(rByte(x+1) And &H80)
‌ 
Dim exA As Double = CDbl(ex1 * &H2 + ex2)

という感じであろう。

※コメント頂いたように、例の場合、たまたまx+1の配列をANDした結果0が得られたが、
実際ANDした結果、H80が返ってくる場合もある。

ソースを汎用的に書き直すのであれば、
 
ex2 = CInt(rByte(x+1) And &H80) / &H80 '8ビット右へシフト

とした方が良い。

また、x+1から取り出す値は1ビット分なので、
計算量を減らすという意味で(減っているか分からないが)、符号付きと同じ考え方で、
 
If (rByte(x+1) And &H80) = &H80 Then
 exA = CDbl(ex1 * &H2 + 1I)
Else
 exA = CDbl(ex1 * &H2)
End If

という感じでも良いかもしれない。


最後はM部であるが、これは簡単。

M部の先頭のみ7FでANDして51を取り出して、残りはそのまま桁をそろえて加算すれば、
先の投稿の通りの数値が得られる。

51 x 100 x 100 + 8B x 100 + DE = 518BDE(5344222)

ソースを書くと、
 
Dim ma1 As Int32 = CInt(rByte(x+1) And &H7F)
Dim maA As Double = CDbl(ma1 * &H100 * &H100 + rByte(x+2) * &H100 + rByte(x+3))

でどうだろうか。

あとはIEEE754の仕様に沿って計算すれば答えが得られる。

※上のソースコードであえてDouble型に型変更しているのは、IEEEの計算過程でDouble型での演算が発生するためである。
(ビルド時に勝手にコンパイラが型変換してくれるので、わざわざこういう処理を書かなくてもよいとは思うが)

ちょっと端折ったかなと感じたので、ソースを加筆しておく。
 
Dim rVal As Double
‌ 
If exA >= 1.0R AndAlso exA <= 254.0Then 'E部 = 1〜254
‌ 
  rVal = (1.0+ maA * 2.0R ^ -23.0R) * 2.0R ^ (exA - 127.0R)
  '上式演算でべき乗計算がDouble型なので、maAとexAをDouble型にした。
‌ 
End If
‌ 
'符号は仕様書のように上式で-1sをしてもよい
If bSng = True Then '符号有り
  rVal *= -1.0'符号を付ける
End If

他も同様に場合分けして演算すればよい。
なお、IEEE754は単精度なので、Single型で計算した方がよいかもしれないが、べき乗計算においてどうしてもDouble型で計算されてしまうので、Double型で演算を行っている。



さて、残るはデータの解読である。
幸いデータの取り出し長は12ビット毎なので、簡単である。

第7節でも使った下記で取り出しを考えてい。

7D D7 D4 7D B7 BF

1オクテットは8ビットである。12ビットにするためには、2オクテット取り出す必要がある。
なので、組み合わせとして、8ビット+4ビット、4ビット+8ビットという形で順に取得していくことになる。

,泙此2オクテット 7D D7をrByteから切り出して、
7Dはそのまま、D7はF0(240)でANDしてD0を取り出す。

続いて0111 1101 1101の12ビットにするためには、
7Dは4桁上げ、D0は4桁下げる必要がある。そして結果を足し算すればよい。
4桁上げ下げすると言うことは、全体として10(16)倍してあげればよいので(※)、

7D x 10(16) + D0 ÷ 10(16) = 7DD(2013)

ということになる。
確認のため7DDを2進にすると0111 1101 1101である。

※2進数で考えるなら、1を10000にしてあげればよい。10000は10進で16である。


ソースにすると
 
Dim pack1 As Int32 = rByte(x)
Dim pack2 As Int32 = rByte(x+1) And &HF0
‌ 
Dim packA As Int32 = pack1 * &H10 + pack2 / &H10

という感じで取り出せる。


同様に次のデータを取り出す。
次の2オクテットはD7 D4に含まれるので、rByteから切り出してあげる。
必要な部分は、D7はF(15)とANDして7を取り出し、D4はそのまま取り出す。

今度は7を8桁上げ、D4はそのままでよいので、
(8桁は16進なら100倍)

7 x 100(256) + D4 = 7D4(2004)

つまり、0111 1101 0100が得られる。


ソースを書くなら
 
Dim pack1 As Int32 = rByte(x+1) And &HF
Dim pack2 As Int32 = rByte(x+2)
‌ 
Dim packA As Int32 = pack1 * &H100 + pack2

という感じになる。


これを´↓´◆ΑΑΔ鳩り返していけばよい。
もし12ビット以外の場合でもこの要領で取り出していけば取り出せるはずである。


こちらも端折ったので、ソースを加筆すると、
最終的に取り出す公式(解法式)が Y x 10D = R + (X1 + X2) x 2E なので、
 
Dim sclB As Double = 2.0R ^ -5.0'5節の16〜17オクテット
Dim sclD As Double = 10.0R ^ 0.0'5節の18〜19オクテット
‌ 
Dim vOrg As Double = (rVal + CDbl(packA) * sclB) / sclD

と解放すればよい。vOrgが最終的に求める値である。


ざっくりと考え方だけ書いたが、ようするにByteToInt関数を応用して改良していくだけなので、がんばってソースコードを書いてみて欲しい。

実際に書くときは、いろいろなテクニックを駆使した方がよりすっきしとしたソースになると思う。





あくまでサンプルソースですので、異常処理などは一切書いていません。
その辺は自分で悩んで処理を追加してください。
また、あらかじめ参照やインポートしているライブラリや端折って書いていることがあります。
その辺は暗黙で書いているので、サンプルをそのままコピペしても異常や注意情報が上がることがあります。
詳しくはネットや本を調べれば分かることなので、その辺も悩んで見てください。



関連投稿
気象データを解く〜データを入手する
気象データを解く〜データを開く、前に用語とか整理
気象データを解く〜データを開く まずは基本的な第0節から
気象データを解く〜データを開く 第1節を読む
気象データを解く〜データを開く 第3節を読む
気象データを解く〜データを開く 第4節を読む
気象データを解く〜データを開く 第5節を読む
気象データを解く〜データを開く 第7節を読む
気象データを解く〜データを開く 続くデータを見ていく
気象データを解く〜データを開く surfファイルを読む

気象データを解く〜データをパソコンで解析する その1
気象データを解く〜データをパソコンで解析する その3
気象データを解く〜データをパソコンで解析する その4
気象データを解く〜データをパソコンで解析する その5


単位変換や気象計算覚え書き
単位変換覚え書き〜風向風速
 
コメント
以前、風向風速のところでコメントしたものです。
MSMのことを調べていたら、またロケッこさんのページに
たどり着きました。気象業務支援センターの技術資料より
圧倒的にわかりやすく、とても助かります!
VBAで処理しようとしてまして、ビット処理でまだ頭を悩ませてますが、ロケッこさんのページを見て、なんとかいけそうな気がしてきました。
  • またかず
  • 2013/05/21 10:50 PM
またかずさま
お久しぶりです。コメントいただきありがとうございます。

確かにセンターの資料はじっくり読まないと肝心なところが分かりづらいかもしれません。

VBAからずいぶん離れてしまったので、アドバイスはできませんが(VB.netなら現役なのですが)、
解読がんばってください。
  • rockecco
  • 2013/05/21 11:02 PM
ロケッこさま

E部分の処理ですが、

Dim ex1 As Int32 = CInt(rByte(x) And &H7F)
Dim ex2 As Int32 = CInt(rByte(x+1) And &H80)
Dim exA As Double = CDbl(ex1 * &H2 + ex2)

となっておりますが、rByte(x+1)の最上位ビットが0の場合はよいですが、1の場合はAND&H80をした後に&H80で割ってからex1*&H2に足すのではないでしょうか?
  • またかず
  • 2013/05/23 10:05 PM
またかずさま

ご指摘ありがとうございます。

確かに、上の例の場合は、このままで通ってしまいますが、
実際にはx+1は8ビット右へシフトした形にしなければならないので、
x+1をH80でANDした結果をさらにH80で割ってあげる(シフト演算になれていれば8ビット右へシフト)必要があります。

または、符号ビットと同じように、ANDした結果、H0かH80が出てくるかで、加算しないか、1を加算するという方法でもよいかもしれません。


後ほど訂正させて頂きます。
  • rockecco
  • 2013/05/23 10:44 PM
ロケッこさま

ロケッこさんの別案のやりかたもありますね。
今回の件は、たまたま計算がオーバーフローしたので、気づきました。
ほんと、ロケッこさんの記事のおかげで助かっております。
それにしても、grib2は私には扱いにくいです。
NetCDF形式の方が、楽なんですが・・・
  • またかず
  • 2013/05/23 11:16 PM
またかずさま

誤りを見つけていただき、こちらこそありがとうございます。

私もはじめはGRIB2は苦戦しました。
とくに、符号付きか、なしかなど細かい仕様は技術資料では不足がちなので、
結局WMOの仕様を読まないと分からないこともありますし。

そういった自分の経験を残しておくことで、
少しでも生データを活用できる足場になってもらえればいいかなと思っております。
  • rockecco
  • 2013/05/24 7:18 AM
コメントする








    
この記事のトラックバックURL
トラックバック

calendar

S M T W T F S
      1
2345678
9101112131415
16171819202122
23242526272829
30      
<< September 2018 >>

search this site.

よく使う、検索される投稿

categories

アマゾン

楽天

selected entries

archives

recent comment

recent trackback

profile


※当ブログはリンクフリーですが、 取材や雑誌等で掲載される場合は、事前にお知らせください

others

mobile

qrcode

powered

無料ブログ作成サービス JUGEM