このwiki内では、グリッドファイル(GrADS 4byteバイナリ、特に断らなければdirect書き込み)を〜.bin、コントロールファイルを〜.ctl、変数名をvarと表記します。用途に応じた、GrADSコマンドにより一連の設定を行う例を掲載しています。
東北大学大学院理学研究科 流体地球物理学講座 GrADSスクリプトのTips
GrADS スクリプトリファレンス
北海道大学 海洋気候物理学研究室 Grads Tips
GrADS Scripting Language
http://akyura.sakura.ne.jp/study/GrADS/Script/oper...
| 論理和(OR) & 論理積(AND) ! 否定(NOT) - 負 % 連結 = 等しい != 等しくない > より大きい >= より大きいか等しい < より小さい <= より小さいか等しい + 加算 - 減算 * 乗算 / 除算
http://akyura.sakura.ne.jp/study/GrADS/Script/stin...
端末(標準出力)に情報を書く:
say 式表記
入力を要求する文字列を書く:
prompt 式表記
「prompt」コマンドはリターンをつけない以外は「say」コマンドと同様の働きをする。
標準入力から入力された文字列や値を読む:
pull 変数
if (式表記) ..... else ..... endif
else if文は存在しない。どうしても必要な場合は以下で代用可能。
if (式表記1) ..... else if (式表記2) ..... else ..... endif endif
http://akyura.sakura.ne.jp/study/GrADS/Script/intr...
一つの文字列から一つの単語を得る関数:
res = subwrd(文字列,整数)
戻値は文字列中の「整数」番目の単語(スペースで区切られる)。もし文字列が短すぎれば、結果は「NULL」文字列。
文字列が含まれる複数行から一つの行を得る:
res = sublin(文字列(複数行を含む,整数)
結果は文字列から「整数」番目の行(リターンで区切られる)。もし文字列の行が少なすぎれば、「NULL」文字列を返す。
文字列の一部分を得る:
res = substr(文字列,start,length)
文字列の「start」番目から長さ「length」の文字列部分を得る。もし文字列が短すぎれば結果が短いか「NULL」文字になる。「start」と「length」は整数でなければならない。
4番目の単語を取得
'd hogehoge' say result wrd=subwrd(result, 4)
2列目の4番目の単語を取得
'd hogehoge' say result lin=sublin(result, 2) wrd=subwrd(lin, 4)
2列目の4番目の単語の1文字目から4文字目までを取得
'd hogehoge' say result lin=sublin(result, 2) wrd=subwrd(lin, 4) str=substr(wrd, 1, 4)
'set parea 0.8 10.8 1.0 8.3' 'set xlopts 1 6 0.22' 'set ylopts 1 6 0.22' 'set grid on 3 1' 'xcbar 0.8 10.8 0.4 0.8 -fwidth 0.18 -fheight 0.19 -edge triangle -line on'
'print hoge.eps' '!convert -rotate 90 -density 600 -resize 600 +antialias hoge.eps hoge.png'
'set line 0 1 6' 'draw recf 1.0 7.0 4.0 7.5' 'set line 1 1 6' 'draw rec 1.0 7.0 4.0 7.5'
'set gxout contour' 'set clopts 1 5 0.22' 'set cthick 5' 'set ccolor 1' 'set cstyle 1' 'set clab on' 'set grads off'
'set gxout vector' 'set clopts 1 5 0.17' 'set cthick 5' 'set ccolor 1' 'set arrscl 0.3 1.0' 'set arrowhead 0.12' 'set grads off' 'd skip(u, 4, 2);v'
そのまま重ね描きしてもよいが、白い色で太く描いた上に、黒で細く描くと、要素ごとに浮き上がって見えるので見やすい。
(釜江 2015 天気 図7より)
'set gxout contour' 'set clab off' 'set cthick 10' 'set ccolor 0' 'd var' 'set cthick 4' 'set ccolor 1' 'd var'
'set gxout line' 'set vrange -10 10' 'set cthick 3' 'set ccolor 1' 'set cstyle 1' 'set cmark 0' 'set tlsupp year' 'set grads off' 'd const(var(t=1), 0.0)' 'set digsiz 0.15' 'set cthick 5' 'set ccolor 1' 'set cstyle 1' 'set cmark 2' 'set tlsupp year'
markplot, trackplotとdraw stringを使うと、凡例を好きな形式で描ける。
'set parea 0.6 4.8 1.0 4.0' 'set xlopts 1 6 0.18' 'set ylopts 1 6 0.18' 'set grid on 3 15' 'set x 1' 'set y 1' 'set t 1 29' 'set gxout bar' 'set ylab on' 'set ylpos 0 r' 'set ylint 100' 'set vrange 0 500' 'set ccolor 15' 'set barbase bottom' 'set bargap 0' 'set baropts filled' 'set grads off' 'd var1' 'set gxout line' 'set vrange -2.5 2.5' 'set ylab on' 'set ylpos 0 l' 'set ylint 1.0' 'set grid off' 'zeroline' 'set digsiz 0.15' 'set cthick 5' 'set ccolor 1' 'set cstyle 1' 'set cmark 2' 'd var2'
GrADSスクリプトライブラリ→zeroline.gs
平均値hogemeanを折れ線で、平均値に分散stddevを±した範囲minus〜plusを四角で、
最小hogemin〜最大hogemaxの範囲をエラーバーで表示
'set gxout errbar' 'set bargap 50' 'set baropts outline' 'set ccolor 4' 'd hogemin;hogemax' 'set gxout bar' 'set ccolor 3' plus = '(hogemean+stddev)' minus = '(hogemean-stddev)' 'd 'minus';'plus 'set gxout line' 'set cmark 0' 'set cthick 6' 'set digsiz 0.05' 'set ccolor 2' 'd hogemean'
OpenGrADS Cookbooks
統計的に有意な領域や、一定の条件を満たす領域を示す場合、シェードの上に何らかのマスキングを重ねる場合が多い。
斜め線のハッチはhatch.gsを使えば可能。 GrADSスクリプトライブラリ → hatch
GrADS2ではハッチや点描の重ね描きが可能。 GrADS2まとめ → GrADS v2.1以降
一つのグリッドに一つのドットを描きたい場合は以下の方法で描くことが可能。
'd var' ## シェード等で物理量をプロット 'set gxout grid' 'set gridln off' 'set cthick 8' 'set digsiz 0.02' 'd maskout(const(var,0), mask)' ## 任意のマスキング処理
実際には小さい"ゼロ"を太字にすることで、点に見せかけている。
http://grads.wikidot.com/plot-layout
'd hogehoge' say result lin=sublin(result, 2) wrd=subwrd(lin, 4) a=substr(wrd, 1, 8) 'draw string 1.0 7.0 'a
http://akyura.sakura.ne.jp/study/GrADS/Script/intr...
#プロット 'set gxout scatter' 'set cmark 3' 'set ccolor 1' 'set digsiz 0.15' 'set t 1 30' 'set vrange -3 3' 'set xlint 1.0' 'set vrange2 -3.5 3.5' 'set ylint 1.0' 'set grads off' 'd var;var2' #相関係数を表示 'set line 0 1 5' 'draw recf 8 7 10 7.5' 'set line 1 1 6' 'draw rec 8 7 10 7.5' 'd tcorr(var, var2, t=1, t=30)' say result lin=sublin(result, 2) wrd=subwrd(lin, 4) cor=substr(wrd, 1, 6) 'set strsiz 0.2 0.2' 'set string 1 tr 4.0 0' 'draw string 8 7.8 Cor='cor #回帰係数を表示 'd tregr(var, var2, t=1, t=30)' say result lin=sublin(result, 2) wrd=subwrd(lin, 4) reg=substr(wrd, 1, 6) 'set strsiz 0.2 0.2' 'set string 1 tr 4.0 0' 'draw string 8 7.4 Reg='reg #回帰直線を描画 'regline var var2 -c 1 -l 1 -t 5' #タイトル表示 'draw title TITLE' 'draw xlab x axis title [unit]' 'draw ylab y axis title [unit]'
GrADSスクリプトライブラリ→regline.gs
欠損値の格子は、背景色(通常は黒、set display color whiteでは白)のまま、色が塗られない。
そこで予め全領域に色を塗っておくことで、欠損域を特定の色にすることができる。
例. 欠損域に(r, g, b)=(40, 20, 0)の色を塗る 'set gxout grfill' 'se clevs -1000 -998' 'set rgb 99 40 20 0' 'set ccols 0 99 0' 'd const(var, -999, -u)' ・・・ 'd var'
t=1 while(t<=??) 'set t 't 'd a1' a1=subwrd(result,4) 'a1t't'=const(a1(t=1),'a1')' 'd b1' b1=subwrd(result,4) 'b1t't'=const(b1(t=1),'b1')' 'd a2' a2=subwrd(result,4) 'a2t't'=const(a2(t=1),'a2')' 'd b2' b2=subwrd(result,4) 'b2t't'=const(b2(t=1),'b2')' t=t+1 endwhile t=1 while(t<=??) 'q w2xy a1t't' b1t't a1t.t=subwrd(result,3) b1t.t=subwrd(result,6) 'q w2xy a2t't' b2t't a2t.t=subwrd(result,3) b2t.t=subwrd(result,6) 'set line 1 1 5' 'draw line 'a1t.t' 'b1t.t' 'a2t.t' 'b2t.t t=t+1 endwhile 'draw polyf 'a1t.t' 'b1t.t' 'a2t.t' 'b2t.t' ... 'a1t.t' 'b1t.t'
各月の気候値(???ヶ月目まで) 'set t 1 12' 'clim=ave(var, t+0, t=???, 12)' 偏差 'modify clim seasonal' 'set t 1 ???' 'anom=var-clim' 気候値の季節平均 'set t 1' 'jja=ave(clim, t=6, t=8)' 'djf=( clim(t=12)+clim(t=1)+clim(t=2) )/3'
- 月平均値を年平均値に変換
t=1 while(t<=???) 'set t 't 'd ave(var,t+0,t+11)' t=t+12 endwhile
同様に6時間値→日平均値変換、etc
- 領域平均時系列
'set x 1' 'set y 1' 'set t 1 ???' 'a=aave(var,x=1,x=??,y=1,y=??)' 'd aave'
1変数の場合
'set x 1 ???' 'set y 1 ??' 'set gxout fwrite' 'set fwrite hoge.bin' t=1 while(t<=????) 'set t 't z=1 while(z<=??) 'set z 'z 'd var' z=z+1 endwhile t=t+1 endwhile 'disable fwrite'
複数の3次元、4次元変数を一つのバイナリに出力するのはさらに複雑。
(例)(360x180x1x12)の変数aaaa, bbbb, cccc, ffffと、
(360x180x10x12)の変数dddd, eeeeが含まれるバイナリを出力
'reinit' 'sdfopen inputdata.nc' 'set x 1 360' 'set y 1 180' 'set gxout fwrite' 'set fwrite -be inputdata.bin' 'set z 1' t=1 while(t<=12) 'set t 't 'd aaaa' 'd bbbb' 'd cccc' z=1 while(z<=10) 'set z 'z 'd dddd' z=z+1 endwhile z=1 while(z<=10) 'set z 'z 'd eeee' z=z+1 endwhile 'set z 1' 'd ffff' t=t+1 endwhile 'disable fwrite'
コントロールファイル
dset ^inputdata.bin title inputdata undef -9999 options big_endian xdef 360 linear 0 1 ydef 180 linear -90 1 zdef 10 levels 1 2 3 4 5 6 7 8 9 10 tdef 12 linear 01JAN0001 1mo vars 6 aaaa 1 0 ・・・ bbbb 1 0 ・・・ cccc 1 0 ・・・ dddd 10 0 ・・・ eeee 10 0 ・・・ ffff 1 0 ・・・ endvars
(例)6hourlyデータからpentad平均(5日平均)を計算
'reinit' 'open 6hourly.ctl' 'set x 1 ???' 'set y 1 ???' 'set gxout fwrite' 'set fwrite -le pentad.grd' t=1 while(t<=1460) 'set t 't z=1 while(z<=23) 'set z 'z 'climu=ave(U,t+0,t+19)' 'd climu' z=z+1 endwhile z=1 while(z<=23) 'set z 'z 'climv=ave(V,t+0,t+19)' 'd climv' z=z+1 endwhile z=1 while(z<=23) 'set z 'z 'climt=ave(T,t+0,t+19)' 'd climt' z=z+1 endwhile t=t+20 endwhile 'disable fwrite'
GrADS2で処理可能なアンサンブルデータのfwrite
複数の変数の場合は、zループの外側、tループの内側に、次以降の変数を並べてゆけばよい。
'set x 1 ???' 'set y 1 ??' 'set gxout fwrite' 'set fwrite hoge.bin' e=1 while(e<=?) 'set e 'e t=1 while(t<=????) 'set t 't z=1 while(z<=??) 'set z 'z 'd var' z=z+1 endwhile ***次の変数(必要であれば)*** t=t+1 endwhile e=e+1 endwhile 'disable fwrite'
変数varの値が正のときはその値を、
負のときは0.0をaに代入する
t=1 while(t<=3) 'd var(t='t')' # 変数resultに "Result value = ***" が格納される r = subwrd(result, 4) # 4番目の文字列が値に相当する if(r>0.0) # 値が正のとき 'a't'=r' # 値をaに代入 else # 値が負のとき 'a't'=0.0' # 0.0をaに代入 endif t=t+1 endwhile 'd a1' say result 'd a2' say result 'd a3' say result
スクリプトの動作を確認するため、途中で止める
'〜' pull hoge '…'
と書くと、〜までのスクリプトの結果が表示され、エンターキーを押すと…以降が実行される。
http://akyura.sakura.ne.jp/study/GrADS/Script/stin...
コメントをかく