栗田です 昨日は所用があったのと、時差が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ファイル１要素の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が使いこなせればもう少しスマートに書けるのですが、 １週間の勉強ではこれが限界でした。さて、ここで質問があるのですが、このファイルでは、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>