気象データ解析ツールであるGrADSについての情報をまとめています。

この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文


if (式表記)
   .....	
else
   .....	
endif

else if文は存在しない。どうしても必要な場合は以下で代用可能。

if (式表記1)
   .....	
else	
 if (式表記2)
    .....	
 else
    .....	
 endif
endif

while文


while(式表記)
	.....
endwhile

whileループを継続するにはcontinueコマンド、 whileループから抜け出すには breakコマンドを使用。

文字列と文字の取得


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 display color white'
'c'

図の縦軸・横軸を見やすく


'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'

解像度の高いpngを出力


'print hoge.eps'
'!convert -rotate 90 -density 600 -resize 600 +antialias hoge.eps hoge.png'

文字列表示


'set strsiz 0.22 0.22'
'set string 1 c 4.0 0' 
'draw string 2.0 8.0 hogehoge' 

白塗りつぶし+黒枠の四角


'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

ドットを打つ(点描、stipple)


統計的に有意な領域や、一定の条件を満たす領域を示す場合、シェードの上に何らかのマスキングを重ねる場合が多い。
斜め線のハッチは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

gxout gridのオプション


'set gxout grid'
'set dignum 2'
'set digsize 0.16'

欠損域に色を塗る


欠損値の格子は、背景色(通常は黒、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'

2変数の時系列から散布図上にトラックをプロット、トラックにより囲われる多角形を塗りつぶす


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'


計算

気候値、偏差、年平均、時系列etc


各月の気候値(???ヶ月目まで)
'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'



出力

(x,y,z,t) 4次元データのfwrite


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'

(x,y,z,t,e) 5次元データの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...

ファイルをX個閉じる


nn=X
while(nn>=1)
'close 'nn
nn=nn-1
endwhile

コメントをかく


「http://」を含む投稿は禁止されています。

利用規約をご確認のうえご記入下さい

管理人/副管理人のみ編集できます