#!/usr/bin/ruby # -*- coding: utf-8 -*- # (2015-11-04 石渡) # 新学術用に試してみた. # % /GFD_Dennou_Work3/momoko/SyncRotEarthRad/script/fig_ts-uv_for-display_flow-vect_coastline.rb --xfact1 0.001 --yfact1 0.001 --uxunit 1 --uyunit 1 --loncnt 110 --range_val 225:300 --wsn 1 --range_time 1030:1031 --xintv 2 --yintv 2 # 途中でカラーマップが変えられない! # v はartificial に 2 倍した. # (2013-6-14 石渡) # 球面プロットでカラーバーを書くとベクトルが表示されない! # 矢印が太すぎるような気がする. # わからん!!! # /GFD_Dennou_Work3/momoko/SyncRotEarthRad/script/fig_ts-uv_for-display_flow-vect.rb --loncnt 45 --range_val 150:250 --wsn 2 DATA_HOME = "/GFD_Dennou_Work3/momoko/SyncRotEarthRad" require "numru/ggraph" include NumRu require "pp" wsn = 4 RPlanet = 6.371 * 10 ** 6 # [m] #Ps = 1e5 # [Pa] Grav = 9.8 # [m/s2] title = "" # SurfTempの描画範囲 # 絶対温度 #val_min = 225 #val_max = 315 # 摂氏 val_min = -58 val_max = 42 val_interval = (val_max - val_min) / 5.0 #YUNIT_vapor = 'kg m-2' #epsfn = "qvap_area_column_mass.eps" # ベクトルの描画間隔 (単位は格子点数) #xintv = 8 #yintv = 6 #xintv = 1 #yintv = 1 xintv = 4 yintv = 3 # 球面プロットと平面プロットの切り替えは # この部分と下の設定部分を変更する. ITR=30 # 球 #ITR=1 # 平面プロット lon_cont = 125 # 球のときに使用 plx=0 ply=0 #xfact1 = 0.01 #yfact1 = 0.01 xfact1 = 1.0 yfact1 = 1.0 # 凡例の単位ベクトルの長さ (デフォルト値) uxunit = 0.03 uyunit = 0.03 time_start = 731 time_end = 1095 #omegas = [0.0, 0.15, 0.5, 1.0] #omegas = [1.0] #omegas = [0.0] omega = 1.0 dir = nil # 簡易 while ARGV.size > 0 case ARGV[0] when '--dir' then ARGV.shift dir = ARGV[0] ARGV.shift when '--omega' then ARGV.shift omega = ARGV[0] ARGV.shift when '--wsn' then ARGV.shift wsn = ARGV[0] ARGV.shift when '--xintv' then ARGV.shift xintv = ARGV[0].to_i ARGV.shift when '--yintv' then ARGV.shift yintv = ARGV[0].to_i ARGV.shift when '--xfact1' then ARGV.shift xfact1 = ARGV[0] ARGV.shift when '--yfact1' then ARGV.shift yfact1 = ARGV[0] ARGV.shift when '--uxunit' then ARGV.shift uxunit = ARGV[0] ARGV.shift when '--uyunit' then ARGV.shift uyunit = ARGV[0] ARGV.shift when '--loncnt' then ARGV.shift lon_cnt = ARGV[0] ARGV.shift when '--range_time' then ARGV.shift tmp = ARGV[0].split(':') time_start = tmp[0].to_f time_end = tmp[1].to_f ARGV.shift when '--range_val' then ARGV.shift tmp = ARGV[0].split(':') val_min = tmp[0].to_f val_max = tmp[1].to_f ARGV.shift val_interval = (val_max - val_min) / 30.0 else raise "ERROR: invalid option: #{ARGV[0]}" end end if dir == nil then # dir = File.join(DATA_HOME, # "101220_diff-DryAir-Vapor_T21L16_sr_Omega" + omega.to_s + "_dt20m") dir = File.join(DATA_HOME,"ExpOcean") end # input var = 'U' path = var+'.nc' gp_u = GPhys::IO.open(path, var).cut('time'=>time_start..time_end, 'sig'=>0.9) var = 'V' path = var+'.nc' gp_v = GPhys::IO.open(path, var).cut('time'=>time_start..time_end, 'sig'=>0.9) var = 'SurfTemp' path = var+'.nc' gp_surftemp = GPhys::IO.open(path, var).cut('time'=>time_start..time_end) var = 'PRCP' path = var+'.nc' gp_prcp = GPhys::IO.open(path, var).cut('time'=>time_start..time_end) var = 'H2OLiq' path = var+'.nc' gp_cloud = GPhys::IO.open(path, var).cut('time'=>time_start..time_end).mean('sig') # 時間平均 gp_u = gp_u.mean('time') gp_v = gp_v.mean('time') gp_surftemp = gp_surftemp.mean('time') gp_prcp = gp_prcp.mean('time') #gp_surftemp = gp_surftemp - 273.0 DCL.sgscmn(07) # 描画 DCL.gropn(wsn) DCL.sgpset('lcntl', false) ; DCL.uzfact(0.7) #GGraph.set_axes('xtickint'=>10, 'xlabelint'=>30) #GGraph.set_axes('ytickint'=>10, 'ylabelint'=>30) # 球の時 (平面プロットの時は 下の 2 行をコメントアウトするだけ) vx0,vx1,vy0,vy1 = 0.1,0.9,0.1,0.9 GGraph.set_fig 'itr'=>ITR, 'viewport'=>[vx0,vx1,vy0,vy1], 'map_axis'=>[lon_cnt, plx, ply] # 海岸線を引く #set_map "coast_world"=>true #GGraph.set_map "coast_world"=>true DCL.uscset('cyspos', 'B' ) # move unit y axis DCL.udlset('LMSG', false) DCL.ugpset('XFACT1', xfact1) DCL.ugpset('YFACT1', yfact1) DCL.ugpset('UXUNIT', uxunit) # size of unit vector by dimentional value DCL.ugpset('UYUNIT', uyunit) DCL::ugrset('VXUNIT', 0.1) # size of unit vector DCL::ugrset('VYUNIT', 0.1) DCL::uxmttl('T', '', 0.0) # title (large character) DCL.ugpset('LNRMAL', false) # scale vector automatically #DCL.ugpset('LMSG', false) # no message DCL.ugpset('LUMSG', false) # no message #DCL.ugpset('INDEX', 959) DCL.ugpset('INDEX', 955) #DCL.ugpset('INDEX', 5) opt_vector = { 'xintv'=>xintv, 'yintv'=>yintv, 'factor'=>2.0, 'flow_vect'=>true, 'annotate'=>false, 'unit_vect'=>false } opt_contour = { 'min'=>val_min, 'max'=>val_max, 'interval'=>val_interval, 'title'=>title, 'annotate'=>false } opt_tone = { 'min'=>val_min, 'max'=>val_max, 'title'=>'', 'nlev'=>60, 'annotate'=>false } opt_cloud_tone = { 'min'=>1e-8, 'max'=>1.0, 'title'=>'', 'levels'=>[3.0e-6, 5.1e-6, 8.0e-6, 1.0], # 'patterns'=>[41,55], 'patterns'=>[15,25,35], 'annotate'=>false } opt_prcp_shade = { 'min'=>1, 'max'=>100, 'title'=>'', 'ltone'=>false, 'levels'=>[3,4,10], 'patterns'=>[21,43,55], 'annotate'=>false } # mm/day に変換 #prcp は [kg m-2 s-1] #rho = 1000kg/m3 #prcp / rho * 1000 * 86400 [mm/day] gp_prcp = gp_prcp * 86400.0 opt_prcp_contour = { 'color'=>true, 'clr_min'=>999, 'clr_max'=>999, 'interval'=>4.0 } # トーンのカラーマップ #DCL.sgscmn(9) # colormap: white-blue-black #DCL.sgscmn(8) # colormap: white-yellow-red-black #DCL.sgscmn(2) # colormap: black-red-yellow-white #DCL.sgscmn(6) # colormap: pastel_rainbow 薄すぎてダメ #DCL.sgscmn(62) # NCL:BlGrYeOrReVi200 #DCL.sgscmn(53) # NCL:ncview_default #DCL.sgscmn(55) # NCL:rainbow #DCL.sgscmn(54) # NCL:rainbow+white+gray #DCL.sgscmn(28) # IDV:Basic->Bright38 #DCL.sgscmn(34) # DCL.sgscmn(28) GGraph.tone( gp_surftemp, true, opt_tone) GGraph.map("grid"=>false,"coast_world"=>true) # 降水 (本当はシェードにしたかったけど. 白にできなかったのであきらめた) DCL.uepset('ltone',false) DCL.uuslni(80) ##GGraph.tone( gp_prcp, false, opt_prcp_shade) #GGraph.contour( gp_prcp, false, opt_prcp_contour) #GGraph.tone( gp_prcp, false, opt_prcp_shade) # 雲水 DCL.sgscmn(27) # NCL:rainbow+white+gray GGraph.tone( gp_cloud, false, opt_cloud_tone) # ベクトルのプロット DCL.sgscmn(77) # colormap: #DCL.sgscmn(7) # colormap: #DCL.sgscmn(77) # colormap: GGraph.vector( gp_u, gp_v, false, opt_vector) # (2013-6-27 石渡) 仕方ないのでカラーバーは最後に出力 #DCL.sgscmn(34) # colormap: black-red-yellow-white #GGraph.color_bar('landscape'=>true, 'voff'=>0.03, 'tickintv'=>0) DCL.grcls exit