#!/usr/bin/ruby # -*- coding: utf-8-with-signature -*- # =begin # = 全球, 昼半球, 夜半球で平均した 3 次元物理量を出力する. # # == 説明 # 3 次元物理量の全球, 昼半球, 夜半球平均を図に出力する. # 整数レベルにある物理量のみ. # sigma 面を指定する. # # == オプション # --time_start 開始時刻 # --time_end 終了時刻 # --sig 鉛直座標値 (sigma) # --wsn 出力先. ps だったら --wsn=2 # # == 使用例 # % omega_QVap.rb --time_start=731 --time_end=1096 --sig=0.1 --wsn=2 # # == 更新履歴 # # * 2013/06/25 石渡 正樹 変更 # * 2011/04/24 納多 哲史 新規作成 # # =end require "getoptlong" # for option_parse require "numru/ggraph" include NumRu require "./load_config.rb" # load_config.rb を利用 #DATA_HOME = "/home/noda/data/" TMPFILE = "dcl.ps" ## オプション解析 parser = GetoptLong.new parser.set_options( ### global option ### ['--time_start', GetoptLong::REQUIRED_ARGUMENT], ['--time_end', GetoptLong::REQUIRED_ARGUMENT], ['--sig', GetoptLong::REQUIRED_ARGUMENT], ['--wsn', GetoptLong::REQUIRED_ARGUMENT] ) parser.each_option do |name, arg| eval "$OPT_#{name.sub(/^--/, '').gsub(/-/, '_')} = '#{arg}'" # strage option value to $OPT_val end wsn = ($OPT_wsn||4) sig = ($OPT_sig||1).to_f time_start = ($OPT_time_start||1).to_f time_end = ($OPT_time_end||9999).to_f MARK_SIZE = 0.01 # 現時点では 2 枚以上書くとエラーまたは図がおかしくなるので 1 枚ずつ vars = [ # {'name'=>'Evap', 'min'=>0, 'max'=>600}, # {'name'=>'Rain', 'min'=>0, 'max'=>500}, # {'name'=>'OLR', 'min'=>270, 'max'=>400}, # {'name'=>'Sens', 'min'=>-20, 'max'=>100}, # {'name'=>'SLR', 'min'=>0, 'max'=>200} # {'name'=>'QVap', 'min'=>240, 'max'=>320}, {'name'=>'QVap', 'min'=>1e-6, 'max'=>1e-4}, # {'name'=>'Evap', 'min'=>-50, 'max'=>600}, # {'name'=>'Rain', 'min'=>-50, 'max'=>600}, # {'name'=>'OLR', 'min'=>-50, 'max'=>600}, # {'name'=>'Sens', 'min'=>-50, 'max'=>600}, # {'name'=>'SLR', 'min'=>-50, 'max'=>600} ] # load_config.rb を利用 #omegas = [0.0, ## 0.01, 0.02, 0.03, # 0.05, ## 0.06, 0.07, # 0.1, 0.15, 0.2, 0.25, # 0.33, 0.4, 0.5, 0.6, 0.67, 0.75, 0.8, 0.85, 1.0] ## , 1.2, 1.5, 2.0, 3.0] #omegas = [0.0, 1.0] na_omegas = NArray.to_na(omegas) nomega = omegas.size # 各変数で回す for v in vars var = v['name'] epsfn = var + '.eps' # if File.exist?(epsfn) then # puts(epsfn + " exists. process skipped.") # next # end all = NArray.sfloat(nomega) west = NArray.sfloat(nomega) east = NArray.sfloat(nomega) ## 各 omega で回す. for i in 0..nomega-1 omega = na_omegas[i] var.to_s + ' ' + omega.to_s ### 読み込み dir = 'SR_Omega' + omega.to_s file = var + '.nc' p path = File.join(DATA_HOME, dir, file) gphys = GPhys::IO.open(path, var) na_lon = GPhys::IO.open(path, 'lon').val na_lat = GPhys::IO.open(path, 'lat').val na_lat_weight = GPhys::IO.open(path, 'lat_weight').val nlon = na_lon.size nlat = na_lat.size ### 平均の計算 #### 時間平均 gphys = gphys.cut('time'=>time_start..time_end,'sig'=>sig).mean('time') #### 空間平均 ##### 経度方向の重み付け. gphys = gphys * na_lat_weight.reshape(1, nlat) / na_lat_weight.sum all[i] = gphys.mean('lon').sum('lat').val west[i] = gphys.cut('lon'=>na_lon[0]..na_lon[nlon/2 -1]).mean('lon').sum('lat').val east[i] = gphys.cut('lon'=>na_lon[nlon/2]..na_lon[nlon-1]).mean('lon').sum('lat').val ## check # p all - (west + east)/2.0 end ## gphys オブジェクトの生成 # va_omega = VArray.new( na_omegas, ## {"long_name"=>"Omega scaled by the Earth's value","units"=>"1"}, # {"long_name"=>"Omega","units"=>"1"}, # "omega" ) # (2013-5-28 石渡編集) # Omega 軸の long_name は Omega/Omega_E にしたかったのだが, # omega_heat-wat_atmos-surf_by_Watt.rb でうまくいかなかったので # しかたなくこちらでも Omega/OmegaE にした. va_omega = VArray.new( na_omegas, {"long_name"=>"Omega/OmegaE"}, "omega" ) axis_omega = Axis.new.set_pos(va_omega) data_all = VArray.new( all, {"long_name"=>gphys.data.long_name, "units"=>gphys.data.units.to_s}, gphys.data.name ) gp_all = GPhys.new( Grid.new(axis_omega), data_all ) data_east = VArray.new( east, {"long_name"=>gphys.data.long_name, "units"=>gphys.data.units.to_s}, gphys.data.name ) gp_east = GPhys.new( Grid.new(axis_omega), data_east ) data_west = VArray.new( west, {"long_name"=>gphys.data.long_name, "units"=>gphys.data.units.to_s}, gphys.data.name ) gp_west = GPhys.new( Grid.new(axis_omega), data_west ) ## 描画 # 描画準備 DCL.gropn(wsn) #GGraph.set_fig('viewport'=>[0.15,0.80,0.15,0.6]) #GGraph.set_fig('viewport'=>[0.3,0.6,0.3,0.5]) GGraph.set_fig('itr'=> 2, 'viewport'=>[0.3,0.6,0.3,0.5]) #GGraph.set_axes('ytitle'=>y_name) DCL.sgpset('lfull', true) # 全画面表示 DCL.sgpset('lcntl', false) DCL.uzfact(0.6) # 初期化? # DCL.grfrm #GGraph.set_axes('xtickint'=>xlabel_int) #GGraph.set_axes('xlabelint'=>xlabel_int) #GGraph.set_axes('ytickint'=>ylabel_int) #GGraph.set_axes('ylabelint'=>ylabel_int) #GGraph.line( 'max'=>val_max, 'min'=>val_min ) # opt = { # 'max'=>v['max'], 'min'=>v['min'], \ # 'legend'=>false, 'annotate'=>false # # 'title'=>y_name # } # (2013-5-28 石渡編集) opt = { 'max'=>v['max'], 'min'=>v['min'], \ 'legend'=>false, 'annotate'=>false, \ 'title'=>'' } ## line # opt['index'] = 7 # opt['type'] = 1; GGraph.line( gp_all, true, opt ) ## mark opt['index'] = 5; opt['size'] = MARK_SIZE opt['type'] = 4 # circle GGraph.mark( gp_all, true, opt ) opt['type'] = 5 # cross GGraph.mark( gp_west, false, opt ) opt['type'] = 6 # square GGraph.mark( gp_east, false, opt ) # GGraph.mark( gp_east, false, opt ) # ad hoc DCL.grcls # remove extra space of ps file if wsn == "2" then if File.exist?(TMPFILE) then puts("MESSAGE: after processing ... #{epsfn} will be generate.") `dclpsrot #{TMPFILE} > #{epsfn}` end end end exit(0)