#!/usr/bin/env ruby =begin = 熱領域平均値時系列作図スクリプト(赤道, 中緯度) (e)熱収支( 潜熱フラックス, 顕熱フラックス, 降水 ) =end require "getopts" # for option_parse require "numru/ggraph" include NumRu include Misc::EMath ## 描画メソッドメイン def plot_line_main(gphys_array) line_index_ary = Array.new; name_ary = Array.new default_index = 12 line_hash = { "index"=> default_index } line_index_ary.push(default_index) name_ary.push(gphys_array[0].data.get_att("long_name")) x = gphys_array[0].coord(0).val y = gphys_array[0].data.val DCL.uulinz(x, y, 1, default_index) gphys_array.size.times{ |num| unless num == 0 line_index = (num*10 + default_index) line_index_ary.push(num*10 +default_index) name_ary.push(gphys_array[num].data.get_att("long_name")) x = gphys_array[num].coord(0).val y = gphys_array[num].data.val DCL.uulinz(x, y, 1, line_index) end } return name_ary, line_index_ary end ##----------------------- # ラインインデックスの種類を表示 def __show_line_index(str_ary,line_index_ary, x=0.15, y=0.12) charsize = 1.0 * DCL.uzpget('rsizec1') dvx = 0.01 dvy = charsize*1.5 raise TypeError,"Array expected" if ! str_ary.is_a?(Array) vxmin,vxmax,vymin,vymax = DCL.sgqvpt vx = x vy = y + charsize/2 str_ary.size.times{|num| DCL::sgtxzv(vx,vy,"--- #{str_ary[num]}", charsize, 0, -1, line_index_ary[num]) vy -= dvy if num == 2 vx = 0.30 vy = 0.12 - charsize/2 end } nil end ################################################################ # make gphys ################################################################ L = UNumeric.new(2.5008e6, 'J.Kg-1') year = [1979, 2003] month = ["01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12"] # 初期値 # 潜熱フラックス lhtfl_long_name = "Latent heat flux" lhtfl_tropical_long_name = lhtfl_long_name + " (S30 - N30)" lhtfl_nh_midlat_long_name = lhtfl_long_name + " (N30 - N60)" lhtfl_sh_midlat_long_name = lhtfl_long_name + " (S30 - S60)" # 顕熱フラックス shtfl_long_name = "Sensible heat flux" shtfl_tropical_long_name = shtfl_long_name + " (S30 - N30)" shtfl_nh_midlat_long_name = shtfl_long_name + " (N30 - N60)" shtfl_sh_midlat_long_name = shtfl_long_name + " (S30 - S60)" # 降水 prate_long_name = "Precipitation" prate_tropical_long_name = prate_long_name + " (S30 - N30)" prate_nh_midlat_long_name = prate_long_name + " (N30 - N60)" prate_sh_midlat_long_name = prate_long_name + " (S30 - S60)" ## make compisit-GPhys object # make array which set the gphys objects. ## ガウス重み path_gw = "../../../GAUSSIAN_WEIGHT/gaussian_weight.nc" gp_gw = GPhys::NetCDF_IO.open(path_gw, 'gw') axmonth = [] i = 0 lhtfl_name = ""; lhtfl_unit = "" shtfl_name = ""; shtfl_unit = "" prate_name = ""; prate_unit = "" lhtfl_law = []; lhtfl_nh_mid = []; lhtfl_sh_mid = [] shtfl_law = []; shtfl_nh_mid = []; shtfl_sh_mid = [] prate_law = []; prate_nh_mid = []; prate_sh_mid = [] year[0].upto(year[1]) do |y| month.each do |m| # NetCDFVar => GPhys path_lhtfl = "../../../LHTFL.NCEP/LHTFL.#{y}.NCEP/LHTFL_#{y}-#{m}_NCEP.nc" # 大気上端上向き長波放射 path_shtfl = "../../../SHTFL.NCEP/SHTFL.#{y}.NCEP/SHTFL_#{y}-#{m}_NCEP.nc" # 地表面正味長波放射 path_prate = "../../../PRATE.NCEP/PRATE.#{y}.NCEP/PRATE_#{y}-#{m}_NCEP.nc" # 地表面正味短波放射 gp_lhtfl = GPhys::NetCDF_IO.open(path_lhtfl, 'lhtfl').mean('lon') gp_shtfl = GPhys::NetCDF_IO.open(path_shtfl, 'shtfl').mean('lon') gp_prate = GPhys::NetCDF_IO.open(path_prate, 'prate').mean('lon')*L # 時系列 VArray 作り ## lhtfl lhtfl_law << ( ( ( gp_lhtfl * gp_gw ).cut( 'lat'=> 30..-30 ) ).sum.val ) lhtfl_nh_mid << ( ( ( gp_lhtfl * gp_gw ).cut( 'lat'=> 60.. 30 ) ).sum.val ) lhtfl_sh_mid << ( ( ( gp_lhtfl * gp_gw ).cut( 'lat'=>-30..-60 ) ).sum.val ) lhtfl_name = gp_lhtfl.data.name lhtfl_unit = gp_lhtfl.data.get_att("units") ## shtfl shtfl_law << ( ( ( gp_shtfl * gp_gw ).cut( 'lat'=> 30..-30 ) ).sum.val ) shtfl_nh_mid << ( ( ( gp_shtfl * gp_gw ).cut( 'lat'=> 60.. 30 ) ).sum.val ) shtfl_sh_mid << ( ( ( gp_shtfl * gp_gw ).cut( 'lat'=>-30..-60 ) ).sum.val ) shtfl_name = gp_shtfl.data.name shtfl_unit = gp_shtfl.data.get_att("units") ## prate prate_law << ( ( ( gp_prate * gp_gw ).cut( 'lat'=> 30..-30 ) ).sum.val ) prate_nh_mid << ( ( ( gp_prate * gp_gw ).cut( 'lat'=> 60.. 30 ) ).sum.val ) prate_sh_mid << ( ( ( gp_prate * gp_gw ).cut( 'lat'=>-30..-60 ) ).sum.val ) prate_name = gp_prate.data.name prate_unit = gp_prate.data.get_att("units") axmonth << DCL::dateg3(year[0],1,1,y.to_i,m.to_i,1) i += 1 p y if i%12 == 0 end end ### for intensity # 日付軸作成 namonth = NArray.to_na(axmonth) # 日付軸の NArray vamonth = VArray.new(namonth).rename!("month") # VArray 化. 名前は "month" vamonth.set_att('long_name','month') # 属性設定. long_name と vamonth.set_att('units','month since 1979 JAN') # units を設定. grid = Grid.new( Axis.new(true).set_cell_guess_bounds(vamonth).set_pos_to_center ) # Axis 化して, セル位置を適当に決めて Grid 化 # lhtfl を GPhys 化 ## VArray -> GPhys gp_lhtfl_law = GPhys.new(grid, VArray.new(NArray.to_na(lhtfl_law)).rename!(lhtfl_name)) gp_lhtfl_nh_mid = GPhys.new(grid, VArray.new(NArray.to_na(lhtfl_nh_mid)).rename!(lhtfl_name)) gp_lhtfl_sh_mid = GPhys.new(grid, VArray.new(NArray.to_na(lhtfl_sh_mid)).rename!(lhtfl_name)) ## データの属性設定 gp_lhtfl_law.data.set_att('long_name', lhtfl_tropical_long_name) gp_lhtfl_nh_mid.data.set_att('long_name', lhtfl_nh_midlat_long_name) gp_lhtfl_sh_mid.data.set_att('long_name', lhtfl_sh_midlat_long_name) gp_lhtfl_law.data.set_att('units', lhtfl_unit) gp_lhtfl_nh_mid.data.set_att('units', lhtfl_unit) gp_lhtfl_sh_mid.data.set_att('units', lhtfl_unit) # shtfl を GPhys 化 ## VArray -> GPhys gp_shtfl_law = GPhys.new(grid, VArray.new(NArray.to_na(shtfl_law)).rename!(shtfl_name)) gp_shtfl_nh_mid = GPhys.new(grid, VArray.new(NArray.to_na(shtfl_nh_mid)).rename!(shtfl_name)) gp_shtfl_sh_mid = GPhys.new(grid, VArray.new(NArray.to_na(shtfl_sh_mid)).rename!(shtfl_name)) ## データの属性設定 gp_shtfl_law.data.set_att('long_name', shtfl_tropical_long_name) gp_shtfl_nh_mid.data.set_att('long_name', shtfl_nh_midlat_long_name) gp_shtfl_sh_mid.data.set_att('long_name', shtfl_sh_midlat_long_name) gp_shtfl_law.data.set_att('units', shtfl_unit) gp_shtfl_nh_mid.data.set_att('units', shtfl_unit) gp_shtfl_sh_mid.data.set_att('units', shtfl_unit) # prate を GPhys 化 ## VArray -> GPhys gp_prate_law = GPhys.new(grid, VArray.new(NArray.to_na(prate_law)).rename!(prate_name)) gp_prate_nh_mid = GPhys.new(grid, VArray.new(NArray.to_na(prate_nh_mid)).rename!(prate_name)) gp_prate_sh_mid = GPhys.new(grid, VArray.new(NArray.to_na(prate_sh_mid)).rename!(prate_name)) ## データの属性設定 gp_prate_law.data.set_att('long_name', prate_tropical_long_name) gp_prate_nh_mid.data.set_att('long_name', prate_nh_midlat_long_name) gp_prate_sh_mid.data.set_att('long_name', prate_sh_midlat_long_name) gp_prate_law.data.set_att('units', prate_unit) gp_prate_nh_mid.data.set_att('units', prate_unit) gp_prate_sh_mid.data.set_att('units', prate_unit) ################################################################ # 描画ルーチン ################################################################ ## # 事前準備 DCL.uscset('cyspos', 'B' ) # y 軸の単位の位置を下方へ rsizel2 = DCL.uzrget('rsizel2') # 現在のラベルサイズを取得 DCL.uzrset('rsizel2', rsizel2*0.42 ) # ラベルサイズをデフォルトの 0.5 倍に ## # お絵描きメイン # 縦軸の最大描画範囲(W/m^2) max = +100 min = 0 # ビューポート設定 vpt = NArray[0.15,0.85,0.23,0.58] p idate = (year[0].to_s+"0101").to_i # 日付 0 日. 19790101 ndate = (year[-1].to_s+"1201").to_i # 日付 最後の日. 20031231 p nd = DCL.dateg1(idate,ndate) # 日日数 #GGraph.set_axes("date?"=>["x", idate, nd]) # 日付軸を有効にする DCL.swpset('lsep', true) # ページ毎別々のファイルに落す DCL.gropn(2) DCL.sgpset('lcntl', false) # DCL.sgpset('lfull',true) # フルスクリーン DCL.uzfact(0.55) # DCL.sgpset('lcorner',false) # コーナーを取っ払う DCL.sgpset('lfprop',true) # DCL.udpset('lmsg',false) # DCL.uscset('cyspos', 'B' ) # y 軸の単位の位置を下方へ ### 時系列 DCL::grfrm # ページ確定 DCL::grswnd(0.0, nd, min, max) DCL::grsvpt(*vpt) DCL::grstrn(1) DCL::grstrf DCL::ucxacl('B', idate, nd) DCL::ucxacl('T', idate, nd) DCL::uyaxdv('L', (max - min)/10, (max - min)/5) DCL::uyaxdv('R', (max - min)/10, (max - min)/5) DCL::uxmttl('T', 'Mean heat flux (1979-2003)', 0.0) DCL::uysttl('L', "heat flux(#{lhtfl_unit.to_s})", 0.0) DCL::uxsttl('B', 'YEAR', 0.0) # DCL::uscset('cyunit', lhtfl_unit.to_s ) # y 軸の単位 names, idxs = plot_line_main([ gp_lhtfl_law, gp_lhtfl_nh_mid, gp_lhtfl_sh_mid, gp_shtfl_law, gp_shtfl_nh_mid, gp_shtfl_sh_mid, gp_prate_law, gp_prate_nh_mid, gp_prate_sh_mid ]) __show_line_index(names,idxs) DCL.grcls ################################################## # 画像ファイル名を変更 File.rename('dcl_001.ps', 'HEAT_1979-2003_TIME-SERIES_NCEP.ps')