* GrADSメモ [#bf4d626d]

#contents

** できたらいいな。 [#z0e70bf5]
ご存知のことがあれば、コメントで教えてくださいm(_ _)m

- ある2つの面的な時系列データの、同じグリッド同士の時系列データの相関係数分布の表示。
- ある面的データの、球面調和関数スペクトルの算出
- テキスト出力のファイルへの保存。 -- [[kei]] &new{2007-08-01 (Wed) 13:33:06};
- きれいな画像の保存(printimの解像度はずいぶん低いので…) -- [[kei]] &new{2007-08-01 (Wed) 13:34:39};

#comment
** Script関係 [#ef94b690]

*** 引数を使う [#c8634a3b]
function main(args)をスクリプトの頭で定義する。
 (example)
 function main(args)
 var1=subwrd(args,1)
 var2=subwrd(args,2) 
 ...
** 注意点アレコレ [#ea5481a4]
- sumとaveの違い
 ave(var,t=1,t=10)*10 と
 sum(var,t=1,t=10) は同じではない!原因は欠損値の取り扱い(下記参照)。
- 時間平均→空間平均にすべし。
 sum(aave(var,lon=0,lon=360,lat=-90,lat=90),t=1,t=10) と
 aave(sum(var,t=1,t=10),lon=0,lon=360,lat=-90,lat=90) は異なる。後者が正解。
- 2変数以上の積の変数を平均(もしくは積分)する際、変数の積→平均(積分)すること。
 F=A*Bのとき 
 (Fの平均)=(Aの平均)*(Bの平均)+(A及びBの偏差の共分散)。
 相関が0のときは上の共分散もゼロ。共分散=(A,Bの相関)*(Aの標準偏差)*(Bの標準偏差)
 
- 一方、2変数以上の和なら順序が逆でもOK。(まあ、当たり前)→欠損値に注意(下記参照)
 F=A+Bのとき
 (Fの平均)=(Aの平均)+(Bの平均) or (A+Bの平均)

- 上記は、偏微分式でも適用可能。すなわち
 F=▽・(A,B)=∂A/∂x+∂B/∂yのとき (例えば収束)
 (Fの平均) = ∂(Aの平均)/∂x+∂(Bの平均)/∂y = (∂A/∂x+∂B/∂yの平均)

- 鉛直積分収束量計算:鉛直積分→水平微分の順。
 数学的には
 ∫(∂(q*u)/∂x)dp と ∂(∫(q*u)dp)/∂x (鉛直積分⇔水平微分)
 は同値だが、GrADSにおいては地表面付近の∂(q*u)/∂x がおかしくなるため
 (欠損があるため)、∂(∫(q*u)dp)/∂x (鉛直積分→水平微分)を使ったほうが良い。

- 欠損値の取り扱い
 変数同士の演算 → 全て欠損値を返す(要注意!)
 sum → 欠損値をゼロとして(欠損値を飛ばして)計算 
 ave → sumの結果を欠損値を除いた個数で割る。

- 欠損値の取り扱い2
 変数同士の演算ですべて欠損値を返すため、
 ave(var1,t=T1,t=T2)-ave(var2,t=T1,t=T2) と ave(var1-var2,t=T1,t=T2)
 は異なることがある。
 ave(var1-var2,t=T1,t=T2) と等しくなるのは
 ave(maskout(var1,var2+1e30,t=T1,t=T2)-ave(maskout(var2,var1+1e30),t=T1,t=T2)
 である。1e30というのは、var1のとりうる範囲が-1e30以上、ってことにしているから。
** 気圧重みづけ鉛直積分 [#bfb00b47]
 d vint(start pressure (hPa),var,end pressure (hPa))
**複数のctlファイルを開けて,それを選ぶ [#s8b27839]
 > set dfile 番号
**ある値以下{以上}をカットして表示. [#t0dc9ce8]
 > d maskout(var,{-}var {値})
**シェード{コンター}で表示 [#r8dfc3a3]
 > set gxout shaded {contour}
** 平均 [#i1efbb7b]
 > ave(var,t=1,t=24)など.ネスティングも可能
** 演算 [#k48da11c]
- 変数同士の四則計算,ave,sum,log,pow(x,y)なんかが使用可能.
** カラーバー表示 [#ifdadf20]
 > run cbar --> fg,bgなし
 > run cbarn --> fg,bgあり
** ヴェクター表示。 [#b6aef020]
 > d u;v
 間引きしたいときは
 > d skip(u,2,2);skip(v,2,2)
** プリントアウト [#h478b559]
 > enable print hoge.gx
 > print
 > disable print
 > ! gxps -c -i hoge.gx -o hoge.ps
 -c -rで黒背景。

''追記'':printimというコマンドを使えば一発。
 printim tmp.gif gif white

** テキスト出力 [#q8b92338]
 set gxout print
 set prnopts %10.3e 1 1 # 小数点以下3桁、一列に表示、ブランクは一文字分
画面に出たものをファイルに保存するやり方がわからん!コピペするしかないのか??スクリプトなら何とかなりそうだけど。。。
** タイトルつける [#l2b308c5]
 > draw title xxxxxxxxxxx
** netCDFを読む [#v6427ef0]
 % gradsnc
 > sdfopen xxx.nc
 > q file など
** バイナリで保存。(データの切り出しなど。) [#yfd495b1]
 > set fwrite xxx.dat
 > set gxout fwrite
 > set x (y,z,t, etc.) x1 x2
  ここで、xを設定しなおさないと、なぜかxmax+1まで保存されてしまう。
  時系列に保存したい場合は>set t t1 t2 など。
 > d var
 > disable fwrite
** 収束->ラクチン! [#a40ec129]
 > hdivg(x方向フラックス,y方向フラックス)
 鉛直積分値があればそのまま使えるし、下のvintを使って各高度の値から鉛直積分を求めても使える。
** 可降水量 [#nda794b2]
 > vint(地表気圧(mb単位), 変数(比湿等), 上がどこまで(定数))
 多分、ある層からある層までを区切る、なんてことはできない!?
 気圧を重力加速度で割る、という作業はすでにされているようだ。
** 背景を白に [#u0ae3517]
 > set display color white
** 中間の値に色を塗らない [#o3eb4be0]
 > set black -0.2 0.2 (-0.2から0.2は塗らない)

** (線の)凡例をつける。 [#wadc5e9d]
- cbar_l.gs/cbar_lineが便利。http://www.iges.org/grads/gadoc/library.html
- デフォルトの線種をそのまま使用したい場合は、
 > cbar_l -t "line1" "line2" "line3" -p
で、画面のクリックした位置にレジェンドを付加できる。線種をカスタマイズしたいときは
 > cbar_line -c 色 -m マーカ -l ライン -t "テキスト" -p
として、各項目に番号を振る。ちなみにGrADS標準は下記の通り。
|CENTER:|CENTER:|CENTER:|CENTER:|c
|Num|Color|Line|Mark|h
|1|1|1|2|
|2|3|1|3|
|3|7|1|4|
|4|2|1|5|
|5|6|1|1|
|6|9|1|2|
|7|10|1|3|
|8|11|1|4|
|9|12|1|5|
|10|15|1|1|

** アニメーション [#wa2e5324]
普通に
 set t 1 10
 d var
でできるが、コマ送りが早すぎなどという問題点がある。ftp://grads.iges.org/grads/scripts/xanim.gs を使うとさらに便利。
 set t 1 10                  # 同じ
 xanim -pause -grfill var    # マウスクリックでコマ送り(かつgrfillで表示)
 xanim -skip 4 -repeat 5 var # 4コマ目ごと、5回繰り返し
など。オプションは変数より前にないとダメらしい。
** 画面内の分割 [#r45672ca]
page.gs (ftp://ftp.cpc.ncep.noaa.gov/wd51we/grads/page.gs)が使える。
 page q1 # 右上
 d var1
 page q2 # 左上
 d var2
という具合。元に戻すには
 page reset
** 移動平均の図化 [#v1a0d530]
 set t 1 100               # 例えば。
 define a=ave(var,t-0,t+9) # 前方10個平均。define省略可。第二引数はtだけだとダメ。
 set missconn on           # 欠損は無視してつなぐ。
 d skip(a,10)

** 領域平均値の時系列グラフ [#o6ebda69]
 set t 1 100
 define a=aave(var,lon=-180,lon=180,lat=-90,lat=90) 
     # 面積重みつき全球平均
 set x 1 # Lineを書くため、Xを固定(値は何でも良い)
 d a
- aaveは、値を返してくれるが直接図化はされないため、一度defineする必要がある。

''追記'':tloopを使うとdefineする必要はない。
 set x 1
 set y 1
 d tloop(aave(var,lon=-180,lon=180,lat=-90,lat=90)
** 統計情報の表示 [#qd5ec259]
 set gxout stat
 d var
  欠損値の数や和、標準偏差、分散などが表示される。
  sublin, subwrdなどで取り出して使用するのが基本。

** 相関係数の算出(時間変動) [#p1c5e4ec]
ある時間軸一次元のデータの、領域全体の時系列変化に対する相関係数を計算する。
 set t 1 100
 define tmp=aave(var,lon=-180,lon=180,lat=-90,lat=90)
 set t 1
 d tcorr(tmp,var,t=1,t=100) 
  # 全球平均値tmpの変動に対しての、全グリッドの時間変動の相関を表示

** 相関係数の算出(空間変動) [#nc91edce]
同領域2変数の空間相関係数を計算する。
 d scorr(var1,var2,lon=-180,lon=180,lat=-90,lat=90)
  # var1とvar2の全球空間分布の一致度(相関)を表示
同様に、回帰直線の傾きは以下のようにして求められる。
 d sregr(var1,var2,lon=-180,lon=180,lat=-90,lat=90)

** 両対数グラフの作成 [#oc9fb30a]
GrADSには、Z軸に関してのみ対数軸が用意されている(set zlog on)が、X,Yに関しては存在しないため少々面倒。

対策は2通りある。ひとつは散布図(set gxout scatter)を使う方法、もうひとつは対数で描画したい横軸をZ軸としてデータ及びctlファイルを作成しset xyrev on/set zlog onを使う方法。前者は簡単だが点と点を結ぶ線を引けない。後者はデータを作り直さなければならず少し厄介だが、きれいな図が描ける。
*** 両対数グラフ/scatter使用時 [#w8f50a04]
散布図にしてXの値・Yの値にlogをかける。
 set gxout scatter
 set vrange X1 X2   #通常の線グラフならY軸に対応するが、散布図ではX軸
 set vrange2 Y1 Y2
 set digsiz 0.1   #デフォルトの点は若干小さいので大きくする。
 d log10(lon);log10(var)
*** 両対数グラフ/zlog使用時 [#q8d563b4]
描きたい図の横軸をZ軸にしたデータとコントロールファイルを作成。(元々のZ軸はX軸にしたりする。)
スペクトルの表示に便利。
#ref(KE3_18_w.gif,right,around,zoom,160x120)
 set x 1       #(実際は)高度を指定
 set z 1 300     #両対数グラフのX軸とするべき軸の範囲を指定
 set zlog on     #Z軸を対数表示
 set xyrev on    #描画領域のXYを入れ替え
 set xlevs 1 2 5 10 20 50 100 200 300 #X軸のアノテーションを指定(指定ナシも可)
 set cmark 0     #マークを描かない
 set cthick 6    #線を太く。
 d log10(var)    #変数varを対数表示。

** 平均日変化グラフの描画(平均季節変化グラフも同様) [#v6f47a7c]
Diurnal PatternやSeasonal Climatologyを算出するにはaveに時間ステップを指定してそれぞれの平均を取ると良い。
 set t 1 4  # tを日変化なら4ステップというように適当に設定。
 set x 1
 set y 1     # tloop用。
 d tloop(ave(aave(var,lon=X1,lon=X2,lat=Y1,lat=Y2),t-0,t+120,4))
    # ミソはaveの設定。はじめに設定したt=1,4のそれぞれから4つ飛び
    # で平均を取る。
あーラクチンだねー。

** グラフの中にグラフを書く [#ec973226]
vpageを使う。ただ、位置などの微調整が難しい。
#ref(KE_9_2.gif,right,around,zoom,160x120)
 # 現在の座標情報を得る(あくまでも目安)
 q gxinfo
 # gxinfoの情報を元に、背景色の四角を書く(透けるのを避けるため)
 # 次のvpageの座標と同じにしても良いが、余白を少なくするには微調整が必要。 
 set line 0
 draw recf xlo ylo xhi yhi
 # 続いてvpageを設定。良い位置にするにはかなり細かい調整が必要。
 set vpage xlo xhi ylo yhi
 # X-Yラベルのフォントサイズを変更(そのままでは縮小される)
 # デフォルトが 色番号=1、太さ=4、大きさ=0.12 らしいので、参考にする。
 set xlopts 1 2 0.3
 set ylopts 1 2 0.3 
 # gradsのロゴなどをはずす
 set grads off
 
 そして普段どおり描く。塗りつぶしの図なんかの場合は背景の四角はいらないね。
 vpage以前の設定は持ち越されないので注意。