#!/usr/bin/env ruby # # 表題: 熱帯降水活動コンポジット解析 # # 履歴: 2003/11/15 やまだ由 # 履歴: 2003/11/18 やまだ由 # 履歴: 2004/03/18 やまだ由 # # ------------------------- require "numru/ggraph" include NumRu include NMath # ------------------------- class GPhys def composite(rain, threshold=0.0005, flag=false ) # gphys : コンポジットをとる変数 (lon,sigmap,time) # rain : コンポジットを評価する rain データ (lon,time) # flag=true : コンポジットの参照点 (lon,time) を返す # flag=false : gphys のコンポジットの結果 (lon,sigmap) を返す # threshold : rain 閾値 # 閾値 # threshold = 0.0005 # 軸の値を保存 grid = @grid.copy grid_lon = grid.axis(0) grid_sigmap = grid.axis(1) grid_time = grid.axis(2) # 軸のサイズを保存 grid_lon_size = grid.axis(0).pos.val.size grid_sigmap_size = grid.axis(1).pos.val.size grid_time_size = grid.axis(2).pos.val.size # 作業領域初期化 rain_cal = NArray.sfloat( rain.coord(0).val.size + 6 ).fill!(0.0) comp_num = 0 ; count = 0 comp_point = NArray.sfloat(grid_lon_size,grid_time_size).fill!(0.0) # データ data = @data.copy.val # コンポジット値初期化 comp = NArray.sfloat(grid_lon_size , grid_sigmap_size ).fill!(0.0) grid_time_size.times{ |time| # コンポジット格子値計算用 rain データ作成 rain_cal[3..-4] = rain.val[true,time] # 重ね合わせ (rain_cal.size-6).times{ |comp_num| # 周囲よりも値が大きければその点を中心経度に移動 if rain_cal[comp_num+3] > threshold then num = comp_num + 3 if rain_cal[num] > rain_cal[num - 3] && rain_cal[num] > rain_cal[num - 2] && rain_cal[num] > rain_cal[num - 1] && rain_cal[num] > rain_cal[num + 1] && rain_cal[num] > rain_cal[num + 2] && rain_cal[num] > rain_cal[num + 3] then comp_point[comp_num,time] = 1.0 print '(', time, ',', comp_num, ')', "\n" puts rain_cal[num] n = grid_lon_size/2 - comp_num tmp = NArray.sfloat(grid_lon_size,grid_sigmap_size).fill!(0.0) nm = grid_lon_size if n == 0 then tmp = data[true,true,time] elsif n > 0 then tmp[0..(n-1),true] = data[(nm-n)..(nm-1),true,time] tmp[n..(nm-1),true] = data[0..(nm-1-n),true,time] else tmp[(nm+n)..(nm-1),true] = data[0..(-1-n),true,time] tmp[0..(nm-1+n),true] = data[(-n)..(nm-1),true,time] end comp = comp + tmp count = count + 1 end end } } # 平均 comp = comp / count # Gphys オブジェクト生成 ---- # 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) # composite した変数 grid = Grid.new(grid_lon, grid_sigmap) comp = VArray.new(comp).rename("#{@data.copy.name}_comp") comp = GPhys.new(grid, comp) if flag return comp_point else return comp end end def composite_online_old(rain, threshold=0.0005, flag=false ) # gphys : コンポジットをとる変数 (lon,sigmap,time) # rain : コンポジットを評価する rain データ (lon,time) # flag=true : コンポジットの参照点 (lon,time) を返す # flag=false : gphys のコンポジットの結果 (lon,sigmap) を返す # threshold : rain 閾値 # 閾値 # threshold = 0.0005 # 軸の値を保存 grid = @grid.copy grid_lon = grid.axis(0) grid_sigmap = grid.axis(1) grid_time = grid.axis(2) # 軸のサイズを保存 grid_lon_size = grid.axis(0).pos.val.size grid_sigmap_size = grid.axis(1).pos.val.size grid_time_size = grid.axis(2).pos.val.size # 作業領域初期化 comp_num = 0 ; count = 0 comp_point = NArray.sfloat(grid_lon_size,grid_time_size).fill!(0.0) # データ data = @data.copy.val # コンポジット値初期化 comp = NArray.sfloat(grid_lon_size , grid_sigmap_size ).fill!(0.0) # grid_time_comp = grid_lon.pos.val/13+(60.0+180.0/13) grid_lon.pos.val.each{ |comp_num| =begin # T159L48_eml # if comp_num > 0 # if comp_num < 100 if comp_num > 96 if comp_num < 97 # if comp_num > 45 # if comp_num < 46 # if comp_num > -40 # if comp_num < 10 # if comp_num > -180 # if comp_num < -100 # 経度と時間 time = comp_num*1.0/13+(61.0+180.0/13) ## time = comp_num*1.0/13+(28.0+180.0/13) ## time = comp_num*1.0/13+(59.0+180.0/13) ## time = comp_num*1.0/13+(60.0+180.0/13) ## time = comp_num*1.0/13+(59.0+180.0/13) ## time = comp_num*1.0/13+(62.0+180.0/13) =end =begin # T159L48_non : wave-CISK # if comp_num > -50 # if comp_num < 50 if comp_num > 0 if comp_num < 70 # 経度と時間 # time = comp_num*11.0/120 +(51 + 180.0*11/120) time = comp_num*11.0/120 +(24 + 180.0*11/120) =end if comp_num > 0 if comp_num < 70 time = comp_num*11.0/120 +(51 + 180.0*11/120) print '(',time, ',', comp_num, ')', "\n" # 座標値に直す time = (time*4).to_i comp_num = comp_num*grid_lon_size/360 + grid_lon_size/2 print '(',time, ',', comp_num, ')', "\n" # コンポジットの参照点 comp_point[comp_num,time] = 1.0 # 作業配列 tmp = NArray.sfloat(grid_lon_size,grid_sigmap_size).fill!(0.0) n = ((grid_lon_size/2) - comp_num ).to_i nm = grid_lon_size if n == 0 then tmp = data[true,true,time] elsif n > 0 then tmp[0..(n-1),true] = data[(nm-n)..(nm-1),true,time] tmp[n..(nm-1),true] = data[0..(nm-1-n),true,time] else tmp[(nm+n)..(nm-1),true] = data[0..(-1-n),true,time] tmp[0..(nm-1+n),true] = data[(-n)..(nm-1),true,time] end comp = comp + tmp count = count + 1 end end } # 平均 comp = comp / count # Gphys オブジェクト生成 ---- # 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) # composite した変数 grid = Grid.new(grid_lon, grid_sigmap) comp = VArray.new(comp).rename("#{@data.copy.name}_comp") comp = GPhys.new(grid, comp) if flag return comp_point else return comp end end def composite_online(rain, threshold=0.0005, flag=false ) # 軸の値を保存 grid = @grid.copy grid_lon = grid.axis(0) grid_sigmap = grid.axis(1) grid_time = grid.axis(2) # 軸のサイズを保存 grid_lon_size = grid.axis(0).pos.val.size grid_sigmap_size = grid.axis(1).pos.val.size grid_time_size = grid.axis(2).pos.val.size # 作業領域初期化 comp_num = 0 ; count = 0 comp_point = NArray.sfloat(grid_lon_size,grid_time_size).fill!(0.0) # データ data = @data.copy.val # コンポジット値初期化 comp = NArray.sfloat(grid_lon_size , grid_sigmap_size ).fill!(0.0) =begin # ------------------------------------- grid_lon.pos.val.each{ |comp_num| if comp_num > 55 if comp_num < 75 time = - comp_num*22.0/140 + 91 comp,comp_point,count = composite_calc(time,comp_num,data,comp,comp_point,count) end end } # ------------------------------------- grid_lon.pos.val.each{ |comp_num| if comp_num > 82 if comp_num < 93 time = - comp_num*22.0/140 + 97 comp,comp_point,count = composite_calc(time,comp_num,data,comp,comp_point,count) end end } # ------------------------------------- grid_lon.pos.val.each{ |comp_num| if comp_num > 90 if comp_num < 103 time = - comp_num*22.0/140 + 98.5 comp,comp_point,count = composite_calc(time,comp_num,data,comp,comp_point,count) end end } =end # ------------------------------------- grid_lon.pos.val.each{ |comp_num| if comp_num > -50 if comp_num < 50 time = comp_num*11.0/120 +(51 + 180.0*11/120) comp,comp_point,count = composite_calc(time,comp_num,data,comp,comp_point,count) end end } # 平均 comp = comp / count # Gphys オブジェクト生成 ---- # 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) # composite した変数 grid = Grid.new(grid_lon, grid_sigmap) comp = VArray.new(comp).rename("#{@data.copy.name}_comp") comp = GPhys.new(grid, comp) if flag return comp_point else return comp end end def composite_calc(time,comp_num,data,comp,comp_point,count) # 軸の値を保存 grid = @grid.copy grid_lon = grid.axis(0) grid_sigmap = grid.axis(1) grid_time = grid.axis(2) # 軸のサイズを保存 grid_lon_size = grid.axis(0).pos.val.size grid_sigmap_size = grid.axis(1).pos.val.size grid_time_size = grid.axis(2).pos.val.size print '(',time, ',', comp_num, ')', "\n" time = (time*4).to_i comp_num = comp_num*grid_lon_size/360 + grid_lon_size/2 print '(',time, ',', comp_num, ')', "\n" comp_point[comp_num,time] = 1.0 # 作業配列 tmp = NArray.sfloat(grid_lon_size,grid_sigmap_size).fill!(0.0) n = ((grid_lon_size/2) - comp_num ).to_i nm = grid_lon_size if n == 0 then tmp = data[true,true,time] elsif n > 0 then tmp[0..(n-1),true] = data[(nm-n)..(nm-1),true,time] tmp[n..(nm-1),true] = data[0..(nm-1-n),true,time] else tmp[(nm+n)..(nm-1),true] = data[0..(-1-n),true,time] tmp[0..(nm-1+n),true] = data[(-n)..(nm-1),true,time] end comp = comp + tmp count = count + 1 return comp,comp_point,count end end # ----------------------------------------------------------------- class Ape_mkfig def mkfig_plot_comp(gphys, gphys_comp) if @fig_num == 3 then @data.mkdump_comp(gphys, gphys_comp) else @data.plot_comp(gphys, gphys_comp) end end def gr_comp_tend file_name = Array.new Dir.foreach($gr2ncfile_path) { |item| file_name.push("#{$gr2ncfile_path}/#{item}") if item =~ /.nc$/ } @data = Ape.new(file_name) # tr_tppn rain = @data.gphys_open("tr_tppn").cut(true,0,true).lon_lotate lost_axes = rain.lost_axes.to_s.sub("y=","lat=") rain = rain[true,-400..-1].set_lost_axes(lost_axes) # t_conv composite gphys = @data.gphys_open("pf_t_conv")[true,true,-400..-1].lon_lotate comp = gphys.composite(rain, 0.0007) comp = comp. set_att("ape_name", "cumulus_heating_composite" ). set_att("units", gphys.data.get_att("units") ). rename("comp_t_conv"). set_lost_axes(lost_axes) # 軸の long_name 直し (pressure level => sigma) z = comp.grid_copy.axis(1).pos. set_att("long_name","sigma"). set_att("units","1") z = Axis.new().set_pos(z) grid = comp.grid_copy.change_axis(1, z) comp = GPhys.new( grid, comp.data ) # 描画 dim = comp.axis(0).pos.val.size mkfig_plot(comp[(dim/3)..(dim*2/3),true]) comp = comp.cut(0,true).rename("comp_t_conv_zonal") mkfig_plot(comp) end def gr_comp_point file_name = Array.new Dir.foreach($gr2ncfile_path) { |item| file_name.push("#{$gr2ncfile_path}/#{item}") if item =~ /.nc$/ } @data = Ape.new(file_name) # tr_tppn rain = @data.gphys_open("tr_tppn").cut(true,0,true).lon_lotate lost_axes = rain.lost_axes.to_s.sub("y=","lat=") rain = rain[true,-401..-1] grid_1 = Axis.new(). set_pos(VArray.new(NArray.sfloat(401).indgen!/4).rename("time"). put_att("units","days")) grid_0 = rain.grid_copy.axis(0) grid = Grid.new(grid_0,grid_1) rain = GPhys.new(grid, rain.data ).set_lost_axes(lost_axes) # tr_u250 u250 = @data.gphys_open("U").cut(true,0.25,true).lon_lotate lost_axis = [u250.data.get_att("lost_axis"), u250.lost_axes.to_s.sub("z=","sigma=").sub(" hPa","")] u250 = u250[true,-401..-1].rename("tr_u250") u250 = ( u250 - u250.mean(0) ) grid_1 = Axis.new(). set_pos(VArray.new(NArray.sfloat(401).indgen!/4).rename("time"). put_att("units","days")) grid_0 = u250.grid_copy.axis(0) grid = Grid.new(grid_0,grid_1) u250 = GPhys.new(grid, u250.data ).set_lost_axes(lost_axis). add_lost_axes("(diff) from (mean) zonal") # t_conv composite gphys = @data.gphys_open("T")[true,true,-401..-1].lon_lotate comp = gphys.composite(rain, 0.0007, true) comp = comp. set_att("ape_name", "composite_point" ). set_att("units", "1" ). rename("comp_point_tr_tppn"). # rename("comp_point"). set_lost_axes(lost_axes) # set_lost_axes("") grid_1 = Axis.new(). set_pos(VArray.new(NArray.sfloat(401).indgen!/4).rename("time"). put_att("units","days")) grid_0 = comp.grid_copy.axis(0) grid = Grid.new(grid_0,grid_1) comp = GPhys.new(grid, comp.data ) # 描画 # mkfig_plot(comp) mkfig_plot_comp(rain, comp) comp = comp.rename("comp_point_tr_u250") mkfig_plot_comp(u250, comp) end def gr_comp_tuom file_name = Array.new Dir.foreach($gr2ncfile_path) { |item| file_name.push("#{$gr2ncfile_path}/#{item}") if item =~ /.nc$/ } @data = Ape.new(file_name) # tr_tppn rain = @data.gphys_open("tr_tppn").cut(true,0,true).lon_lotate lost_axes = rain.lost_axes.to_s.sub("y=","lat=") rain = rain[true,-401..-1].set_lost_axes(lost_axes) # T composite gphys = @data.gphys_open("T")[true,true,-401..-1].lon_lotate t_comp = gphys.composite(rain, 0.0007) t_comp = t_comp. set_att("ape_name", "T_(U,_-OMG)_composite"). set_att("units", "K, (m s-1, Pa s-1)" ) t_comp = (t_comp - t_comp.mean(0)). set_lost_axes(lost_axes). add_lost_axes("(diff) from (mean) zonal"). rename("comp_tuom") # Tv composite gphys = @data.gphys_open("TV")[true,true,-401..-1].lon_lotate tv_comp = gphys.composite(rain, 0.0007) tv_comp = tv_comp. set_att("ape_name", "Tv_(U,_-OMG)_composite"). set_att("units", "K, (m s-1, Pa s-1)" ) tv_comp = (tv_comp - tv_comp.mean(0)). set_lost_axes(""). set_lost_axes(lost_axes). add_lost_axes("(diff) from (mean) zonal"). rename("comp_tvuv") # Q composite gphys = @data.gphys_open("Q")[true,true,-401..-1].lon_lotate q_comp = gphys.composite(rain, 0.0007) q_comp = q_comp. set_att("ape_name", "Q_(U,_-OMG)_composite"). set_att("units", "kg/kg, (m s-1, Pa s-1)" ) q_comp = (q_comp - q_comp.mean(0)). set_lost_axes(""). set_lost_axes(lost_axes). add_lost_axes("(diff) from (mean) zonal"). rename("comp_quv") # U composite gphys = @data.gphys_open("U")[true,true,-400..-1].lon_lotate u_comp = gphys.composite(rain, 0.0007) u_comp = u_comp - u_comp.mean(0) # OMG composite gphys = @data.gphys_open("OMG")[true,true,-400..-1].lon_lotate om_comp = gphys.composite(rain, 0.0007) om_comp = om_comp - om_comp.mean(0) # ベクトルの間引 u_data = NArray.sfloat( u_comp.axis(0).pos.val.size , u_comp.axis(1).pos.val.size ).fill!(0.0) om_data = NArray.sfloat( om_comp.axis(0).pos.val.size , om_comp.axis(1).pos.val.size ).fill!(0.0) num_l = u_comp.axis(0).pos.val.size/120 # num_l = u_comp.axis(0).pos.val.size/120*3 num_z = u_comp.axis(1).pos.val.size/48 u_comp.axis(0).pos.val.size.times{ |num| if (num % num_l) == 0 u_data[num,true] = u_comp.data.val[num,true] om_data[num,true] = om_comp.data.val[num,true] elsif u_data[num,true] = 0.0 om_data[num,true] = 0.0 end } num_z =1 if num_z < 1 u_comp.axis(1).pos.val.size.times{ |num| if (num % num_z) == 0 u_data[true,num] = u_data[true,num] om_data[true,num] = om_data[true,num] elsif u_data[true,num] = 0.0 om_data[true,num] = 0.0 end } # 軸の long_name 直し (pressure level => sigma) z = t_comp.grid_copy.axis(1).pos. set_att("long_name","sigma"). set_att("units","1") z = Axis.new().set_pos(z) grid = t_comp.grid_copy.change_axis(1, z) # 配列数取得 dim = t_comp.axis(0).pos.val.size =begin # 描画 (T) t_comp = GPhys.new( grid, t_comp.data ) mkfig_plot(t_comp, u_data, -om_data) # 描画 (Tv) tv_comp = GPhys.new( grid, tv_comp.data ) mkfig_plot(tv_comp, u_data, -om_data) # 描画 (Q) q_comp = GPhys.new( grid, q_comp.data ) mkfig_plot(q_comp, u_data, -om_data) =end # 描画 (T) t_comp = GPhys.new( grid, t_comp.data ) mkfig_plot(t_comp[(dim/3)..(dim*2/3),true], u_data[(dim/3)..(dim*2/3),true], -om_data[(dim/3)..(dim*2/3),true]) # 描画 (Tv) tv_comp = GPhys.new( grid, tv_comp.data ) mkfig_plot(tv_comp[(dim/3)..(dim*2/3),true], u_data[(dim/3)..(dim*2/3),true], -om_data[(dim/3)..(dim*2/3),true]) # 描画 (Q) q_comp = GPhys.new( grid, q_comp.data ) mkfig_plot(q_comp[(dim/3)..(dim*2/3),true], u_data[(dim/3)..(dim*2/3),true], -om_data[(dim/3)..(dim*2/3),true]) end def gr_comp_tuom_online file_name = Array.new Dir.foreach($gr2ncfile_path) { |item| file_name.push("#{$gr2ncfile_path}/#{item}") if item =~ /.nc$/ } @data = Ape.new(file_name) # tr_tppn rain = @data.gphys_open("tr_tppn").cut(true,0,true).lon_lotate lost_axes = rain.lost_axes.to_s.sub("y=","lat=") rain = rain[true,-401..-1].set_lost_axes(lost_axes) # rain = rain[true,-401..-1].set_lost_axes("") # T composite gphys = @data.gphys_open("T")[true,true,-401..-1].lon_lotate t_comp = gphys.composite_online(rain, 0.0007) t_comp = t_comp. set_att("ape_name", "T_(U,_-OMG)_composite"). set_att("units", "K, (m s-1, Pa s-1)" ) t_comp = (t_comp - t_comp.mean(0)). set_lost_axes(""). set_lost_axes(lost_axes). add_lost_axes("(diff) from (mean) zonal"). rename("comp_tuom") # Tv composite gphys = @data.gphys_open("TV")[true,true,-401..-1].lon_lotate tv_comp = gphys.composite_online(rain, 0.0007) tv_comp = tv_comp. set_att("ape_name", "Tv_(U,_-OMG)_composite"). set_att("units", "K, (m s-1, Pa s-1)" ) tv_comp = (tv_comp - tv_comp.mean(0)). set_lost_axes(""). set_lost_axes(lost_axes). add_lost_axes("(diff) from (mean) zonal"). rename("comp_tvuv") # Q composite gphys = @data.gphys_open("Q")[true,true,-401..-1].lon_lotate q_comp = gphys.composite_online(rain, 0.0007) q_comp = q_comp. set_att("ape_name", "Q_(U,_-OMG)_composite"). set_att("units", "kg/kg, (m s-1, Pa s-1)" ) q_comp = (q_comp - q_comp.mean(0)). set_lost_axes(""). set_lost_axes(lost_axes). add_lost_axes("(diff) from (mean) zonal"). rename("comp_quv") # U composite gphys = @data.gphys_open("U")[true,true,-401..-1].lon_lotate u_comp = gphys.composite_online(rain, 0.0007) u_comp = u_comp - u_comp.mean(0) # OMG composite gphys = @data.gphys_open("OMG")[true,true,-401..-1].lon_lotate om_comp = gphys.composite_online(rain, 0.0007) om_comp = om_comp - om_comp.mean(0) # ベクトルの間引 u_data = NArray.sfloat( u_comp.axis(0).pos.val.size , u_comp.axis(1).pos.val.size ).fill!(0.0) om_data = NArray.sfloat( om_comp.axis(0).pos.val.size , om_comp.axis(1).pos.val.size ).fill!(0.0) num_l = u_comp.axis(0).pos.val.size/120 num_z = u_comp.axis(1).pos.val.size/48 # num_l = u_comp.axis(0).pos.val.size/120*3 # num_z = u_comp.axis(1).pos.val.size/48 u_comp.axis(0).pos.val.size.times{ |num| if (num % num_l) == 0 u_data[num,true] = u_comp.data.val[num,true] om_data[num,true] = om_comp.data.val[num,true] elsif u_data[num,true] = 0.0 om_data[num,true] = 0.0 end } num_z =1 if num_z < 1 u_comp.axis(1).pos.val.size.times{ |num| if (num % num_z) == 0 u_data[true,num] = u_data[true,num] om_data[true,num] = om_data[true,num] elsif u_data[true,num] = 0.0 om_data[true,num] = 0.0 end } # 軸の long_name 直し (pressure level => sigma) z = t_comp.grid_copy.axis(1).pos. set_att("long_name","sigma"). set_att("units","1") z = Axis.new().set_pos(z) grid = t_comp.grid_copy.change_axis(1, z) # 配列数取得 dim = t_comp.axis(0).pos.val.size #=begin # 描画 (T) t_comp = GPhys.new( grid, t_comp.data ) # t_comp = GPhys.new( grid, t_comp.data ).rename("comp_tuom_mono") mkfig_plot(t_comp[(dim/3)..(dim*2/3),true], u_data[(dim/3)..(dim*2/3),true], -om_data[(dim/3)..(dim*2/3),true]) # 描画 (Tv) tv_comp = GPhys.new( grid, tv_comp.data ) mkfig_plot(tv_comp[(dim/3)..(dim*2/3),true], u_data[(dim/3)..(dim*2/3),true], -om_data[(dim/3)..(dim*2/3),true]) # 描画 (Q) q_comp = GPhys.new( grid, q_comp.data ) mkfig_plot(q_comp[(dim/3)..(dim*2/3),true], u_data[(dim/3)..(dim*2/3),true], -om_data[(dim/3)..(dim*2/3),true]) #=end # 描画 (T) t_comp = GPhys.new( grid, t_comp.data ) # mkfig_plot(t_comp, u_data, -om_data) # 描画 (Tv) tv_comp = GPhys.new( grid, tv_comp.data ) # mkfig_plot(tv_comp, u_data, -om_data) # 描画 (Q) q_comp = GPhys.new( grid, q_comp.data ) # mkfig_plot(q_comp, u_data, -om_data) end def gr_comp_point_online file_name = Array.new Dir.foreach($gr2ncfile_path) { |item| file_name.push("#{$gr2ncfile_path}/#{item}") if item =~ /.nc$/ } @data = Ape.new(file_name) # tr_tppn rain = @data.gphys_open("tr_tppn").cut(true,0,true).lon_lotate lost_axes = rain.lost_axes.to_s.sub("y=","lat=") rain = rain[true,-401..-1] grid_1 = Axis.new(). set_pos(VArray.new(NArray.sfloat(401).indgen!/4).rename("time"). put_att("units","days")) grid_0 = rain.grid_copy.axis(0) grid = Grid.new(grid_0,grid_1) rain = GPhys.new(grid, rain.data ).set_lost_axes(lost_axes) # tr_u250 u250 = @data.gphys_open("U").cut(true,0.25,true).lon_lotate lost_axis = [u250.data.get_att("lost_axis"), u250.lost_axes.to_s.sub("z=","sigma=").sub(" hPa","")] u250 = u250[true,-401..-1].rename("tr_u250") u250 = ( u250 - u250.mean(0) ) grid_1 = Axis.new(). set_pos(VArray.new(NArray.sfloat(401).indgen!/4).rename("time"). put_att("units","days")) grid_0 = u250.grid_copy.axis(0) grid = Grid.new(grid_0,grid_1) u250 = GPhys.new(grid, u250.data ).set_lost_axes(lost_axis). add_lost_axes("(diff) from (mean) zonal") # tr_u850 u850 = @data.gphys_open("U").cut(true,0.85,true).lon_lotate lost_axis = [u850.data.get_att("lost_axis"), u850.lost_axes.to_s.sub("z=","sigma=").sub(" hPa","")] u850 = u850[true,-401..-1].rename("tr_u850") u850 = ( u850 - u850.mean(0) ) grid_1 = Axis.new(). set_pos(VArray.new(NArray.sfloat(401).indgen!/4).rename("time"). put_att("units","days")) grid_0 = u850.grid_copy.axis(0) grid = Grid.new(grid_0,grid_1) u850 = GPhys.new(grid, u850.data ).set_lost_axes(lost_axis). add_lost_axes("(diff) from (mean) zonal") # t_conv composite gphys = @data.gphys_open("T")[true,true,-401..-1].lon_lotate grid_2 = Axis.new(). set_pos(VArray.new(NArray.sfloat(401).indgen!/4).rename("time"). put_att("units","days")) grid_0 = gphys.grid_copy.axis(0) grid_1 = gphys.grid_copy.axis(1) grid = Grid.new(grid_0,grid_1,grid_2) gphys = GPhys.new(grid, gphys.data ) comp = gphys.composite_online(rain, 0.0007, true) comp = comp. set_att("ape_name", "composite_point" ). set_att("units", "1" ). # rename("comp_point"). rename("comp_point_tr_tppn"). # set_lost_axes(lost_axes) set_lost_axes("") grid_1 = Axis.new(). set_pos(VArray.new(NArray.sfloat(401).indgen!/4).rename("time"). put_att("units","days")) grid_0 = comp.grid_copy.axis(0) grid = Grid.new(grid_0,grid_1) comp = GPhys.new(grid, comp.data ) # 描画 # mkfig_plot(comp) mkfig_plot_comp(rain, comp) comp = comp.rename("comp_point_tr_u250") mkfig_plot_comp(u250, comp) comp = comp.rename("comp_point_tr_u850") mkfig_plot_comp(u850, comp) end end # ----------------------------------------------------------------- class Ape def plot_comp(gphys, gphys_u=nil, gphys_v=nil) gropn(1) if $gropn == nil if gphys.class == Array plot_set(gphys[0]) else plot_set(gphys) end plot_main_comp(gphys, gphys_u) # タイトルを標準出力へ if gphys.class == Array then if gphys[0].data.get_att("ape_name") then puts "#{gphys[0].name} (#{gphys[0].data.get_att("ape_name")})" end elsif gphys.data.get_att("ape_name") then puts "#{gphys.name} (#{gphys.data.get_att("ape_name")})" else puts "#{gphys.name}" end end def mkdump_comp(gphys, gphys_u=nil, gphys_v=nil) gropn(3) if $gropn == nil plot_comp(gphys, gphys_u) file_name = "#{$file_label}_#{gphys_u.name}.gif" `rm dcl_* tmp.pnm` plot(gphys_u) `for file in *.xwd;do xwdtopnm ${file} | pnmcut 2 2 904 654 > tmp.pnm;done` `ppmtogif tmp.pnm > #{$fig_path}#{file_name}` `rm dcl_* tmp.pnm` $file_name = $pre_file end def plot_main_comp(gphys, gphys_u=nil) # attribute の long_name を消す (タイトル描画位置の都合) if gphys.data.get_att("long_name") gphys = gphys.set_att("long_name","") end GGraph.set_fig($fig_set_hash) # 描画 if gphys.axnames.size == 1 || gphys.name =~ "gt_" GGraph.line( gphys, true, $line_hash ) elsif $tone_flag == false GGraph.contour( gphys, true, $cont_hash ) elsif $tone_hash['levels'][0] == $tone_hash['levels'][1] GGraph.tone( gphys, true ) GGraph.contour( gphys, false ) unless $cont_flag == false else GGraph.tone( gphys, true , $tone_hash ) GGraph.contour( gphys, false, $cont_hash ) unless $cont_flag == false end plot_set(gphys_u) GGraph.contour( gphys_u, false, $cont_hash ) # タイトル tani = "(#{gphys.data.get_att("units")})" DCL::uxsttl("t", tani , -1.0 ) if gphys.data.get_att("ape_name") then DCL::uxmttl("t", gphys.data.get_att("ape_name").gsub("_", " "), -1.0 ) end $file_label = "#{@file}@#{gphys.name}" if $file_label == "filename" # nc ファイル名 if $gropn == 2 charsize = 0.79 - $file_label.size * 0.0115 DCL::sgtxzv(charsize,0.62,$file_label,0.8 * DCL.uzpget('rsizec1'),0,-1,1) else charsize = 0.81 - $file_label.size * 0.0095 DCL::sgtxzv(charsize,0.6,$file_label,0.8 * DCL.uzpget('rsizec1'),0,-1,1) end # トーンバー plot_set(gphys) DCL::Util::color_bar($cbar_conf) end end