#!/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 ################################################################ year = [1979, 2003] month = ["01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12"] # 初期値 # 大気上端上向き長波放射 ulwrf_long_name = "OLR Top" ulwrf_tropical_long_name = ulwrf_long_name + " (S30 - N30)" ulwrf_nh_midlat_long_name = ulwrf_long_name + " (N30 - N60)" ulwrf_sh_midlat_long_name = ulwrf_long_name + " (S30 - S60)" # 大気上端正味短波放射 nswrf_long_name = "NSR Top" nswrf_tropical_long_name = nswrf_long_name + " (S30 - N30)" nswrf_nh_midlat_long_name = nswrf_long_name + " (N30 - N60)" nswrf_sh_midlat_long_name = nswrf_long_name + " (S30 - S60)" # 地表面正味長波放射 nlwrs_long_name = "NLR Surface" nlwrs_tropical_long_name = nlwrs_long_name + " (S30 - N30)" nlwrs_nh_midlat_long_name = nlwrs_long_name + " (N30 - N60)" nlwrs_sh_midlat_long_name = nlwrs_long_name + " (S30 - S60)" # 地表面正味短波放射 nswrs_long_name = "NSR Surface" nswrs_tropical_long_name = nswrs_long_name + " (S30 - N30)" nswrs_nh_midlat_long_name = nswrs_long_name + " (N30 - N60)" nswrs_sh_midlat_long_name = nswrs_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 ulwrf_name = ""; ulwrf_unit = "" nswrf_name = ""; nswrf_unit = "" nlwrs_name = ""; nlwrs_unit = "" nswrs_name = ""; nswrs_unit = "" ulwrf_law = []; ulwrf_nh_mid = []; ulwrf_sh_mid = [] nswrf_law = []; nswrf_nh_mid = []; nswrf_sh_mid = [] nlwrs_law = []; nlwrs_nh_mid = []; nlwrs_sh_mid = [] nswrs_law = []; nswrs_nh_mid = []; nswrs_sh_mid = [] year[0].upto(year[1]) do |y| month.each do |m| # NetCDFVar => GPhys path_ulwrf = "../../../ULWRF.NCEP/ULWRF.#{y}.NCEP/ULWRF_#{y}-#{m}_NCEP.nc" # 大気上端上向き長波放射 path_dswrf = "../../../DSWRF.NCEP/DSWRF.#{y}.NCEP/DSWRF_#{y}-#{m}_NCEP.nc" # 大気上端下向き短波放射 path_uswrf = "../../../USWRF.NCEP/USWRF.#{y}.NCEP/USWRF_#{y}-#{m}_NCEP.nc" # 大気上端上向き短波放射 path_nlwrs = "../../../NLWRS.NCEP/NLWRS.#{y}.NCEP/NLWRS_#{y}-#{m}_NCEP.nc" # 地表面正味長波放射 path_nswrs = "../../../NSWRS.NCEP/NSWRS.#{y}.NCEP/NSWRS_#{y}-#{m}_NCEP.nc" # 地表面正味短波放射 gp_ulwrf = GPhys::NetCDF_IO.open(path_ulwrf, 'ulwrf').mean('lon') gp_dswrf = GPhys::NetCDF_IO.open(path_dswrf, 'dswrf').mean('lon') gp_uswrf = GPhys::NetCDF_IO.open(path_uswrf, 'uswrf').mean('lon') gp_nswrf = gp_uswrf - gp_dswrf gp_nlwrs = GPhys::NetCDF_IO.open(path_nlwrs, 'nlwrs').mean('lon') gp_nswrs = GPhys::NetCDF_IO.open(path_nswrs, 'nswrs').mean('lon') # 時系列 VArray 作り ## ulwrf ulwrf_law << ( ( ( gp_ulwrf * gp_gw ).cut( 'lat'=> 30..-30 ) ).sum.val ) ulwrf_nh_mid << ( ( ( gp_ulwrf * gp_gw ).cut( 'lat'=> 60.. 30 ) ).sum.val ) ulwrf_sh_mid << ( ( ( gp_ulwrf * gp_gw ).cut( 'lat'=>-30..-60 ) ).sum.val ) ulwrf_name = gp_ulwrf.data.name ulwrf_unit = gp_ulwrf.data.get_att("units") ## nswrf nswrf_law << ( ( ( gp_nswrf * gp_gw ).cut( 'lat'=> 30..-30 ) ).sum.val ) nswrf_nh_mid << ( ( ( gp_nswrf * gp_gw ).cut( 'lat'=> 60.. 30 ) ).sum.val ) nswrf_sh_mid << ( ( ( gp_nswrf * gp_gw ).cut( 'lat'=>-30..-60 ) ).sum.val ) nswrf_name = gp_nswrf.data.name nswrf_unit = gp_nswrf.data.get_att("units") ## nlwrs nlwrs_law << ( ( ( gp_nlwrs * gp_gw ).cut( 'lat'=> 30..-30 ) ).sum.val ) nlwrs_nh_mid << ( ( ( gp_nlwrs * gp_gw ).cut( 'lat'=> 60.. 30 ) ).sum.val ) nlwrs_sh_mid << ( ( ( gp_nlwrs * gp_gw ).cut( 'lat'=>-30..-60 ) ).sum.val ) nlwrs_name = gp_nlwrs.data.name nlwrs_unit = gp_nlwrs.data.get_att("units") ## nswrs nswrs_law << ( ( ( gp_nswrs * gp_gw ).cut( 'lat'=> 30..-30 ) ).sum.val ) nswrs_nh_mid << ( ( ( gp_nswrs * gp_gw ).cut( 'lat'=> 60.. 30 ) ).sum.val ) nswrs_sh_mid << ( ( ( gp_nswrs * gp_gw ).cut( 'lat'=>-30..-60 ) ).sum.val ) nswrs_name = gp_nswrs.data.name nswrs_unit = gp_nswrs.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 化 # ulwrf を GPhys 化 ## VArray -> GPhys gp_ulwrf_law = GPhys.new(grid, VArray.new(NArray.to_na(ulwrf_law)).rename!(ulwrf_name)) gp_ulwrf_nh_mid = GPhys.new(grid, VArray.new(NArray.to_na(ulwrf_nh_mid)).rename!(ulwrf_name)) gp_ulwrf_sh_mid = GPhys.new(grid, VArray.new(NArray.to_na(ulwrf_sh_mid)).rename!(ulwrf_name)) ## データの属性設定 gp_ulwrf_law.data.set_att('long_name', ulwrf_tropical_long_name) gp_ulwrf_nh_mid.data.set_att('long_name', ulwrf_nh_midlat_long_name) gp_ulwrf_sh_mid.data.set_att('long_name', ulwrf_sh_midlat_long_name) gp_ulwrf_law.data.set_att('units', ulwrf_unit) gp_ulwrf_nh_mid.data.set_att('units', ulwrf_unit) gp_ulwrf_sh_mid.data.set_att('units', ulwrf_unit) # nswrf を GPhys 化 ## VArray -> GPhys gp_nswrf_law = GPhys.new(grid, VArray.new(NArray.to_na(nswrf_law)).rename!(nswrf_name)) gp_nswrf_nh_mid = GPhys.new(grid, VArray.new(NArray.to_na(nswrf_nh_mid)).rename!(nswrf_name)) gp_nswrf_sh_mid = GPhys.new(grid, VArray.new(NArray.to_na(nswrf_sh_mid)).rename!(nswrf_name)) ## データの属性設定 gp_nswrf_law.data.set_att('long_name', nswrf_tropical_long_name) gp_nswrf_nh_mid.data.set_att('long_name', nswrf_nh_midlat_long_name) gp_nswrf_sh_mid.data.set_att('long_name', nswrf_sh_midlat_long_name) gp_nswrf_law.data.set_att('units', nswrf_unit) gp_nswrf_nh_mid.data.set_att('units', nswrf_unit) gp_nswrf_sh_mid.data.set_att('units', nswrf_unit) # nlwrs を GPhys 化 ## VArray -> GPhys gp_nlwrs_law = GPhys.new(grid, VArray.new(NArray.to_na(nlwrs_law)).rename!(nlwrs_name)) gp_nlwrs_nh_mid = GPhys.new(grid, VArray.new(NArray.to_na(nlwrs_nh_mid)).rename!(nlwrs_name)) gp_nlwrs_sh_mid = GPhys.new(grid, VArray.new(NArray.to_na(nlwrs_sh_mid)).rename!(nlwrs_name)) ## データの属性設定 gp_nlwrs_law.data.set_att('long_name', nlwrs_tropical_long_name) gp_nlwrs_nh_mid.data.set_att('long_name', nlwrs_nh_midlat_long_name) gp_nlwrs_sh_mid.data.set_att('long_name', nlwrs_sh_midlat_long_name) gp_nlwrs_law.data.set_att('units', nlwrs_unit) gp_nlwrs_nh_mid.data.set_att('units', nlwrs_unit) gp_nlwrs_sh_mid.data.set_att('units', nlwrs_unit) # nswrs を GPhys 化 ## VArray -> GPhys gp_nswrs_law = GPhys.new(grid, VArray.new(NArray.to_na(nswrs_law)).rename!(nswrs_name)) gp_nswrs_nh_mid = GPhys.new(grid, VArray.new(NArray.to_na(nswrs_nh_mid)).rename!(nswrs_name)) gp_nswrs_sh_mid = GPhys.new(grid, VArray.new(NArray.to_na(nswrs_sh_mid)).rename!(nswrs_name)) ## データの属性設定 gp_nswrs_law.data.set_att('long_name', nswrs_tropical_long_name) gp_nswrs_nh_mid.data.set_att('long_name', nswrs_nh_midlat_long_name) gp_nswrs_sh_mid.data.set_att('long_name', nswrs_sh_midlat_long_name) gp_nswrs_law.data.set_att('units', nswrs_unit) gp_nswrs_nh_mid.data.set_att('units', nswrs_unit) gp_nswrs_sh_mid.data.set_att('units', nswrs_unit) ################################################################ # 描画ルーチン ################################################################ ## # 事前準備 DCL.uscset('cyspos', 'B' ) # y 軸の単位の位置を下方へ rsizel2 = DCL.uzrget('rsizel2') # 現在のラベルサイズを取得 DCL.uzrset('rsizel2', rsizel2*0.42 ) # ラベルサイズをデフォルトの 0.5 倍に ## # お絵描きメイン # 縦軸の最大描画範囲(W/m^2) max = +200 min = -200 # ビューポート設定 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/5, max/2.5) DCL::uyaxdv('R', max/5, max/2.5) DCL::uxmttl('T', 'Mean radiation flux (1979-2003)', 0.0) DCL::uysttl('L', "radiation(#{ulwrf_unit.to_s})", 0.0) DCL::uxsttl('B', 'YEAR', 0.0) # DCL::uscset('cyunit', ulwrf_unit.to_s ) # y 軸の単位 names, idxs = plot_line_main([ gp_ulwrf_law, gp_ulwrf_nh_mid, gp_ulwrf_sh_mid, gp_nswrf_law, gp_nswrf_nh_mid, gp_nswrf_sh_mid, gp_nlwrs_law, gp_nlwrs_nh_mid, gp_nlwrs_sh_mid, gp_nswrs_law, gp_nswrs_nh_mid, gp_nswrs_sh_mid ]) __show_line_index(names,idxs) DCL.grcls ################################################## # 画像ファイル名を変更 File.rename('dcl_001.ps', 'RADIATION_1979-2003_TIME-SERIES_NCEP.ps')