#!/usr/bin/env ruby # # 表題: 熱帯降水活動コンポジット解析 # # 履歴: 2003/11/15 やまだ由 # 履歴: 2003/11/18 やまだ由 # # ------------------------- require "numru/ggraph" require "colorbar.rb" include NumRu include NMath # ------------------------- # 実験名 # $cum_param = "adj" $cum_param = "kuo" $exp = "con" # $exp = "mradlAa" # $exp = "mradlAb" # $exp = "mradlAc" # $exp = "mradlAd" # $exp = "fqradl-2K" # $exp = "fqradl-1K" # $exp = "fqradl-05K" # 閾値 threshold = 3000.0 # NetCDF ファイルの場所 $filepath = "/work11/aqua3d/yukiko/mradl/composite/data/" $filepath_out = "/home/yukiko/work/aqua-mradl-ronbun/composite_nc" END{ # mradl_optdcl # ディレクトリ名 $dir = "#{$cum_param}700#{$exp}" # 出力ファイル名 $outfile = "#{$dir}_composite_tuv_y32.nc" composite($dir,threshold) } # ------------------------- def composite(dir,threshold) # コンポジットを評価する rain 初期化 rain_data = GPhys::NetCDF_IO.open("#{$filepath}/#{$dir}/rain_sum_seq.nc",'rain') rain_data_cal = NArray.sfloat(128+6).fill!(0.0) # コンポジットをとる変数の指定 item = ["u", "sigdot", "t" ] # 変数初期化 data = {} ; $comp = {} ; tmp = {} comp_num = 0 ; count = 0 comp_point = NArray.sfloat(128,400).fill!(0.0) item.each{ |i| # 変数の読み込み data[i] = GPhys::NetCDF_IO.open("#{$filepath}/#{$dir}/#{i}_out_seq.nc",i) data[i] = data[i].val[true,-16..-1,true,0] # コンポジット値初期化 $comp[i] = NArray.sfloat(128,16).fill!(0.0) } 400.times{ |time| # コンポジット格子値計算用 rain データ作成 rain_data_cal[3..-4] = rain_data.val[true,time,0] # 重ね合わせ (rain_data_cal.size-6).times{ |comp_num| # 周囲よりも値が大きければその点を中心経度に移動 if rain_data_cal[comp_num+3] > threshold then num = comp_num + 3 if rain_data_cal[num] > rain_data_cal[num - 3] && rain_data_cal[num] > rain_data_cal[num - 2] && rain_data_cal[num] > rain_data_cal[num - 1] && rain_data_cal[num] > rain_data_cal[num + 1] && rain_data_cal[num] > rain_data_cal[num + 2] && rain_data_cal[num] > rain_data_cal[num + 3] then comp_point[comp_num,time] = 1.0 print '(', time, ',', comp_num, ')', "\n" puts rain_data_cal[num] # n = 63 - comp_num n = 64 - comp_num item.each{ |i| tmp[i] = NArray.sfloat(128,16).fill!(0.0) if n == 0 then tmp[i] = data[i][true,true,time] elsif n > 0 then tmp[i][0..(n-1),true] = data[i][(128-n)..127,true,time] tmp[i][n..127,true] = data[i][0..(127-n),true,time] else tmp[i][(128+n)..127,true] = data[i][0..(-1-n),true,time] tmp[i][0..(127+n),true] = data[i][(-n)..127,true,time] end $comp[i] = $comp[i] + tmp[i] } count = count + 1 end end } } # 平均 ------------------------------------------ puts count item.each{ |i| $comp[i] = $comp[i] / count # t, u, sigdot は平均からの偏差 composite_sub(dir, i) if i == "t" || i == "u" || i == "sigdot" } # Gphys オブジェクト生成 & nc ファイル生成 ---- # 出力ファイル ncfile = NetCDF.create("#{$filepath_out}/#{$outfile}") # 軸の値を保存 data_u = GPhys::NetCDF_IO.open("#{$filepath}/#{$dir}/u_out_seq.nc",'u') grid_lon = data_u.axis(0) grid_sigmap = data_u.axis(1) grid_time = data_u.axis(2) # composite point grid = Grid.new(grid_lon, grid_time) comp_point = VArray.new(comp_point).rename("comp_point") comp_point = GPhys.new(grid, comp_point) GPhys::NetCDF_IO.write( ncfile, comp_point ) # composite した変数 grid = Grid.new(grid_lon, grid_sigmap) item.each{ |i| if i == "t" || i == "u" || i == "sigdot" then $comp[i] = VArray.new($comp[i]).rename("#{i}_anm") else $comp[i] = VArray.new($comp[i]).rename(i) end $comp[i] = GPhys.new(grid, $comp[i]) GPhys::NetCDF_IO.write( ncfile, $comp[i] ) } ncfile.close end # ------------------------- def composite_sub(dir, i) # 変数初期化 data_mean = {} # 平均からの偏差 # data_mean[i] = # GPhys::NetCDF_IO.open("#{$filepath}/#{$dir}/#{i}_y32_sum.nc",i) # data_mean[i] = data_mean[i].val[true,-16..-1,0].mean(0) data_mean[i] = $comp[i].mean(0) 128.times{ |num| $comp[i][num,true] = $comp[i][num,true] - data_mean[i] } end # dcl --------------------------------------------- # いちまいめ ----- def dcl_s(num,dir) #ncfile 名 ncfile = $outfile # 変数の読み込み comp_point = GPhys::NetCDF_IO.open("#{$filepath}composite_nc/#{ncfile}", 'comp_point') # 設定 ---------- DCL.gropn(num) DCL::sgpset('lfull',true) # 全画面表示 DCL.sgpset('lcorner',false) # コーナーマークを書かない DCL.uzfact(1.0) if num == 2 DCL.uzfact(0.5) if num == 1 DCL::slmgn(0.0, 0.0, 0.0, 0.15) # 描画マージン指定 GGraph.set_fig('viewport'=>[0.2,0.87,0.12,0.6]) # タイトル #DCL::slsttl('#date #time YAMADA.YU', 'T', 1.0, 1.0, 0.01, 1) DCL::slsttl("#{dir} y=32 composite sampling point", 'T', -0.3, -0.8, 0.02, 2) # composite 図に用いた点 GGraph.contour( comp_point ) unless $graph == 1 then fig = "figs/sampling-#{$exp}-#{$cum_param}-comp" mradl_dcldump(fig) if num == 1 mradl_dclps(fig) if num == 2 end # お片付け ---- DCL.grcls end # にまいめ ------- def dcl_t(num,dir) # 絵を描く変数の指定 item = ["u", "sigdot", "t"] $comp = {} #ncfile 名 ncfile = $outfile item.each{ |i| # 変数の読み込み $comp[i] = GPhys::NetCDF_IO.open("#{$filepath}composite_nc/#{ncfile}", "#{i}_anm") } # 設定 ---------- DCL.gropn(num) DCL::sgpset('lfull',true) # 全画面表示 DCL.sgpset('lcorner',false) # コーナーマークを書かない DCL.uzfact(1.0) if num == 2 DCL.uzfact(0.5) if num == 1 DCL::slmgn(0.0, 0.0, 0.0, 0.15) # 描画マージン指定 GGraph.set_fig('viewport'=>[0.2,0.87,0.12,0.6]) # タイトル #DCL::slsttl('#date #time YAMADA.YU', 'T', 1.0, 1.0, 0.01, 1) DCL::slsttl("#{dir} y=32 composite [t (u,-sigdot)]", 'T', -0.3, -0.8, 0.02, 2) # 経度幅, トーン, コンター設定 lon_width = 18 long_min = 64 - lon_width long_max = 64 + lon_width levels = [-80E-2,-40E-2,0,40E-2,80E-2,120E-2] GGraph.set_tone_levels( 'levels'=>levels, 'patterns'=>[25999,35999,45999,55999,70999,75999,85999] ) GGraph.set_linear_contour_options( 'min'=>-1.2, 'nlev'=>13, 'max'=>1.2 ) # 温度 トーン GGraph.tone( $comp["t"][long_min..long_max,true] ) # 風 (u,-sigdot) ベクトル DCL::uglset('LNRMAL', false) #DCL::ugrset('XFACT1', 0.001) #DCL::ugrset('YFACT1', 10000.0) DCL::ugrset('XFACT1', 3*10E-6*100) DCL::ugrset('YFACT1', 6.0*1000) DCL::uglset('LUNIT', true) DCL::ugrset('VXUNIT', 0.05) DCL::ugrset('VYUNIT', 0.05) DCL::ugrset('VXUOFF', 0.01) DCL::ugvect($comp["u"].val[long_min..long_max,true], - $comp["sigdot"].val[long_min..long_max,true]) GGraph.contour( $comp["t"][long_min..long_max,true] , false) #-- トーンバー ---- x = [-140E-2,-120E-2,-80E-2,-40E-2, 0, 40E-2, 80E-2,120E-2,140E-2] ipat = [25999, 25999, 35999, 45999, 55999, 70999, 75999, 85999] x = NArray.to_na(x) ipat = NArray.to_na(ipat) DCL::ueitlv DCL::uestln(x,ipat) DCL::Util::color_bar("vx0"=>0.05,"vy0"=>0.05,"tick1"=>1,"tick2"=>2,"vxlength"=>0.25,"label_size"=>0.01) unless $graph == 1 then # ps, gif 作成 fig = "figs/t-wind-xzeq-#{$exp}-#{$cum_param}-comp" mradl_dcldump(fig) if num == 1 mradl_dclps(fig) if num == 2 end # お片付け ---- DCL.grcls end # ---------------------------------------------------------- # gif, eps 変換 def mradl_dclps(fig) # gif, ps ファイル生成 puts `mv dcl.ps #{fig}.ps` end def mradl_dcldump(fig) # gif,ファイル生成 puts `xwd -out hoge.xwd` puts `xwdtopnm hoge.xwd > hoge.pnm; ppmtogif hoge.pnm > #{fig}.gif` puts `rm hoge.xwd hoge.pnm` end # ---------------------------------------------------------- # 実行オプション解析 def mradl_optdcl opt = ARGV $graph = 1 ; $cum_param = 0 ; $exp = 0 opt.size.times{ |i| if opt[i] == "-sump" then $graph = 2 end if opt[i] == "-comp" then $graph = 3 end if opt[i] == "-cum" then $cum_param = opt[i+1] end if opt[i] == "-exp" then $exp = opt[i+1] end } end