[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でダウンロードできますので、興味が
あればのぞいてみ
てください。また、制作者がすぐ近くにおりますので、質問があれ
ば、こちら
で情報を仕入れることもできます。