[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[dennou-ruby:003053] Fw: 電脳Rubyで2次元gribファイル読み込み
- To: dennou-ruby@xxxxxxxxxxx
- Subject: [dennou-ruby:003053] Fw: 電脳Rubyで2次元gribファイル読み込み
- From: Takeshi Horinouchi <horinout@xxxxxxxxxxxxxxxxx>
- Date: Tue, 03 Mar 2009 23:27:21 +0900
堀之内です．
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 --------------------
堀之内 武
北海道大学 地球環境科学研究院 地球圏科学部門
require 'numru/gphys'
require 'narray_miss'
include NumRu
rmiss=-999
# -------------------------------------------------------------------
def yymmdd_to_ddd( y, m, d )
dmon=[0,31,28,31,30,31,30,31,31,30,31,30,31]
dmon[2] = 29 if y%4==0 && (y%100 != 0 || y%400 == 0)
ddd = d
while m > 0
ddd = ddd + dmon[m-1]
m = m -1
end
return ddd
end
def yyddd_to_mmdd ( y,d )
dmon=[0,31,28,31,30,31,30,31,31,30,31,30,31]
dmon[2] = 29 if y%4==0 && (y%100 != 0 || y%400 == 0)
m = 1
while d > dmon[m] && m < 12
d = d - dmon[m]
m = m + 1
end
return m, d
end
# -------------------------------------------------------------------
yr = 2004
fname = 'tmp'
varname = 'TMP'
iflag = 1
jend = yymmdd_to_ddd( yr, 12, 31 )
for i in 1...jend do
mo,dy = yyddd_to_mmdd( yr, i )
yy = sprintf("%2.2d",yr)
mm = sprintf("%2.2d",mo)
dd = sprintf("%2.2d",dy)
file = fname + yy + mm + dd + "00.grib"
if File.exist?( file ) then
gphys = GPhys::IO.open( file,varname )
value = gphys.val
if iflag == 1 then
x = gphys.coord(0)
y = gphys.coord(1)
z = gphys.coord(2)
t = VArray.new( NArray.sfloat(jend-1).indgen,
{"units"=>"days since"+yr.to_s+"-1-1 0" },
"time" )
lon = Axis.new.set_pos(x)
lat = Axis.new.set_pos(y)
lev = Axis.new.set_pos(z)
time= Axis.new.set_pos(t)
long = gphys.data.get_att('long_name')
unit = gphys.data.get_att('units')
name = gphys.name
array = NArray.sfloat(x.length,y.length,z.length,jend-1).fill!(rmiss)
iflag = 0
end
for k in 0...z.length do
for m in 0...y.length do
for l in 0...x.length do
array[l,m,k,i-1] = value[l,m,k]
end
end
end
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