今度はプログラムを使って図を描きます. プログラムを取っておくことで同じ図が何度でも描けるようになります. また,少し複雑なことをするには最初からプログラムを作ると良いでしょう.
目次
可視化の第二歩 の「いろんな断面」の最初の作図をプログラムとして電子ファイルにしてみましょう. irb で行った,
tmp = gpopen 'air.2012-01.nc/air' set_fig "itr"=>32 # 正距方位図法 set_map "coast_world"=>true tone_and_contour tmp[true,true,3,9]
という入力を再現するには,それらが入ったテキストファイルを作ります.
ここではファイル名を sample_tone_1.rb とします
(Rubyプログラムは .rb で終わるようにします).
その際少し追加が必要です.
スタートアップファイル irbrc_ggraph.rb
は冒頭で陽に require
しましょう
*1.
また,描画終了処理のため,DCL.grcls
を呼びます:
ソースコード sample_tone_1.rb:
require File.expand_path("~/irbrc_ggraph.rb") # irb用スタートアップファイル tmp = gpopen 'air.2012-01.nc/air' set_fig "itr"=>32 # 正距方位図法 set_map "coast_world"=>true tone_and_contour tmp[true,true,3,9] DCL.grcls # 描画終了処理.末尾に必ずよぶ
実行:
$ ruby sample_tone_1.rb
結果: (クリックでフルサイズ表示)
上記のように irb 用スタートアップファイルを用いず, GGraph の機能だけを使って同じ図を描いてみましょう. プログラムの中身は次のようにしました.
ソースコード sample_tone_1b.rb:
require "numru/ggraph" include NumRu # DCLの準備 (中身は好きに編集する) def prep_dcl(iws) # iws : DCL出力デバイス.1:画面, 2:FILE(default:PDF) DCL.sgscmn(10) # カラーマップ番号を10番に. DCL.gropn(iws) DCL.sgpset('isub', 96) # 下付き添字を表す制御文字を '_' から '`' に # (アンダースコアはよく使われるので,エラー防止のため) DCL.glpset('lmiss',true) # DCLの欠損値処理を on に. end tmp = GPhys::IO.open 'air.2012-01.nc', 'air' prep_dcl(1) GGraph.set_fig "itr"=>32 # 正距方位図法 GGraph.set_map "coast_world"=>true GGraph.tone_and_contour tmp[true,true,3,9] GGraph.color_bar DCL.grcls # 描画終了処理.
前の例からの変更のポイントは次です.
prep_dcl
として,それを呼ぶようにした
*2.gpopen
は GPhys のもともとのメソッド GPhys::IO.open
に変えた
(GPhys::IOのマニュアル 参照).include GGraph
はやめ,
GGraph のメソッドが陽に見えるようにした(それで set_fig
等の GGraph メソッドの前に陽に GGraph.
がつくようになった.これは趣味の問題.しなくてもいい).プログラムなら複雑な描画も実現できます (最初は簡単な描画を実装し,実行して調整するというのを繰り返して,徐々に充実させていくとよいでしょう). また,入力パラメタを設けることで用途を広げることもできます.
以下のプログラムは,2012年の北半球の総観場をアニメで表示します. ソース冒頭にコメントとして書いてあるように, コマンドライン引数としてオプションが指定できます. つまり,
$ ruby sample_synop.rb
とオプション無しで実行することも,
$ ruby sample_synop.rb 2 0 21 30
などと,オプションをつけて実行することもできます.
ソースコード sample_synop.rb:
# coding: cp932 # = 2012年1月の北半球の総観場(NCEP再解析) # 使用法: # ruby sample_synop.rb [オプション] # オプション(スペースで区切って並べる): # 1番目: DCLの装置番号 (default: 1) # 2番目: 描画間隔(秒) (default: 0.5) 正ならアニメ,0以下ならクリックで次を表示 # 3番目: 最初の日付(1月何日か) (default: 1) # 4番目: 最後の日付(1月何日か) (default: 31) # 使用例: # ruby sample_synop.rb # ruby sample_synop.rb 4 1 # ruby sample_synop.rb 4 0 21 30 require "numru/ggraph" include NumRu #< コマンドラインオプション > iws = ( ARGV[0] || 1 ).to_i # 第1引数は装置番号 (1,2,or 4) tsleep = [ ( ARGV[1] || 0.5 ).to_f, 0.0 ].max # 第2引数は描画間隔(秒) wait = ( tsleep <= 0.0 ) #->true/false; 0以下(true)ならマウスクリックを待つ fday = ( ARGV[2] || 1 ).to_i # 第3引数は最初の日付(1月何日か) lday = ( ARGV[3] || 31 ).to_i # 第4引数は最後の日付(1月何日か) #< 使用データ > tmp = GPhys::IO.open 'air.2012-01.nc', 'air' hgt = GPhys::IO.open 'hgt.2012-01.nc', 'hgt' #< DCL設定(メソッド) > def prep_dcl(iws=1,wait=false) DCL.sgscmn(4) # カラーマップ番号指定 DCL.swpset("iwidth",900) # 画面の幅 DCL.swpset("iheight",450) # 画面の高さ DCL.swpset("lwait",wait) # 次の描画の前にマウスクリックを待つ DCL.swpset("lalt",true) # 裏で描画(パラパラアニメ用) DCL.gropn(iws) # DCL描画装置初期化 (iws=1 or 2) DCL.sldiv("y",2,1) # 画面分割 (描画順("y"oko/"t"ate),数:横,数:縦) DCL.sgpset('isub', 96) # 下付き添字を表す制御文字を '_' から '`' に # (アンダースコアはよく使われるので,エラー防止のため) DCL.glpset('lmiss',true) # DCLの欠損値処理を on に. end #< 描画 > prep_dcl(iws,wait) for d in fday..lday time = Date.new(2012,1,d) #<fig 1> GGraph.set_fig "itr"=>32,"viewport"=>[0.05,0.81,0.12,0.88] # 32:正距方位図法 GGraph.set_map "coast_world"=>true GGraph.tone tmp.cut("level"=>850,"time"=>time),true,"title"=>"T & Z", "min"=>230,"max"=>300,"int"=>5 GGraph.contour hgt.cut("level"=>850,"time"=>time),false,"int"=>50 GGraph.color_bar #<fig 2> GGraph.set_fig "itr"=>1,"viewport"=>[0.16,0.81,0.2,0.85] GGraph.set_axes "xmaplabel"=>"lon", "xlabelint"=>90, "xtickint"=>30 GGraph.tone tmp.cut("lat"=>45,"time"=>time,"level"=>1000..100).eddy(0),true, "title"=>"T' & Z'","min"=>-24,"max"=>24,"int"=>4 GGraph.contour hgt.cut("lat"=>45,"time"=>time,"level"=>1000..100).eddy(0), false,"int"=>100 GGraph.color_bar sleep tsleep end DCL.grcls # 描画終了処理.
結果の一枚: (クリックでフルサイズ表示)
このプログラムは,2012年1月の日々の気温(色)・等圧面高度(コンター)の 水平断面(850hPa)と鉛直断面(45N)をアニメで表示します (鉛直断面は経度平均からの偏差.オプションで期間や描画間隔等が指定できます). 見どころは,移動性高低気圧発達時(左図でみる)に 高度の偏差の鉛直断面(右図)が西に傾くかや, (地衡流からわかる)南風/北風に対応して高温/低温偏差になるか,などです.
プログラムについて,本チュートリアルでこれまで出てこなかったところを解説します.
Ruby ではコマンドライン引数は,ARGV
という組み込みの配列に文字列として入ります(スペースで区切って順番に).
iws = ( ARGV[0] || 1 ).to_i
では最初の要素 ARGV[0]
があればそれを,なければ 1 を,整数化して変数 iws
に格納します
*3.
コマンドライン引数を使わず, (初歩的な計算機演習でよくあるように) プログラムを実行するとパラメタ入力を求めるようにもできますが, 引数にするほうが繰り返しの実行が楽に出来ます. なお,本格的なコマンドラインオプションをサポートするには, Ruby の組み込みライブラリ optparse を使うとよいでしょう(optparse は ARGV をコマンドラインオプションとして解釈するためのライブラリです).
DCLの設定では,DCL.sldiv("y",2,1)
で画面分割します
(第一引数は "y" なら横→縦並び,"t"なら,縦→横並び).
横幅の広い画面を確保するため,DCL.swpset
で画面の縦横のピクセル数を指定します(出力先がPostScriptの
場合もここで指定した縦横サイズの比が確保されます).
GPhys の cut
メソッドでは,座標軸が時刻の場合
(単位が "days since 2001-01-01" などとなっていてカレンダーと対応がつく場合),
直接の座標値のほか,Date または DateTime
でも指定できます(time = Date.new(2012,1,d)
, tmp.cut("level"=>850,"time"=>time)
のところ).
GPhys の eddy
メソッド
(GGraph.tone tmp.cut("lat"=>45,"time"=>time,"level"=>1000..100).eddy(0)
のところ)では,平均からのずれ(偏差)を計算します.eddy(0)
で,第1次元(ここでは経度)に沿った平均(いわゆる帯状平均)からのずれになります.
GGraph.set_axes
では "xmaplabel"
等により,
経度にふさわしい座標軸を描きます.(その他のGGraphメソッドは既出.
GGraphのマニュアル (日本語訳) も参照してください.)
GGraph.set_fig "map_axis"=>[180,-90,0]
で行います."map_axis"
は地図投影の「接点」を指定します.
その3パラメタは順に,接点の経度, 緯度,経線の回転角
です(DCL.umpsnt の UXC, UYC, ROT です:
DCLのGRPH2マニュアル参照).[次へ進む]
*1冒頭に require
行を入れない場合,
ruby -r ~/irbrc_ggraph.rb sample_tone_1.rb と実行します.
*2 メソッドにせず,そのままベタに書いてもいいのですが,
まとまりがわかりやすいようそうしました
(読者が使いまわすことを念頭に,このプログラムには不必要なものも含まれてます).
irbrc_ggraph.rb
の中をみて,他にも好きなように取り入れるといいでしょう.
自分用のオリジナルな prep_dcl
を独立したファイルにして require
するのも良いです.
*3 || は "or" をとる演算子で,その左側を実行したら偽
(false か nil) の場合,右側を実行してその結果を返します.
1はもともと整数なので整数化する必要はありませんが,||
を活かすための方策です.