[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[dennou-ruby:003067] Re: Fw: 電脳Rubyで2次元gribファイル読み込み
- To: dennou-ruby@xxxxxxxxxxx
- Subject: [dennou-ruby:003067] Re: Fw: 電脳Rubyで2次元gribファイル読み込み
- From: Naoyuki Kurita <nkurita@xxxxxxxxxxxxx>
- Date: Wed, 4 Mar 2009 19:40:34 -0500
栗田です
堀之内さんの以下のご助言に従って改良をしてみました。
行ったことは、以下にまとめます。
1)grid_copyを使って、軸情報をコピーする
2)insert_axisを使って、時間軸を追加する
3)shapeを使って変数を取り出す
4)Dir.globを使って、存在するファイルだけ
  上書きする
できあがったプログラムは、取り込むファイルが1つ
(例 tmp_2008010100.grib)だけであれば、正しく
NetCDF変換できるのですが、ファイルが複数個
(例 tmp_2008010100.grib, tmp_2008010200.grib)に
なると、3次元データvalueを4次元のArray 
に代入する
37行目で以下のエラーが出ます。
line 37: array[false,jday] = value #jday<-time変数
grib2nc.rb:37: undefined method `[]=' for nil:NilClass (NoMethodError)
	from grib2nc.rb:19:in `each'
	from grib2nc.rb:19
代入するデータや、jdayがおかしいのかと思い出力させてみても
データはそのまま入っています。
$ ruby grib2nc.rb
time -> 0
NArray.float(288,145,23):  <- value
[ [ [ 241.052, 241.052, 241.052, 241.052, 241.052, 241.052,  
241.052, ... ],
    [ 242.271, 242.239, 242.239, 242.208, 242.177, 242.177,  
242.146, ... ],
    [ 244.958, 244.927, 244.896, 244.896, 244.833, 244.802,  
244.771, ... ],
    [ 247.583, 247.614, 247.646, 247.677, 247.677, 247.708,  
247.708, ... ],
    [ 251.677, 251.552, 251.489, 251.521, 251.552, 251.614,  
251.646, ... ],
    [ 261.271, 260.958, 260.739, 260.552, 260.427, 260.271,  
260.114, ... ],
    [ 270.521, 270.271, 270.083, 269.927, 269.771, 269.614,  
269.458, ... ],
    [ 273.489, 273.427, 273.364, 273.239, 273.021, 272.771,  
272.396, ... ],
    [ 273.896, 274.114, 274.208, 274.146, 273.896, 273.427,  
272.802, ... ],
    [ 273.708, 274.208, 274.614, 274.896, 274.958, 274.771,  
274.302, ... ],
 ...
time->1
NArray.float(288,145,23):  <-value
[ [ [ 241.052, 241.052, 241.052, 241.052, 241.052, 241.052,  
241.052, ... ],
    [ 242.271, 242.239, 242.239, 242.208, 242.177, 242.177,  
242.146, ... ],
    [ 244.958, 244.927, 244.896, 244.896, 244.833, 244.802,  
244.771, ... ],
    [ 247.583, 247.614, 247.646, 247.677, 247.677, 247.708,  
247.708, ... ],
    [ 251.677, 251.552, 251.489, 251.521, 251.552, 251.614,  
251.646, ... ],
    [ 261.271, 260.958, 260.739, 260.552, 260.427, 260.271,  
260.114, ... ],
    [ 270.521, 270.271, 270.083, 269.927, 269.771, 269.614,  
269.458, ... ],
    [ 273.489, 273.427, 273.364, 273.239, 273.021, 272.771,  
272.396, ... ],
    [ 273.896, 274.114, 274.208, 274.146, 273.896, 273.427,  
272.802, ... ],
    [ 273.708, 274.208, 274.614, 274.896, 274.958, 274.771,  
274.302, ... ],
 ...
grib2nc.rb:37: undefined method `[]=' for nil:NilClass (NoMethodError)
	from grib2nc.rb:19:in `each'
	from grib2nc.rb:19
時間軸をforではなく、Dir.globを使って、ファイル名 
から時間情報を取り
出したところからこのようなエラーが出ましたので、何かそれに関係して
いると思うのですが、ちょっと検討がつきません。ご助言をお願いでき
ませんでしょうか?
また、プログラム内で、tmp_yyyymmddhh.gribというファイル 
から、
時間情報を取り出しているのですが、
fname=tmp_2008010100.grib
p fname.slice(-15..-6)
 => "2008010100"
p fname.slice(-15..-6).split(/(....)(..)(..)(..)/)
 => ["", "2008", "01", "01", "00"]
と最初に””が入ってしまいます。どうもPerlと勝手が違う 
ようで、
正規表現が間違っているのが大きな原因だと思いますが、、、、
お忙しいところお手数おかけして申し訳ございませんが
お助け願います。余計な事に手を出して失敗するなんて
情けない。。。。
require 'numru/gphys'
require 'narray_miss'
require 'date'
include NumRu
rmiss   = -999
iflag   = 1
yr      = 2008
var     = 'tmp'
day1 = Date.new(yr,1,1)
nd   = ( Date.new(yr+1,1,1)-day1 ).to_i
t    = VArray.new( NArray.sfloat(nd).indgen,
             {"units"=>"days since"+yr.to_s+"-1-1 0" },
	         "time" )
time = Axis.new.set_pos(t)        
Dir.glob(var+"_"+yr.to_s+"??????.grib").each{ |fname|
    (dummy,y,m,d,h) = fname.slice(-15..-6).split(/(....)(..)(..)(..)/)
p    jday = ( Date.new(y.to_i,m.to_i,d.to_i)-day1 ).to_i
    varname = GPhys::IO.var_names( fname )[0]
    gphys   = GPhys::IO.open( fname,varname )
p    value   = gphys.val 
    if iflag == 1 then 
         axis  = gphys.grid_copy
          long = gphys.data.get_att('long_name')
          unit = gphys.data.get_att('units')
          name = gphys.name
	  (x, y, z) = gphys.shape
	  array = NArray.sfloat(x,y,z,nd).fill!(rmiss)	
	  iflag = 0
   end
   array[false,jday] = value
}
data = VArray.new( array, {"long_name"=>long, "units"=>unit}, name )
p gphys = GPhys.new( axis.insert_axis!(-1,time), data )
 
file = NetCDF.create("test.nc")
GPhys::NetCDF_IO.write(file,gphys)
file.close
On 2009/03/04, at 9:26, Takeshi Horinouchi wrote:
栗田さま
動いてよかったです.プログラムはだいぶ短くなりましたね.
ここから先は趣味の問題みたいなもんですが,
さらに短く,かつ見通しもよりよく(これが大事)できると思います.
なにせ,形式変換プログラムが簡単にかけるのは GPhys の
売りのひとつなので,ついついベストを追求したくなります.
データが手元になく確認できないので方針だけになりますが.
GRIBファイルをあらわす gphys の格子コピーを作成して
insert_axis 等を使います.すると,lon = gphys.axis(0)
等は不要になります.また,GPhysやNArray等の shape
メソッドを活用するともっと簡単になるところもあります.
あと,年間のループをまわした上でファイルがあるときだけ
対応する代わりに Dir.glob で存在するふぁいるだけ
拾ってループをまわすこともできます.
# なにせ趣味の領域なので,栗田さんにやってみてくださ 
いという
  つもりはありませんが.
堀之内
ご助言に従って、複数の変数が日付毎に格納されている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でダウンロードできますので、興味が
あればのぞいてみ
てください。また、制作者がすぐ近くにおりますので、質問があれ
ば、こちら
で情報を仕入れることもできます。