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

このwiki内では、グリッドファイル(GrADS 4byteバイナリ、特に断らなければdirect書き込み)を〜.bin、コントロールファイルを〜.ctl、変数名をvarと表記します。
GrADSの基本的な使い方は、以下のサイトに系統的に整理されています。

東北大学大学院理学研究科 流体地球物理学講座 公開情報/GrADS
GrADS コマンドリファレンス
IT memo
GrADS リファレンスマニュアル
GrADS Documentation Index




重要なバグ

fwrite時のset x


fwriteの際にはxとyのsetが必須。q dimsで表示されるxの格子点数は1つ大きい場合がある。


open 〜.ctl
set x 1 144
set y 1 73
set t 1 20
set gxout fwrite
set fwrite -le hoge.bin
d var

×
open 〜.ctl
set t 1 20
set gxout fwrite
set fwrite -le hoge.bin
d var


時間平均と領域平均


時間平均と領域平均を行う際は、時間平均が先。

IT memo


a=ave(var, t=1, t=10)
d aave(a, lon=A, lon=B, lat=C, lat=D)

×
a=aave(var, lon=A, lon=B, lat=C, lat=D)
d ave(a, t=1, t=10)

画面出力


描画方法の設定 set gxout



描画方法(2次元)




d var (2 demension)
 set gxout contour (default) : 等値線
 set gxout shaded            : 陰影
 set gxout shade2            : 陰影(横線が入らないように塗りつぶし;V2.0.0以降)
 set gxout grfill            : グリッド塗りつぶし
 set gxout grid              : 数値の表示
d var;var2 (2 demension)
 set gxout vector (default)  : ベクトル
 set gxout stream            : 流線
 set gxout barb              : 矢羽

描画方法(1次元)



d var (1 demension)
 set gxout line (default)    : 折れ線
 
d var;var2 (1 demension)
 set gxout bar               : 棒、四角
 set gxout errbar            : エラーバー
 set gxout linefill          : 折れ線間の塗りつぶし

図、値の表示 d

d var

dはdisplayの略
ファイル番号を指定して表示
d var.ファイル番号

鉛直レベルを指定して表示
d var(z=2)

時間ステップを指定して表示
d var(t=2)

鉛直・時間を指定して表示
d var(z=2, t=10)

※tとzだけでなく、「固定されている次元のステップを指定して表示」することができる。
x, y, z, t 方向に次元をもつデータの場合、デフォルトではz=1, t=1に設定されているため、zとtを指定して描くことができる。
xやyを固定すると、xやyのステップを指定して描くことができる。

 set x 1
 set y 1
 set z 1 10
 set t 1 12
 d var(x=1, y=5)

鉛直断面図


経度X度における東西鉛直断面図

'set lon X' (またはset x ?)
'set y ? ??' (またはset lat ? ??)
'set z ? ??' (またはset lev ? ??)
'd var'

東西平均値の南北断面図

'set x 1'
'set y ? ??' (またはset lat ? ??)
'set z ? ??' (またはset lev ? ??)
'd ave(var, x=1, x=??)'

鉛直プロファイルの時間変化

'set x ?'
'set y ?'
'set z ? ??' (またはset lev ? ??)
'set t ? ??' (またはset time ??? ???)
'd var'


2変数を用いた表示(ベクトル、2変数の幅の表示)

風ベクトルを表示(デフォルト)
d var1;var2

矢羽を表示

 set gxout barb
 d var1;var2

流線を表示

 set gxout stream
 d var1;var2

2変数の間を棒グラフで表示(1次元)

 set x 1
 set y 1 12
 set gxout bar
 d var1;var2

2変数の間をエラーバーで表示(1次元)

 set x 1
 set y 1 12
 set gxout errbar
 d var1;var2

2変数の間を塗りつぶし(1次元)

 set x 1
 set y 1 12
 set gxout linefill
 set lfcols 1 15
 d var1;var2

3変数を用いた表示(ベクトルの色塗り分け)

風ベクトルを色で塗り分け
d var1;var2;var3

図の消去 c

c

cはclearの略
色などの設定は初期化せず、図だけ消去
c norset

東北大学大学院理学研究科 流体地球物理学講座 公開情報/GrADS

画像ファイル出力


出力方法による比較


printimよりprintのほうが、shadedよりshade2を使用したほうがきれいに出力できる。


print eps出力


IT memo

eps形式にしたい場合
print 〜.eps

90度回転させない(landscapeモード、横長画面のとき)
print -R 〜.eps



本来、eps出力は、一度grads metafileを作成し、gxepsで変換する、
 enable print tmp
 print tmp
 !gxeps -c -i tmp -o hoge.eps
 !rm -f tmp

の手順で行いますが、
enable printのコマンドを省略し、print文で直接epsファイルを出力することが可能。

http://cola.gmu.edu/grads/gadoc/gradcomdprint.html


print+convertによる高解像の画像ファイル出力

 print 〜.eps
 !convert -rotate 90 -density 600 -resize 600 〜.eps 〜.png

printコマンドでepsを作成してから、ImageMagickのconvertコマンドを使って
高解像のpng, gifを作成する。
-rotateで90度回転、-densityでピクセル数の設定、-resizeでサイズの設定。

会津大学UNIXウィキ

printim png, gif出力


png形式にしたい場合
printim 〜.png

gif形式にしたい場合
printim 〜.gif gif

mybkg.pngの背景の上に色番号0を透過させて重ねた画像を出力
printim image.png -b mybkg.png -t 0

GrADS Commands: printim

ただし、printimでは画質の低い画像しか作れません(特に数字が読みづらい!)



図の邪魔な横線を消す


  • GrADS 2.0.0から使用可能になった、set gxout shade2を使う。

GrADS2まとめ

 set gxout shade2
 d var
 print hoge.eps

上記「出力方法による比較」を参照。

http://cola.gmu.edu/grads/changelog-2.0.txt
http://oceansciencehack.blogspot.jp/2012/02/gradsp...

  • convertを用いる

png、gif等の形式であれば、epsで出力してから変換。それでも横線が目立つ場合は+antialiasオプションをつける。
convert +antialias 〜.eps 〜.png

  • gimpを用いる

gimpを起動、画像ファイルを「アンチエイリアスなし」で開き、別名保存する。

IT memo

北大佐々木さん



軸ラベルが汚い



一つの図に、複数の要素を重ね書きした場合(複数回displayコマンドを使用した場合)、
図のx軸、y軸ラベルが汚い、ラベルや図の枠が太くなるといった問題が生じる。
加えて、ファイルサイズが不必要に大きくなってしまう。

gradsではデフォルトの設定で、xlab、ylab、frame、gridがonになっており、
特に変更しない限り、軸ラベル、図の枠線、グリッドラインは、「d」が行われるたびに重ね書きされる。
例えば、set gxout lineで5本の折れ線を重ね書きしていれば、これらも5回重ね書きされている。
その状態でepsに出力すると、ラベルや線が汚く表示され、ファイルサイズも大きくなる。

※下の図の例の場合、上の図のファイルサイズは下の図の2倍以上。




これを防ぐには、最初の「d」の前で設定をoff、最後の「d」の前で設定をonにすればよい。

 'set frame off'
 'set grid off'
 'set xlab off'
 'set ylab off'
 'd var'       ##最初のdisplay文
・・・
 'set frame on'
 'set grid on'
 'set xlab on'
 'set ylab on'
 'd var'       ##最後のdisplay文


データ出力


x,y,z,t方向に幅のあるデータのfwrite


gsファイルにて

 'open 〜'
 'set x 1 〜'
 'set y 1 〜'
 'set z 1 〜'
 'set gxout fwrite'
 'set fwrite 〜.bin'
 t=1
 while(t<=〜)
 'set t 't''
 'd 〜'
 t=t+1
 endwhile
 'disable fwrite'
  • bオプションをつけて実行するとよい
grads -blc hoge.gs

東北大学大学院理学研究科 流体地球物理学講座 公開情報/GrADS

fwriteでsequentialファイルを出力


http://cola.gmu.edu/grads/gadoc/gradcomdqfwrite.ht...
set fwrite <-be or -le> <-sq or -st> <-ap or -cl> fname

big endian direct書き込みの場合
set fwrite -be hoge.bin

little endian direct書き込みの場合
set fwrite -le hoge.bin

big endian sequential書き込みの場合
set fwrite -be -sq hoge.bin

出力したファイルを読むコントロールファイルは
options big_endian
options little_endian
options big_endian sequential

など

計算、関数


平均 ave, mean


次元(x, y, z, t)または軸(lon, lat, lev, time)についての平均
d ave(var, x=1, x=144)

12ごとの平均(monthly気候値)

 set t 1 12
 d ave(var, t+0, t=120, 12) 

経度0〜360°平均を取る場合は、0°の値と360°の値が重複して計算されないよう、-bをつける。

 set lon 0
 set lat -90 90
 a=ave(var, lon=0, lon=360, -b)
 d a

重み(格子の大きさの違い)を考慮しない平均
d mean(var, x=1, x=144)

領域平均 aave


GrADS コマンドリファレンス

経度A〜B°、緯度C〜D°領域の平均(重みも考慮)

 a=aave(var, lon=A, lon=B, lat=C, lat=D)
 set x 1
 d a

全球平均

 a=aave(var, global)
 d a

重み(格子の大きさの違い)を考慮しない平均
d amean(var, x=1, x=144, y=1, y=73)

合計

 sum(var, x=1, x=144)                   ある次元の合計
 sumg(var, x=1, x=144)                   ある次元の重みを考慮しない合計
 asum(var, x=1, x=144, y=1, y=73)        ある面の合計
 asumg(var, x=1, x=144, y=1, y=73)        ある面の重みを考慮しない合計

積分 gint


varのaからbまでの積分値
gint(var, a, b)

鉛直積分 vint


気圧面がhPa表記のとき。最下層から10hPaまでの鉛直積分
d vint(lev(z=1), var, 10)

Pa表記のとき
d vint(lev(z=1), var, 1000)/100

e.g. 比湿[kg/kg]をvintすると可降水量[mm=kg/m2]になる。vintは各層の値にdp/g[Pa s2/m = kg/m2]をかけているため。ただし対応するctl(nc)の気圧面がPa表記の場合は、hPa表記に直すか、上記のように100で割る必要がある。100000, 1000はあくまでhPaとして受け取られているので、vintをすると100倍のhPa s2/m をかけてしまうため。kg/m2 にするには100で割る必要がある。

!!まちがい!!
d vint(var(z=1)-var(z=1)+1000, var, 10)

のやり方では、z=1におけるvarが欠損値の場合、開始気圧も欠損値扱いになってしまい、計算ができない。

周囲9点重みつき平均 smth9

d smth9(var)

hcurlやhdivgはsmth9をしたほうが、隣同士のグリッドの差を平滑化できる
d smth9(hdivg(ua, va))

増田さんのページ

中央差分 cdiff


GrADS Documentation Index

x(またはy, z, t)軸方向の中央差分
cdiff(var, x)

※x=1とx=nxでも計算を行いたい場合は、x軸に1格子ずらしたctlファイルを二つ作り、
差分をとるという手がある。

渦度 hcurl


GrADS Documentation Index

東西風uaと南北風vaから渦度を計算
hcurl(ua, va)

鉛直断面図

 set lon 0 360
 set lat -90 90
 set lev 1000 100
 a=hcurl(ua, va)
 set lon 0
 d a

収束・発散 hdivg


GrADS Documentation Index

東西風uaと南北風vaから収束発散を計算
hdivg(ua, va)

相関係数 tcorr, scorr


GrADS Documentation Index

時間方向

 set t 1
 d tcorr(var1, var2, t=始め, t=終わり)

時間のみの関数var1と東西・南北・時間の関数var2の変動の相関係数分布

 set x 1
 set y 1
 set t 1 100
 a=var1
 set lon 0 360
 set lat -90 90
 set t 1 100
 b=var2
 set t 1
 d tcorr(a, b, t=1, t=100)  ※

※ 時間軸の情報は一致していないといけない(同時相関)
言い替えれば時間軸情報をずらせばlag相関が可能

時間のみの関数var1のctlファイルと、東西・南北・時間の関数var2のctlファイルを開いて、その変動の相関係数をfwriteするときは、少々面倒な手順となる。

 open var1.ctl
 set x 1
 set y 1
 set t 1 100
 a=var1
 close 1      ※1
 open var2.ctl
 set lon 0 360
 set lat -90 90
 set t 1 100
  b=var2
 set t 1
 set gxout fwrite
 set fwrite hoge.bin
 set x 1 144    ※2
 set y 1 73
 d tcorr(a, b, t=1, t=100)
 disable fwrite

※1 格子情報の異なるctlファイルなので、a=で定義したのち、ファイルを閉じないと格子情報の衝突が起こる。
※2 fwriteの際はdの前にxとyの再設定が必要になる。

空間方向
d scorr(var1, var2, lon=a, lon=b, lat=c, lat=d)

IT memo

回帰係数 tregr, sregr


GrADS Documentation Index

時間方向

 set t 1
 d tregr(var1, var2, t=始め, t=終わり)

空間方向
d sregr(var1, var2, lon=a, lon=b, lat=c, lat=d)

該当する最初の高度 fndlvl


GrADS Documentation Index

varが初めてvar2と等しくなる高度をlev=1000〜100の間で探す
d fndlvl(var, var2, lev=1000, lev=100)

varが初めて100になる高度をlev=1000〜100の間で探す
d fndlvl(var, const(var, 100), lev=1000, lev=100)

20 degree isotherm depth(水温20℃等温度深度)
d20=fndlvl(temp,const(temp,20),lev=4000,lev=5)

2乗和の平方根 mag


GrADS Documentation Index

var1とvar2の2乗和の平方根
mag(var1, var2)

定数、欠損値の強制設定 const


東北大学大学院理学研究科 流体地球物理学講座 公開情報/GrADS

定数0.1としてのvar
const(var, 0.1)

欠損値を0としたvar
const(var, 0, -u)

欠損値を含め全ての値を0としたvar
const(var, 0, -a)

絶対値 abs


varの絶対値
abs(var)

eのx乗 exp


GrADS Documentation Index

eのvar乗
exp(var)

累乗 pow


GrADS Documentation Index

varのa乗 (aは変数でも定数でもよい)
pow(var, a)

平方根 sqrt


GrADS Documentation Index

varの平方根
sqrt(var)

対数 log log10


 log(var)  	自然対数
 log10(var)  	常用対数

三角関数 cos(var) sin(var) tan(var) acos(var) asin(var) atan2(a, b)


 〜(var)    var(ラジアン)の三角関数
 〜(var*3.1416/180)    var(°)の三角関数
 atan2(a, b)    a/bのarctan(ラジアン)

緯度のcos
cos(lat*3.1416/180)

等経度で緯度差dlatの間の距離 [m]
2*6.37*1e+6*3.1416*dlat/360

緯度latにおける等緯度面の円周 [m]
2*6.37*1e+6*cos(lat/180*3.1416)*3.1416

等緯度で経度差dlonの間の距離(緯度latにおける) [m]
2*6.37*1e+6*cos(lat/180*3.1416)*3.1416*dlon/360



同心円状の値を作成

cosをlat, lonの変数として用いて円状の値を作成。ピークの値、中心位置、幅を設定可能。

 'a=cos((lon-90)*3/180*3.1416)+cos(lat*3/180*3.1416)'
 'set gxout shaded'
 'color 0.3 1.8 0.3 -kind white->aqua->lime->yellow->red'
 'd 1*const(maskout(maskout(maskout(maskout(maskout(a,a)
 ,120-lon),lon-60),30+lat),30-lat),0.0,-u)'
 'cbar'




最大値 max, 最小値 min


GrADS Documentation Index

varの最大値、最小値

 max(var, 始点, 終点)
 min(var, 始点, 終点)

最大値、最小値をとる位置 maxloc, minloc


GrADS Documentation Index

t=1でvarが最大値をとるx

 set lat -90 90
 set t 1
 d maxloc(var, lon=0, lon=360)

lat=0でvarが最大値をとるtを2間隔で

 set lon 0 360
 set lat 0
 d maxloc(var, t=1, t=12, 2)

z=1からz=18の間でvarが最大となるzを出力するgradsスクリプト
'reinit'
'open hoge.ctl'
'set z 1'
'set t 1'

y=1
while(y<=64)
'set y 'y
x=1
while(x<=128)
'set x 'x

'd maxloc(hoge(x='x',y='y'),z=1,z=18)'
*say result
lin=sublin(result,20)
*say lin
wrd=subwrd(lin,4)
*say x' 'y' 'wrd
'x'x'y'y'=const(hoge(z=18),'wrd')'

x=x+1
endwhile
say y
y=y+1
endwhile

'set gxout fwrite'
'set fwrite -be hogemaxminloc.bin'
y=1
while(y<=64)
'set y 'y
x=1
while(x<=128)
'set x 'x
'd x'x'y'y''
x=x+1
endwhile
y=y+1
endwhile
'disable fwrite'

該当するとき描かない maskout


GrADS Documentation Index

mask>=0のときだけvarを描く
d maskout(var, mask)

西風だけ描く
d maskout(ua, ua)

絶対値が3以上のときだけ描く
d maskout(var, abs(var)-3.0)

地形をマスク(地上気圧がその面の気圧よりも低いとき、その気圧面は欠損とする)。ただしlevとpsは単位を合わせる。
d maskout(var, ps-lev)

maskoutの領域に色をつける → GrADSスクリプトでよく使う文まとめ→欠損域に色を塗る
複数種類のシェードを同時に描く

東経60〜120°、緯度-20〜40°の範囲のみ、色を変える(特に意味はない)

 d ave(precip, t=6, t=8)
 a=maskout(ave(precip, t=6, t=8), (-abs(lon-120)+60))
 b=maskout(a, -lat+40)
 d maskout(b, lat+20)



shadedは値が存在するところには"white"なり何かしらの色を塗ってしまうため、
シェードにシェードを重ねることはできない。
ところが欠損域には色が"塗られない"ため、maskoutで特定の領域のみシェードを描き、
複数のシェードを一枚に重ねることもできる。

なお、printimの-tオプション、-bオプション、-fオプションを使えば、特定の色を透過させ、複数の図を重ねることもできる。
http://cola.gmu.edu/grads/gadoc/gradcomdprintim.ht...

特定の値、特定の領域のみを扱う計算


maskoutとconstを組み合わせることで、
「ある変数(値、緯度など)を閾値として計算に含める・含めない・0として扱う」
などが可能。
扱い欠損域でない欠損域である
d varundef以外undefで定義された値
d maskout(var, var-2.0)undef以外の2.0以上undef, 2.0以下
d const(var, 0.0)undef以外全て0.0undef
d const(var, 0.0, -u)undef以外、undefは0.0なし
d const(var, 0.0, -a)undef含めて全て0.0なし
d const(maskout(var, var-2.0), 0.0)undef以外の2.0以上は0.0undefと2.0以下
d const(maskout(var, var-2.0), 0.0, -u)undef以外の2.0以上はそのまま、undefと2.0以下は0.0なし
d const(maskout(var, var-2.0), 0.0, -a)undef含めて全て0.0なし
d maskout(maskout(var, lon-120), 150-lon)undef以外の120<lon<150の範囲undef、lon<120、lon>150

任意の長方形
rec1=maskout(maskout(maskout(maskout(const(lev,0.0),145-lon),lon-135),lat-30),lat-50)
→135<lon<145、30<lat<50の長方形

任意の長方形の組み合わせ
 rec1=・・・   *長方形内は0、外は欠損値のデータ
 rec2=・・・
 rec3=・・・
 
 aa=const(rec1,-1.0,-u)+const(rec2,-1.0,-u)+const(rec3,-1.0,-u)
 area1=const(maskout(aa,aa+0.1),0.0)

→任意の領域内は0、外は欠損値のデータ

 open var.ctl   *任意のデータをopen
 d var+area1   *任意のデータを、領域に該当するグリッドはその値、外は欠損値のデータに変換



地図


GrADS コマンドリファレンス

地図を描くためには、あらかじめ環境変数GADDIR(詳細)で指定したディレクトリに、
ftp://cola.gmu.edu/grads/data.tar.Z
からダウンロード、解凍した地図情報ファイル(lowres, mres, hires)を置いておく必要がある。


地図を描く/描かない

set mpdraw on/off

地図の解像度{低, 中, 高}


デフォルトは低
set mpdset {lowres mres hires}

地図の色、線種、太さを指定

set map 1 1 3

海岸線を太くして強調





d var
printim hoge.png



set grid off
set frame off
set map 1 1 8
d var
set gxout contour
set clab off
set cthick 8
set clevs 0
set ccolor 1
d lat
print hoge.eps

地図投影法の指定


正距円筒図法/直交座標/ポーラーステレオ北極/南極/ランベルト/モルワイデ/正投影/ロビンソン/なし
set mproj latlon/scaled/nps/sps/lambert/mollweide/orthogr/robinson/off

GrADS Tutorial

地図投影法 投影カタログ





国境、州境線 {描く, 描かない}


デフォルトは描く、但し低解像度では描かれない
set poli {on off}

日本の県境を描く

新しい地図


mpdsetをmres, hiresに設定し、poliをonに設定した時に描かれる国境線はかなり古い(トルクメニスタンが無いなど)。
2001年当時の地図情報ファイル:
ftp://cola.gmu.edu/grads/boundaries/newmap

詳細な海岸線情報ファイル:
ftp://cola.gmu.edu/grads/boundaries/worldmap

上記ファイルをGADDIR(詳細)で指定したディレクトリに置き、

 set mpdset newmap
 set poli on

など
地図情報ファイル海岸線解像度国境線
lowres(デフォルト)低い-低い
mres
hires高い高い高い
newmapやや高いやや高い、新しい-
worldmapかなり高い--

中部日本


中東



白地図


スライド作成等で、白地図が必要なとき用に。トリミングしてお使い下さい。

フルサイズ pngファイル1 pngファイル2 pngファイル3(右クリックで保存)




        




文字


タイトルを書く

draw title hogehoge

小さい文字でタイトルを書く
draw title \hogehoge\

上:Title 下:\Title\





文字列を書く


GrADS コマンドリファレンス

set string 色 座標の位置 太さ 回転
set strsiz 幅 高さ
draw string x y 文字列

例. 大きなタイトル表示

set parea 1.0 10.0 1.0 7.0
d var
set string 1 bc 6.0 0
set strsiz 0.5 0.5
draw string 5.5 7.0 hogehoge

フォントの変更

上付き、下付き、特殊文字


東北大学大学院理学研究科 流体地球物理学講座 公開情報/GrADS
GrADS pointer

上付き
`a〜`n

下付き
`b〜`n

特殊文字
`番号〜`0

特殊文字の対応表の表示
run font.gs 番号



`3.`0C

Wm^-2
Wm`a-2

解説
run font.gs 3

とすると、右上に、「.」が「°」に対応しているので、フォント番号3の「.」は「°」であることがわかる。
「23℃」を表示したい場合は、「23」「フォント番号3」「.」「フォント番号0」「C」つまり
23`3.`0C

とする。

特殊文字の追加


フォント一覧表



フォント番号6以降は、font6.datなどを自作すればフォントを追加可能。

http://cola.gmu.edu/grads/gadoc/gradcomdsetfont.ht...

キリル文字、ラテン文字の追加
  • インストール

1. font6.datfont7.datfont8.datをダウンロード。

2. font0.dat〜font5.datと同じ場所に置く(例えば/usr/local/grads/data/)
  • 使用方法

 set font 6
 draw string x y hogehoge

例:El Nino

 set font 0
 draw string 1 1 El Ni`7n`0o



その他


複数のファイルからの作図


複数のctlやncを開き、ベクトルとコンターなどを一緒に書こうとすると、エラーとなる場合がある。
これは開いている複数のctlファイルの格子情報(例えば気圧面など)が異なるときに起こる。
二つ目以降のファイルから図を描く前に先に開いていたctlファイルを閉じ、x,y,z,t,gxoutなどをsetし直せば書ける。



open hoge1.ctl
set x,y,z,gxoutなど
d var
close 1
open hoge2.ctl
set x,y,z,gxoutなど
d var

×

open hoge1.ctl
open hoge2.ctl
set x,y,z,gxoutなど
d var
d var.2


コントロールファイルの表示 q ctlinfo

q ctlinfo

コントロールファイルの書式が表示される。NetCDFファイルからbinにデータを書き出し、
対応するコントロールファイルを作成するときに、これを用いればよい。

q (query)

コマンド履歴の表示 his


GradsWiki

Unixコマンドとしてのhistoryと同様、GrADS上でも打ちこんだコマンドの履歴が表示される。
history

または
his

今まで打ったコマンドをスクリプトに書き留める。
his hoge.gs

テキストに書き留める。
his hoge.txt

数字から始まるファイル名は指定できない。



コマンドライン、シェルでGrADSスクリプトを実行


コマンドラインからgsを走らせる
grads -lc hoge.gs

シェルの中でgsを作成したい場合は、

#!/bin/bash

cat>hoge.gs<<EOF
'reinit'
・・・
'quit'
EOF

grads -lc hoge.gs

のように”cat>hoge.gs<<EOF”と”EOF”で囲えばよい。

quitを書いておけば、gradsの終了も行える。
図を描き、表示結果を確認してから終了したいときは、quitの前の行にpull hogeを書けば、
表示されたのち、コマンドラインでエンターを打つとquitされる

fwriteなどで繰り返しの作業の場合は、gradsの起動オプションにbをつけると、黒画面の立ち上げを省略できて効率的。

#!/bin/bash
mon=1
while [ ${mon} -le 12 ] ;do

cat>hoge.gs<<EOF
'reinit'
'open ${mon}.ctl'
・・・
'quit'
EOF

grads -blc hoge.gs
rm -f hoge.gs
mon=$((mon+1))
done



図を描く位置・範囲 parea


東北大学大学院理学研究科 流体地球物理学講座 公開情報/GrADS

図の枠を描く位置を(x始点 x終点 y始点 y終点)で決める。最大の範囲はset vpageで変更できるが、デフォルトでは
0.0 11.0 0.0 8.5
となっている(最大)。
図の枠の位置だけを決めるため、中央に寄せればラベルやカラーバーを画面に収めることができる。

set xlopts 1 4 0.2
set ylopts 1 4 0.2  # ラベルを大きくすると 
d var  # 左や下にはみ出してしまう…
c
set parea 1.5 9.5 1.0 7.5
d var  # 図の左右に1.5ずつ、上下に1.0ずつスペースを入れたため、ラベルが画面に収まる

これを利用して、1つの画面に複数の図を描くことができる。

set parea 0.5 5.5 0 8.5
d var(t=1)
set parea 6.0 10.5 0 8.5
d var(t=2)

mul.gsを使うと楽。

弱い風ベクトルの非表示

d maskout(ua, mag(ua, va)-2);va

maskoutとmagを使用する。
maskoutはmaskout(変数,条件) で条件<0のときマスクアウトされる。
magはmag(変数1, 変数2) で変数1と変数2の√2乗和が計算される。
つまり上記で"2m/s以下の風をマスクアウトしたベクトル表示"ができる。

エンディアンの統一


GrADSの問題ではないが、fwriteでbinを書き出すときと、
ctlファイルで読み取るときのエンディアンを統一させないといけない。

エンディアンについてはここやここ

リトルエンディアン
bin書き出しにて

set gxout fwrite
set fwrite -le hoge.bin

ctlファイルにて
options little_endian

ビッグエンディアン
bin書き出しにて

set gxout fwrite
set fwrite -be hoge.bin

ctlファイルにて
options big_endian

気候値の時系列図


気候値の時系列(例えば、1月から12月までの気候値)を作図しようとすると、
GrADSでは図の左端にt=1, 右端にt=12が描かれるため、t=12から1にかけての変化が描けない。
そこで、modifyを使って周期的なデータにする。

climという気候値なら

set t 1 12
modify clim seasonal

これでclimのt=13〜にはt=1〜が周期的に入る。

set t 1 13  # 図の右端にt=12からt=13(つまりt=1)にかけての変化が描かれる
d clim

尚、時間軸ラベルを月だけにしたいなら

set tlsupp year
d clim

アニメーション軸の変更


デフォルトではt方向にのみアニメーションされるが、
x, y, z方向にアニメーションすることもできる。

デフォルト(t方向)

set x 1 144     #省略可
set y 1 73      #省略可
set t 1 12
set loopdim t   #省略可
d var           #t方向のアニメーション

xy断面のz方向アニメーション

set x 1 144     #省略可
set y 1 73      #省略可
set z 1 9
set loopdim z
d var           #z方向のアニメーション

xz断面のy方向アニメーション

set x 1 144     #省略可
set z 1 9
set y 1 73      #省略可
set loopdim y
d var           #y方向のアニメーション

何も描かない


set frame off  # 図の枠を描かない
set xlab off  # x軸のラベルを描かない
set ylab off  # y軸のラベルを描かない
set ccolor 0  # コンターの色を背景と同じにする
set clab off  # コンターのラベルを描かない
set grid off  # グリッドを描かない
set mpdraw off  # 背景の地図を描かない
set grads off  # GrADSのロゴとタイムスタンプを描かない
d var # すると、何も描かれない(笑)

何かだけ描きたいときに使えるかも?

特殊なコマンド


GrADSスクリプトライブラリ

gradsスクリプトが置いてあるディレクトリ/usr/local/grads-??/lib/に、

コマンド.gs

という名前で中にコマンドを書いておけば(+ .bashrcに
export GASCRP=/usr/local/grads-??/lib と書いておけば)、grads上で
コマンド

を打つだけで中に書いてあることが実行される




















コメントをかく


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

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

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