#!/usr/bin/env ruby # # 表題: 熱帯降水活動コンポジット解析 # # 履歴: 2003/11/15 やまだ由 # 履歴: 2003/11/18 やまだ由 # # ------------------------- require "numru/ggraph" require "colorbar.rb" include NumRu include NMath # ------------------------- # 実験名 $cum_param = "adj" $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 $cum_param = "adj" $exp = "con" $graph = 1 # ディレクトリ名 $dir = "#{$cum_param}700#{$exp}" # 出力ファイル名 # $outfile = "#{$dir}_composite_rain.nc" # $outfile = "#{$dir}_composite_rain_westward.nc" $outfile = "#{$dir}_composite_rain_eastward.nc" composite(threshold) } # ------------------------- def composite(threshold) # 変数初期化 data = {} ; $comp = {} comp_num = 0 ; count = 0 comp_point = NArray.sfloat(128,400).fill!(0.0) # コンポジットを評価する rain 初期化 rain_data = GPhys::NetCDF_IO.open("#{$filepath}/#{$dir}/rain_sum_seq.nc",'rain') rain_data_cal = NArray.sfloat(128+6).fill!(0.0) grid_lon = rain_data.axis(0) grid_time = rain_data.axis(1) # 軸の値を保存 data_ps = GPhys::NetCDF_IO.open("#{$filepath}/#{$dir}/ps_out.nc",'ps') grid_lat = data_ps.axis(1) # コンポジットをとる変数の指定 # item = ["sigdot","t", "qcndc","ps" ] # item = ["sigdot","t","u" ] item_xy = ["rain"] item_xy.each{ |i| if i == "rain" then data[i] = GPhys::NetCDF_IO.open("#{$filepath}/#{$dir}/#{i}_sum.nc",i) else data[i] = GPhys::NetCDF_IO.open("#{$filepath}/#{$dir}/#{i}_z1_out.nc",i) end data[i] = data[i].val[true,true,true] p i # コンポジット値初期化 $comp[i] = NArray.sfloat(128,64).fill!(0.0) } # ---------------------------------------------------------------- comp_sel = [ #-----------------------------------------------# #---- adj-con 東進 -----------------------------# #-----------------------------------------------# #=begin [287,18], [288,17], [289,21], [290,20], [295,26], [296,26], [297,26], [327,61],[328,60], [331,65],[332,64],[333,64],[334,64], [335,65],[336,65], [344,72],[345,72], [345,77],[346,71], [346,76],[347,76], [348,75],[349,75], [349,81],[350,81],[351,81] #=end #-----------------------------------------------# #---- adj-con 西進 -----------------------------# #-----------------------------------------------# =begin [55,92], [54,93], [53,93], [52,94], [51,95], [47,98],[46,99],[45,99], [25,29], [26,28],[27,27], [10,49],[11,48],[12,47],[13,47],[14,46], [107,1], [106,2],[105,3],[104,3],[103,4] , [158,111],[157,112],[156,113], # ナナメ [139,119],[138,120], [137,121],[136,121],[135,122],[134,123], [245,94],[244,95],[243,97],[242,98] , [233,100], [231,101],[230,102] =end # 弱い降水の西進 =begin [182,69], [188,69], [189,68], [190,68], [191,67], [192,66], [193,68], [194,67], [195,67], [196,66], [197,66], [198,65], [199,65], [200,64], [201,63], [202,62], [203,62], [204,61], [205,61], [206,60], [207,62], [208,62], [209,61] =end ] # comp_num: 経度座標, time: 時間座標 comp_sel.each{ |num| comp_num = num[1] time = num[0] comp_point[comp_num,time] = 1.0 # n = 63 - comp_num n = 64 - comp_num item_xy.each{ |i| tmp = NArray.sfloat(128,64).fill!(0.0) if n == 0 then tmp = data[i][true,true,time] elsif n > 0 then tmp[0..(n-1),true] = data[i][(128-n)..127,true,time] tmp[n..127,true] = data[i][0..(127-n),true,time] else tmp[(128+n)..127,true] = data[i][0..(-1-n),true,time] tmp[0..(127+n),true] = data[i][(-n)..127,true,time] end $comp[i] = $comp[i] + tmp } count = count + 1 } # 平均 ------------------------------------------ puts count item_xy.each{ |i| $comp[i] = $comp[i] / count # v, u, ps は平均からの偏差 composite_sub_xy(i) if i == "u" || i == "v" || i == "ps" } # Gphys オブジェクト生成 & nc ファイル生成 ---- # 出力ファイル ncfile = NetCDF.create("#{$filepath_out}/#{$outfile}") # 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_lat) item_xy.each{ |i| if i == "rain" then $comp[i] = VArray.new($comp[i]).rename(i) else $comp[i] = VArray.new($comp[i]).rename("#{i}_anm") end $comp[i] = GPhys.new(grid, $comp[i]) GPhys::NetCDF_IO.write( ncfile, $comp[i] ) } ncfile.close end # ------------------------- def composite_sub_xy(i) # 変数初期化 data_mean = {} 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,exp) #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_rain(num,exp) # 絵を描く変数の指定 item = ["rain" ] $comp = {} #ncfile 名 ncfile = $outfile item.each{ |i| # 変数の読み込み $comp[i] = GPhys::NetCDF_IO.open("#{$filepath}composite_nc/#{ncfile}", i) } # 設定 ---------- 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} composite [rain]", 'T', -0.3, -0.8, 0.02, 2) # 経度幅, トーン, コンター設定 lon_width = 18 long_min = 64 - lon_width long_max = 64 - 1 + lon_width lat_width = 9 lat_min = 32 - lat_width lat_max = 32 - 1 + lat_width levels = [0,200,400,600,800.1000,1200,1400] GGraph.set_tone_levels( 'levels'=>levels, 'patterns'=>[25999,35999,40999,55999,70999,75999,85999] ) GGraph.set_linear_contour_options( 'min'=>0.0, 'nlev'=>10.0 , 'max'=>1000.0 ) # rain トーン GGraph.tone( $comp["rain"][long_min..long_max,lat_min..lat_max] ) GGraph.contour( $comp["rain"][long_min..long_max,lat_min..lat_max] , false) #-- トーンバー ---- x = [ 0, 200, 400, 600, 800, 1000, 1200, 1400 ] ipat = [25999, 35999, 40999, 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/rain-#{$exp}-#{$cum_param}-comp" # fig = "figs/rain-#{$exp}-#{$cum_param}-east-comp" fig = "figs/rain-#{$exp}-#{$cum_param}-west-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 dcl_gph_wind # tmp_u = NArray.sfloat(128,64).fill!(0.0) # tmp_v = NArray.sfloat(128,64).fill!(0.0) 128.times{ |num| unless num % 4 == 0 then # 現在未使用の変数…. tmp_u[num,true] = $comp["u"][num,true] tmp_v[num,true] = $comp["v"][num,true] end } $comp["u"] = u $comp["v"] = v 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