[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[dennou-ruby:003061] Re: Fw: 電脳Rubyで2次元gribファイル読み込み



栗田です

昨日は所用があったのと、時差が14時間ありますので(来週
からは夏時間なので13時間になります)、レスが遅くなり申
し訳ございません。また、たくさんのフォローありがとうご
ざいました。

ご助言に従って、複数の変数が日付毎に格納されているJRA25
データを、モデル出力と同じように、1変数ずつの4次元 NetCDF データに変換するプログラムをリバイスしてみました (grib2nc.rb)。

使用しているデータは、JRA25データのanal_p/anl_p25 と同じ
要素が入った、水平1.25度格子データです。

http://jra.kishou.go.jp/JRA-25/elements.html

これは、6時間値毎に1ファイルになっておりますので、各 ファイル から指定する要素を抜き出し、1ファイル1要素のgrib ファイルを 1年分つくり、それをRubyプログラムで4次元データに するという流れ
になっています(添付したgrib2nc.sh参照)。

Rubyプログラムでは、例えば
TMP2008010100.grib
TMP2008010200.grib
:
:

といったように複数個に分けられているGrib形式の気温デー タを読み込み これを4次元Arrayにした後に、NetCDFで書き出 すということをしてい
ます。また、たまに欠損時間がありますので、それを考慮して、最初に、
365日分の全データに欠損値を入れておき、

Line 38:  array = NArray.sfloat(x,y,z,nd).fill!(rmiss)

そこにデータがある日にちだけデータを上書きして行くことを行ってい
ます。ちょっとTMP2008010100.grib, TMP2008010200.gribと いうファ
イルを作って

ruby grib2nc.rb

としてみたところ、出力データもちゃんと出力できました。

<GPhys grid=<4D grid <axis pos=<'lon' in 'TMP'  288>>
	<axis pos=<'lat' in 'TMP'  145>>
	<axis pos=<'level' in 'TMP'  23>>
	<axis pos=<'time' sfloat[1] val=[0.0,]>>>
data=<'TMP' sfloat[288, 145, 23, 1] val = [241.051971435547,241.051971435547,241.051971435547,241.051971435547 ,...]>>

HGTを登録できれば、Rubyを使ってJRAデータもGrib->NetCDF が便利に できそうです。以前は、Gradsでも、lats4dコマンドを 使って、 Grib->NetCDF変換が便利にできたのですが、最新のGrads2.0 ではサポート されなくなってしまいまいましたので、電脳rubyで簡単にで きると非常に嬉
しいです。

また、NetCDFデータの可視化ですが、現在いるゴダード宇宙 科学研究所 では、panoplyというJavaソフトがあり、コンパイルせ ずに、ただファイル を置いた場所にパスを通すだけで全球データのクイックルックがで きます。

http://www.giss.nasa.gov/tools/panoply/

上記のサイトからfreeでダウンロードできますので、興味が あればのぞいてみ てください。また、制作者がすぐ近くにおりますので、質問があれ ば、こちら
で情報を仕入れることもできます。

require 'numru/gphys'
require 'narray_miss'
require 'date'
include NumRu

rmiss   = -999
iflag   = 1

yr      = 2008
fname   = 'tmp'

day1 = Date.new(yr,1,1)
#nd   = ( Date.new(yr+1,1,1)-day1 ).to_i
nd   = ( Date.new(yr,1,2)-day1 ).to_i

for i in 0..nd do
    date = day1 + i
    file = date.strftime("#{fname}%Y%m%d00.grib")
    if  File.exist?( file ) then
       varname = GPhys::IO.var_names( file )[0]
       gphys = GPhys::IO.open( file,varname )
       value = gphys.val 

       if iflag == 1 then 
          lon = gphys.axis(0)
          lat = gphys.axis(1)
          lev = gphys.axis(2)
          t = VArray.new( NArray.sfloat(nd).indgen,
                    {"units"=>"days since"+yr.to_s+"-1-1 0" },
		    "time" )
          time= Axis.new.set_pos(t)        

          long = gphys.data.get_att('long_name')
          unit = gphys.data.get_att('units')
          name = gphys.name
       	  x    = gphys.coord(0).length
       	  y    = gphys.coord(1).length
       	  z    = gphys.coord(2).length
	  array = NArray.sfloat(x,y,z,nd).fill!(rmiss)	
	  iflag = 0
       end

       array[false,i-1] = value

    end
end

data = VArray.new( array, {"long_name"=>long, "units"=>unit}, name )
gphys = GPhys.new( Grid.new(lon,lat,lev,time), data )
 
file = NetCDF.create("test.nc")
GPhys::NetCDF_IO.write(file,gphys)
file.close




Attachment: grib2nc.sh
Description: Binary data



On 2009/03/03, at 9:27, Takeshi Horinouchi wrote:

堀之内です.

JAMSTECの栗田さん(現在in NY)からこのメールの最後につけた
メールを貰いました.うっかり返信がおそくなってしまい,さき
ほどごめんなさいの連絡をしたところ,

また、試行錯誤した結果、gribからNetCDF変換は、 変数登録
されているものしかできないという気がしていますが、ソース
の解読ができていないので、まだ確信はありません。3次元デ
ータでも、何故かジオポテンシャル高度(HGT)を変換しようと
すると以下と同様のエラーが出ます。おそらく、NCEP等の
変数名にあわせるのかと推察しています。

というフォローも貰いました.(なお,現在は電脳ruby ML
に再登録ずみなので,このメールは届くと思います.)

引用の前に,とりあえずの返信を書きます.
確かに,grib の場合,登録された変数しか扱えません.
GRIBの場合,ご存知のように物理量の種類は番号で
指定されますが,GRIB標準以外に独自拡張もあったり
して,Debianパッケージの標準のインストール場所でいえば
/usr/lib/ruby/1.8/numru/gphys/ に grib_params.rb
というファイルがあり,その中に対応が定義されてます.
GRIBライブラリは西澤さんが作ってくれましたが,
番号に対する名前は,wgrib からとったのだと思います.
ruby上で,みたいファイルにどんな名前の変数が入っている
とみなされるかをみるには,

$ irb
irb(main):001:0> require "numru/gphys"
=> true
irb(main):002:0> include NumRu
=> Object
irb(main):003:0> GPhys::IO.var_names "T.jan.grib"
=> ["TMP"]

などとします.ここまではきっと解読済みですね.

さて,問題のエラー

irb > gphys = GPhys::IO.open( PRMSL2004010100.grib,"PRMSL" )

/usr/lib/ruby/1.8/numru/gphys/grid.rb:213:in `initialize': each
argument must be an Axis (ArgumentError)
	from /usr/lib/ruby/1.8/numru/gphys/grid.rb:209:in `each'
	from /usr/lib/ruby/1.8/numru/gphys/grid.rb:209:in `initialize'
	from prs_grib2nc.rb:78:in `new'
	from prs_grib2nc.rb:78

ですが,座標軸ででてるので,上に書いたような問題では
ありませんね.実際,MSL という変数名は登録済みですし.

もしかしたら,Ruby の GRIB ライブラリで扱えない タイプの
座標系だったりしないでしょうか.確か緯度経度座標にしか
対応してなかったと思いますので,地図投影座標や
international exchange grid (極近くほど格子点が
間引かれる) だったり,すると扱えないことになります.
その場合,残念ながら,根本的な対応をしないと
扱えないということになりますが,さて,実際はどうでしょうか.

後半のほうは,プログラムをまだ読んでないので,またあらためて
ということで.

===========================================================
堀之内様

さて、本日は、電脳Rubyに関してお伺いしたいことがあり
メールいたしました。本来なら、MLに投稿すべきなのです
が、米国から投稿しているためか、何故か登録を受け付けて
もらえませんでしたので、直接メールさせて頂きます。

【動機】
これまで観測データ(text)とGCMを使って仕事をしており
客観解析データもgtool形式に直して仕事しておりました。
しかし、最近、諸般の事情でNetCDFのデータを大量に扱
うことになり、HDDを節約するためにも、これまでの
gtool形式ではなく、NetCDFに統一しようと準備を進めて
います。また、NetCDFを使った解析を行う場合、Perlだと
ファイル読み書きに時間がかかり、ストレスフルなので、
NArrayやGphysという便利なツールが満載のRubyに乗り換
えようかと考えております。沼さんが他界される直前の2001
年の春に話をしたときも、これからは "Rubyだよ" と
言われて
おり、あれから8年も経過してしまいましたが、時間をみつけ
てここ1週間前から勉強を始めました。

【問題発覚】
まずは、grib形式のJRA25データをGphysを使っ
て、Netcdf
変換しようと思い、試しに、地表面気圧データ
( PRMSL2004010100.grib)を読み込んでみました。この中身は、

$ wgrib -s PRMSL2004010100.grib
=>  1:0:d=04010100:PRMSL:MSL:anl:NAve=0

となります。これをコマンドラインで読み込むと、データを
OPENするところでエラーストップします。

irb > require 'numru/gphys'
irb > include NumRu
irb > gphys = GPhys::IO.open( PRMSL2004010100.grib,"PRMSL" )

/usr/lib/ruby/1.8/numru/gphys/grid.rb:213:in `initialize': each
argument must be an Axis (ArgumentError)
	from /usr/lib/ruby/1.8/numru/gphys/grid.rb:209:in `each'
	from /usr/lib/ruby/1.8/numru/gphys/grid.rb:209:in `initialize'
	from prs_grib2nc.rb:78:in `new'
	from prs_grib2nc.rb:78

試しに、TMP2004010100.gribという3次元データに変えて
みたところ、こちらは問題なく読めます。どうもGrib.open
のところで読み込みに失敗しているようですが、Ruby素人
なので、原因究明ができませんでした。他の2次元データも
同様で、bin形式にして、ctlファイルを作ってトラ イしまし
たが、結果は同様でした。これをgrads等を用いてNetCDF
に変換すればファイルは問題なくOPENできますので、どうも
gribやgrads形式では、2次元データサポートしていな
いように
見受けられますが、何か裏技があるのでしょうか?Gphysの
バージョンは、ちょっとわからないのですが、今月初めに
Debian etch用のパッケージをapt-getしてインストールし ました。


【素人プログラム】

上記にも関係しますが、JRA25を一要素ずつ時系列毎に並べた
COARDS形式の4次元NetCDFデータを作成するためのプロ
グラム
を作ってみました(添付ファイル)。これは、wgribで各要素の
データを1日毎のファイルに分割し、JRA25のデータが
ない日は
欠測値(-999)を入れて時系列データを作るようにしてあ ります。
Rubyが使いこなせればもう少しスマートに書けるのですが、
1週間の勉強ではこれが限界でした。

さて、ここで質問があるのですが、このファイルでは、66 行目から
72行目のところで、3次元データ(value)を読み込ん
で、時間軸を加え
た4次元データ(array)を出力しているのですが、どう
やってもRuby
らしく書くことができませんでした。Gphysの機能を使って、
もっと
スマートに書けると思いますが、堀之内さんならどのようにされま
すか?向学の為にご助言頂ければ幸いです。

お忙しい中、長文で失礼いたしました。また、今後ともどこかで
お会いした際には、どうかよろしくお願いいたします。

#本来なら書物をもっと読んでメールすべきですが、
#現在米国に滞在中でして、日本語の本が入手困難な
#状況です。なんとかネットを駆使して勉強しており
#ますが、至らない点が多々あると思います。その点
#はどうかご容赦ください。

*************************************************
栗田直幸 KURITA Naoyuki
独立行政法人 海洋研究開発機構 (JAMSTEC)
地球環境観測研究センター (IORGC)
水循環観測研究プログラム
〒237-0061 神奈川県横須賀市夏島町2-15
電話:046-867-9822
FAX:046-867-9255
E-mail:nkurita@xxxxxxxxxxxxx
http://www.jamstec.go.jp
*************************************************

--------------------- Original Message Ends --------------------

堀之内 武
北海道大学 地球環境科学研究院 地球圏科学部門
<grib2nc.rb>