Backup diff of IT memo/GrADS memo vs current(No. 51) - Kei YOSHIMURA in AORI/IIS U-Tokyo

The added line is THIS COLOR.
The deleted line is THIS COLOR.
#author("2020-03-08T09:58:54+09:00","kei","kei")
* GrADSメモ [#bf4d626d]
個人的なメモ書きです。別の（楽な）方法があればぜひ教えてください。誤りなどの指摘もよろしくお願いします。

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

- [Done]ある2つの面的な時系列データの、同じグリッド同士の時系列データの相関係数分布の表示。
-- 始めまして。筑波大のKamaeと言います。このページの情報はいつも参考にさせて頂いています。
&br;>ある2つの面的な時系列データの、同じグリッド同士の時系列データの相関係数分布の表示。
&br;というのは、Bin Guanさんのfcorr.gsでできると思います。
参考になれば幸いです。
http://mywiki.jp/s0721114/GrADS-Note/GrADS%83X%83N%83%8A%83v%83g%83%89%83C%83u%83%89%83%8A/#12 -- [[Kamae]] &new{2008-11-10 (Mon) 04:41:46};
-- おお、これは知りませんでした。情報提供ありがとうございます。このライブラリは強力そうです。Kamaeさんのページにも今後頻繁にご厄介になっていくと思います。よろしくお願いいたします。 -- [[kei]] &new{2008-11-10 (Mon) 11:52:14};
-- [[Bin Guan's GrADS Script Library:http://www.atmos.umd.edu/~bguan/grads/GrADS_Scripts.html]]　[Kamaeさんから情報提供]
- ある面的データの、球面調和関数スペクトルの算出
- [Done]テキスト出力のファイルへの保存。 -- [[kei]] &new{2007-08-01 (Wed) 13:33:06};
-- 自作した。[[./textout]]
-- UDFを使うともっと気の利いたものが作れそう。たとえば[[ここ:http://web.sfc.keio.ac.jp/~masudako/class/computer/grads/write.html]]のような。
-- Fさんから[[情報提供>IT memo/GrADS memo#q8b92338]]。問題解決。
- きれいな画像の保存（printimの解像度はずいぶん低いので…） -- [[kei]] &new{2007-08-01 (Wed) 13:34:39};
-- [[ここ:http://wind.geophys.tohoku.ac.jp/index.php?%B8%F8%B3%AB%BE%F0%CA%F3%2FGrADS%2FGrADS%A4%CETips#content_1_11]]で書かれているように、epsとして保存するのがベター。 --  &new{2007-11-07 (Wed) 10:50:21};
- [Done]鉛直断面図に、地形マスクをかける。 --  &new{2007-08-15 (Wed) 14:55:58};
-- できたっぽい。[[ここ>IT memo/GrADS memo#zca6364e]]参照 --  &new{2007-11-07 (Wed) 02:17:21};
- 気圧面データから高度毎の分布を描く。 -- [[kei]] &new{2007-11-07 (Wed) 02:37:56};
- bufrファイルの図化。bufrファイルを本当にGrADSで読めている、という決定的現場を見てみたい。困っている人の[[様子>IT memo/GrADS memo/Station data#ib90ae4a]] -- [[kei]] &new{2008-07-16 (Wed) 17:57:32};
-- 格子情報は読める場合があることがわかったが、観測値は扱えない。よって、何らかの形で観測値を吸い出して、再度station data用のフォーマットに書き換えて、という処理によって、可視化に成功した。要するにprepbufr（NCEPのBufr）はGrADSでは読めません。僕には。 -- [[kei]] &new{2009-03-13 (Fri) 16:29:15};
- station dataの散布図。及び、station地点の空間演算。イメージは、こんな感じ。
 stnaave(stnvar1-stnvar2,lon=x1,lon=x2,lat=y1,lat=y2)
できたら便利ですよねえ。（注意：上のコマンドは空想のものです。） -- [[kei]] &new{2008-07-16 (Wed) 17:59:23};
- 既に解決済みかもしれませんが、参考になれば幸いです。
den@academic-express.com
&br;テキスト出力の項で
set gxout print
set prnopts
d var
ok=write('test.txt',result,<append>)
で画面に表示されるものがtest.txtに保存されます。
（画面に表示されるものがresultという特別な変数に保存されています） -- [[田　寛之]] &new{2009-04-13 (Mon) 21:29:39};
--- ありがとうございます。助かります～。 -- [[kei]] &new{2009-04-14 (Tue) 06:36:50};

#comment2_kcaptcha
** 初級編 [#ne1495fb]
*** 入手とインストール [#q08e343d]
- http://www.iges.org/grads/downloads.html から自分にあったものをダウンロードする。バイナリでとってくればコンパイルいらずなので便利。Windowsでも一応動くが、重い。
- 好きな場所にバイナリを保存し、環境変数GADDIR, GASCRP（, GAUDFT）をセットする。GAUDFTは、ユーザー定義関数を作ったときに必要になる。csh系では、~/.cshrcにて以下の変数を定義する。
 setenv GADDIR /home/kei/GrADS/lib
 setenv GASCRP /home/kei/GrADS/lib
 setenv GAUDFT /home/kei/GrADS/udft
*** 基本コマンド [#ka13f62e]
- open '''controlfile'''
-- netCDFの場合は sdfopen '''ncfile'''
- d '''expr'''
- set x (y, z, t) '''num1''' '''num2'''
- set lon (lat, lev) '''val1''' '''val2'''

*** 出力形式 [#r8dfc3a3]
 set gxout outputtype
- よく使うのは、contour, grfill, shaded, vector, scatter　等
- ベクトル表示
 d u;v
 d skip(u,2,2);skip(v,2,2)　# 間引き

*** 演算 [#k48da11c]
- 変数同士の四則計算 
 d var1*var2
- 組み込み関数（ave，sum，log，pow(x,y), sqrtなんかが使用可能）
 d ave(var,t=1,t=24) # ネスティングも可能
-気圧重みづけ鉛直積分（可降水量など）
 d vint(sfcprs(hPa),expr,end pressure (hPa))
 d vint(pressfc/100,spfhprs, 10) # 例
-収束
 d hdivg(xexpr,yexpr)
*** 複数のctlファイルからデフォルトを選ぶ [#s8b27839]
 set dfile num

空間や時間解像度の異なる異なるctlファイルをデフォルトに選んだ後lonlatを指定し直せば、図化できるようになる。
*** ある値以下（以上）をマスクして表示． [#t0dc9ce8]
 d maskout(var,var-10) # 10以上の値を表示
*** 背景を白に [#u0ae3517]
 set display color white
*** 中間の値に色を塗らない [#o3eb4be0]
 set black -0.2 0.2　# -0.2から0.2は塗らない
- 背景色（色番号１）を塗っているので透過はしないことに注意。
*** カラーバー表示 [#ifdadf20]
 run cbar # fg,bgなし
 run cbarn # fg,bgあり
- オプション
 run cbarn 大きさ　水平0／鉛直1　x位置　y位置
*** タイトル、軸ラベルをつける [#l2b308c5]
 draw title zs
 draw xlab (ylab) xxxxx
*** プリントアウト [#h478b559]
 enable print hoge.gx
 print
 disable print
 ! gxps -c -i hoge.gx -o hoge.ps
  # -c -rで黒背景。

''追記''：printimというコマンドを使えば一発（スクリーンショットの取り込み）。
 printim tmp.gif gif white

''さらに追記''：print XXX.epsとすると、自動的に上記の作業をやってくれる。-Rオプションをつけると正常な向きになる。
 print -R foo.eps
*** バイナリで保存。（データの切り出しなど。） [#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

*** 鉛直p速度から鉛直風速への換算 [#b7df4171]
-w=-omega*Rd*(T+0.61*q)/p/g~
where omega(Pa/s), w(m/s), T(K), q(kg/kg), p(Pa), Rd=287(J/K/kg), g=9.8(kg/m/s/s)

In GrADS,
 > d -vvelprs*287*(tmpprs+0.61*spfhprs)/lev/100/9.8

*** 空気密度（ρ）の算出 [#jd435a7b]
- rho=p/Rd/(T+0.61*q)~
where rho(kg/m3), p(Pa), Rd=287(J/K/kg), T(K), q(kg/kg)

In GrADS,
 > d lev*100/287/(tmpprs+0.61*spfhprs)
ちなみに気温0℃、気圧1000hPaでの空気密度は1.275(kg/m3)。

*** 地表面気圧の算出 [#kb493925]
- Psfc=P0*((T0-0.0065*Zsfc)/T0)^(g/Rd/0.0065)~
where Psfc(Pa), P0(Pa), T0(K), Zsfc(m), g=9.8(m/s2), Rd=287(J/K/kg)

In GrADS,
 > d prmslmsl*pow((tmpprs-0.0065*hgtsfc)/tmpprs,9.8/287/0.0065)
誤差があることに注意。
** 中級編（scriptを使ったり、デフォルト以外の関数を使ったり） [#m99c2ebc]
*** scriptで引数を使う [#c8634a3b]
function main(args)をスクリプトの頭で定義する。
 (example)
 function main(args)
 var1=subwrd(args,1)
 var2=subwrd(args,2) 
 ...

*** テキスト出力 [#q8b92338]
 set gxout print
 set prnopts %10.3e 1 1 # 小数点以下3桁、一列に表示、ブランクは一文字分
画面に出たものをファイルに保存するやり方がわからん！コピペするしかないのか？？スクリプトなら何とかなりそうだけど。。。

結局、fwriteでバイナリ出力したものをテキストでファイルに書き出すコードを作った。[[./textout]]

''<追加情報>''

 *NCEP/NCARの再解析からあるグリッドのデータをテキスト形式
 *で縦一列に書き出す方法
 'sdfopen hgt.1993.nc'
 'set lat 30'
 'set lon 125'
 'set t 1 365'
 'set lev 200'
 'set gxout print'
 'set prnopts %7.2f 1 1'
 'd hgt'
 rc=write(hgt.1993.30N125E.asci,result)
 rc=close(hgt.1993.30N125E.ascii)
 これで、200hPa,30N,125Eの高度場の値がhgt.1993.30N125E.asciiと
 いうファイルに書き込まれます。ヘッダーが1行入ってしまいますが.....
 (by 藤波さん＠HyARC)
*** （線の）凡例をつける。 [#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 ＃ ４コマ目ごと、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]
同領域２変数の空間相関係数を計算する。
 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)

これらを利用して、散布図に回帰直線を引きたい場合は、
 slp=sregr(var1,var2,globe)
 av1=aave(var1,globe)
 av2=aave(var2,globe)
 int=slp*-av1+av2
 d var1;slp*var1+int
という感じでできる。
*** 両対数グラフの作成 [#oc9fb30a]
GrADSには、Z軸に関してのみ対数軸が用意されている（set zlog on）が、X,Yに関しては存在しないため少々面倒。

対策は2通りある。ひとつは散布図(set gxout scatter)を使う方法、もうひとつは対数で描画したい横軸をZ軸としてデータ及びctlファイルを作成しset xyrev on/set zlog onを使う方法。前者は簡単だが点と点を結ぶ線を引けない。後者はデータを作り直さなければならず少し厄介だが、きれいな図が描ける。
- 両対数グラフ／scatter使用時
散布図にして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使用時
描きたい図の横軸を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つ飛び
 　　　＃ で平均を取る。
あーラクチンだねー。

''追記''

上記は要はClimatologyを計算していることに他ならない。それを用いてAnomalyを求めたいときは、
 define xxx=ave(var,t-0,t+120,12)
 modify xxx seasonal （もしくはdiurnal）
とする。
*** グラフの中にグラフを書く [#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以前の設定は持ち越されないので注意。

*** テキストファイルの読み込み: read(filename) [#kaee15c9]
 rec=read(filename)
 var1=sublin(rec,2)
 rec=close(filename)
などと記述すると、recという変数の一行目にfilenameというファイルのステータスが入り、2行目以降にファイルの内容が取り込まれる。その中から変数を引き込みたいときは、sublin, subwrdを用いる。（上記の例だとfilenameというファイルの一行目がvar1という変数として引かれている）

80行までの小さなASCIIファイルを読み込める。この機能を使うと、GrADS外のコマンドとの引数受け渡しがある程度可能になる。読み終わったファイルはcloseが必要。

関連：[[テキストファイルの出力:http://meteora.ucsd.edu/~kyoshimura/?IT%20memo%2FGrADS%20memo#q8b92338]]

*** 外部コマンドの呼び出し [#h0424e4e]
スクリプトやコマンドラインにおいて、びっくりマークから始めると外部コマンドを呼び出せる。リダイレクトも使える。スクリプトで使用する際に、＄や”をそのまま用いたい場合は、＼（バックスラッシュor￥マーク）でエスケープしてあげる必要がある。以下、awkの例。
 tt=1
 '! awk "NR=='tt+1'{printf(\"%4.1f\n\",\$3}" ./slpanm_eval.txt > tmp.txt'
** 上級編（UDF使ったり、station data使ったり） [#o4686eb8]
*** UDF設定方法 [#n8d866dd]
- UDFTファイルの場所を環境変数GAUDFTに記述
- 関数が増えるたびにUDFTファイルに追記していく。

*** 空間内挿（lterp） [#yed267f9]
- 便利な空間内挿関数。時間を動かせないのでそうしたい場合は少し面倒。もともとetaについてくる関数っぽい？ 
- Googleで探してもソースが出てこないので、ここに貼り付けてしまう。オープンソースでしょう。たぶん。&attachref(lterp.tar);
- Googleで探してもソースが出てこないので、ここに貼り付けてしまう。オープンソースでしょう。たぶん。&ref(lterp.tar);
- 使い方
 d lterp(var1.1,var2.2) ＃ var1.1: 内挿前　var2.2: 内挿先
 d lterp(var1.1,var2.2)-var2.2　＃　などという演算が可能となる。楽。
- 注意点
-- 鉛直内挿ではない。
-- 鉛直レベルが同じだと、鉛直方向にも同時作業してくれる。
-- 解像度が異なるデータで、ある点の時系列を同時に描画したいときには、前もってdefineする必要がある。
*** 任意の断面図 [#ea4a6ac9]
- gr2stn, collect, coll2grを使用。が、鉛直分布「しか」描けないなど、かなり限定的にしか使えない。
-スクリプト例
 'reinit'
 'open /home/kyoshimura/IsoRSM/cal10km/200702/r_pgb_prs.ctl'
 'set grads off'
 'set zlog on'
 'set x 1'
 'set y 1'
 'set lev 1000 200'
 lon1 = -123.22
 lon2 = -120.82
 lat1 = 38.61
 lat2 = 39.20
 lon = lon1
 'collect 1 free'
 while (lon <= lon2)
   lat = lat1 + (lat2-lat1)*(lon-lon1) / (lon2-lon1)
   'collect 1 gr2stn(pwat,'lon','lat')'
   lon = lon + 0.1
 endwhile
 
 'set x 14 16'
 'set xaxis 'lon1' 'lon2
 'set clab on'
 'set gxout shaded'
 'd coll2gr(1,-u)'
 'run cbarn'

*** 鉛直分布図などに地形マスクをかける [#zca6364e]
#ref(x-sec_06Z13FEB2007_38.571.gif,right,around,zoom,160x120);
- 参考図（→右の図）
- やり方
 set z 1 10
 define tmp=pressfc ＃　全層に表面気圧を入れる
などとやってから、 数多くの正値の太いコンターでd lev-tmp/100を表示させる。
 set clevs 0 5 10 15... ＃　コンターの間に空白がなくなるように、必要なだけ延々と書く。　
 set ccolor 15　＃　灰色
 set cthick 6　＃　太いほうが良い
 set clab off　＃　コンターラベルはなし
 d -(tmp/100-lev)　＃　描画
-失敗例
-- 鉛直図を描いた後、set z 1　などとして線などで標高を描こうとしても、軸がずれてしまう。
-- d maskout(lev-tmp/100,lev-tmp/100)などとやると、角ばったマスクになる。
-- gxout shaded でマスクしようとしても、「透過色」を指定できないのでうまくいかない。
*** 鉛直分布図と非鉛直分布図を同時に描く。 [#u109dfd9]
- 参考図（→上の図）
- やり方：vpageをうまく使って調整する。
- 例：&attachref(CZC-ALT.gs);参照
- 例：&ref(CZC-ALT.gs);参照

*** Gradient/Laplacianの計算 [#cdaf01c9]
hdivgとhcurlをつかって、スカラー量qのdq/dx、dq/dy及びラプラシアン(d2q/dx2+d2q/dy2)を計算できる。簡単にできるスクリプトはこちら。
- &attachref(lap.gs);
- &attachref(gradx.gs);
- &attachref(grady.gs);
- &ref(lap.gs);
- &ref(gradx.gs);
- &ref(grady.gs);
*** EOFの計算 [#r7208b9c]
#ref(gsm_AO_pat1.gif,right,around,zoom,160x120);
大変親切なことに、UCLAのMunnichさんの手によってeofgradsというとっても便利なツールがある。これが非常に便利。多謝。
http://www.atmos.ucla.edu/~munnich/Grads/EOF/

+ マシンに合ったバイナリを入手後、UTFT等を適宜設定し（上記のサイトに丁寧に説明されているため省略）、使用可能にする。
+ 基本的に、
 > set x 1 144
 > eof expression
で、expres_eof.ctl/datというファイルとexpres_pc.ctl/datというファイルが生成され、それぞれ固有ベクトル（EOF空間分布、Z方向に成分）と標準化時係数（T方向に時間）をあらわす。グローバルデータの場合、上記のようにxをグリッド数で再定義する必要がある。
+あとは、適宜それらのコントロールファイルを開き、図化すればよいだけ。

*** Station Dataの取り扱い [#k2a62912]
別ページへ[[IT memo/GrADS memo/Station data]]
*** アンサンブルデータの取り扱い NEW!! [#s790fcfe]
別ページへ[[IT memo/GrADS memo/GrADS2.0]]

*** Using OpenGrADS NEW!! [#f6633f39]
See separate page [[IT memo/GrADS memo/OpenGrADS]]

*** ハドレー循環（子午面質量輸送図）の描き方 [#wf1371bd]
#ref(hadley_IsoGSM_s1.gif,right,around,zoom,160x120)
- GrADS script for drawing a Hadley Circulation (meridional mass stream function).
- vintの3つ目の引数に変数を入れられないため、少々面倒。
- そのため、一層ずつ子午面での質量流線関数を計算している。
- 以下、hadley.gsの中身。
 '! rm -f hadley.dat'
 'reinit'
 
 'open /home/kyoshimura/IsoGSM1/pgb_mon-s1.ctl'
 
 tt=1
 'set t 'tt
 'set lon 0'
 'zps=ave(pressfc/100,lon=-180,lon=180)'
 'set z 1 17'
 'zv=ave(vgrdprs,lon=-180,lon=180)'
 
 zz=2
 while zz <= 17
 'set z 'zz
 'q dim'
 res=sublin(result,4)
 zlv=subwrd(res,6)
 'set fwrite tmp.dat'
 'set gxout fwrite'
 'd (vint(zps,zv,10)-vint(zps,zv,'zlv'))*2*3.14*6.37e6*cos(lat*3.1415/180)'
 'disable fwrite'
 '! cat tmp.dat >> hadley.dat'
 zz=zz+1
 endwhile
  
 'reinit'
 
 'open hadley.ctl'
 'set zlog on'
 'set xlopts 1 3 0.16'
 'set ylopts 1 3 0.16'
 'set xflip on'
 'set ylevs 900 700 500 300 150 100 70 40 20 10'
 'set gxout shaded'
 'set grads off'
 'set clevs -300 -250 -200 -150 -120 -90 -60 -30 -10 10 30 60 90 120 150 200 250 300'
 'rgbset2'
 'd msf/1e9'
 'set gxout contour'
 'set clevs -300 -250 -200 -150 -120 -90 -60 -30 -10 10 30 60 90 120 150 200 250 300'
 'set ccolor 1'
 'd msf/1e9'
 'draw ylab pressure (hPa)'
 'draw title Meridional m.s.f. (10e9kg/s), IsoGSM'
 'run cbarn'
- hadley.ctl 
 dset ^hadley.dat
 options little_endian
 undef -999.
 xdef  1 linear    0.000  2.500
 ydef    73 linear  -90.000  2.500
 zdef    16 levels 
      925     850     700     600     500     400     300     250     200
      150     100      70      50      30      20      10
 tdef    1 linear              00Z1JAN1979     1mo
 vars    1
 msf 16 16 mass stream function
 endvars
*** 相対湿度の計算（自作UDF） [#w82f96e7]
http://www.cactus2000.de/js/calchum.pdf　に基づいて相対湿度を計算するUDFを作った。元の論文は、Lowe, P.R. and J.M. Ficke (1974)らしい。

''インストール''
+&ref(sp2rh.tar); をudfディレクトリ内で解凍し、
+適当なコンパイラでsp2rh.fをコンパイル（エンディアンには注意）し、
+udftにsp2rh.udftをコピペする（パス等は要修正）。

''使い方''
 d sp2rh(temp,pressure,spec_hum)
 # where temp[K], pressure[mb], spec_hum[kg/kg]. 
 # if you want to use constant value, 
 # use dummy variable as "pressfc-pressfc+100". (similar to vint)
** NetCDFファイルの読み込み [#k92a27aa]
*** COARDS規約のファイル [#if3f099e]
 sdfopen XXX.nc
するだけ。
*** 読めそうだけどsdfopenが効かない場合 [#d30c4d73]
簡素版コントロールファイルをつくってxdfopenする。
 xdfopen test.ctl
[[SWINGデータ:http://atoc.colorado.edu/~dcn/SWING/database/clim_s1a/monthly_data.php]]用のサンプル→&attachref(s1a.echam4.ctl);
[[SWINGデータ:http://atoc.colorado.edu/~dcn/SWING/database/clim_s1a/monthly_data.php]]用のサンプル→&ref(s1a.echam4.ctl);
- templateオプションが効くが、たまにSegmentation Faultしてしまう。（特に、演算による時系列アニメで、ファイルをまたぐ場合）

-xdfopen用ctlファイルの作り方
 DSET ^foo.nc  # netcdfファイルを指定
 UNDEF 1e20    # undef値を指定
 XDEF xc       # netcdfファイルのX軸変数を指定（lon, longitude, xcなど）
 YDEF yc       # 同様にY軸（lat, latitude, ycなど）
 ZDEF plev     # 同様にZ軸（plev, lev, heightなど）。無い場合は行自体削除。
 TDEF time     # 時間軸変数を指定。
 VARS 1        # 変数の数
 ta=>ta 25 25 ta # NCファイル内の変数の名前=>GrADSでの名前、レベル数
 ENDVARS
XDEF、TDEFなどで変数を指定した後、
 XDEF xc 3 226.0 1.25   # 226度から1.25間隔で3点
 YDEF yc 3 levels 10 11.25 12.5 #10, 11.25, 12.5度の3点
などと指定することができる。
** Office使いのための論文用画質維持について [#t14eb7c4]
- 複数の図を一枚のパネルにする際、パワポで切り貼りしてワードに貼り付けるという手法は、論文だけでなくそのまま発表にも使えて便利。
- が、論文最終段階で、そのパネルを画像（Tiff or eps）として提出する場合、画質がかなり落ちる。（画面上の画質と同じになる（？））
- 対処法としては、
-- パワポで、複数の図からなるパネル（背景含む）を全部選択し、右クリックで「画像として保存」を選択、「拡張メタファイル（.emf）」で保存。
-- 次に、emfファイルを開く。（デフォルト設定でWindowsの画像Viewerが起動する）
-- 最後に、tiffファイルとして保存する。
- ただ、かなりむだに容量が大きくなるのが困り者。

[2009/5/16 更新]
-Powerpointで作った図のEPS化には、WMF2EPSというソフトが便利。デフォルトではVistaに対応していないが、下のリンクの手続きで対応可能。
--http://www.wmf2eps.de.vu/
--http://www.iplab.is.tsukuba.ac.jp/~zono/tips/wmf2eps_on_vista.html
--http://www.bunmeisha.co.jp/LaTeX2e/latex2e_eps.html
-難点
--拡張メタファイルではトリミングが無視される。そのため、トリミングを使った図はPNGなどに「形式を選択して貼り付け」し、EPS化するしかないようだ。
** EPS/PSで保存した場合の横線の謎[#ad1fa1be]
#ref(tmp.gif,right,around,zoom,160x120)
-http://www.sci.hokudai.ac.jp/~sasa/tips.html
-http://mausam.hyarc.nagoya-u.ac.jp/~hatsuki/grads/node16.html
-http://www.atmos.ucla.edu/~munnich/grads/gxgif.html

で触れられているように、gxepsやgxpsで作った図には変な横線が入ることがある。ベクトルデータの調査の結果、shadedの図では、上下に一グリッド分の細長いポリゴンをぺたぺた重ね塗りしているようで、上下のポリゴンになぜか隙間ができているのが原因の模様。はっきり言ってどうしようもない（多分）。

対策としては、
- shadedを使わない。contourかgrfillを使う。（grfillだと横線だけでなく格子線になるが、比較的目立たない。）
- EPS/PSを使わない。gifなどで論文にしてしまう。GrADSのEPSファイルって、結局ほぼラスターデータだから（文字やshadeが相当小さな単位でのベクトルデータ）、容量や画質はそんなに変わらん。
-- この記述はあいまいで、printimを使って出力されたスクリーンショットのgif画像を使う、ことを意味していると思われる。

くらいかなあ。

''<追加情報>''

gimpでセーブしなおすと横線は消えるらしい。（by T.Yamada@GSFC)

''<追加情報2>''

-http://www.kamiguchi.org/kenji/2006/03/grads-eps.html

の情報によると、アンチエイリアシングなるものの効果によるらしい。で、gsのアンチエイリアス用のオプションである
 -dTextAlphaBits=n
 -dGraphicsAlphaBits=n
のnを両方1に（TextのほうはGrADSに関しては実質効かない）して画像ファイルに変換すれば、確かに不快な横線は消える。が、今度は文字が粗くなったり線の太さがバラバラになったりする。どっちをとるか、だね。

&ref(with_antialias.gif,center,zoom,160x120);  
&ref(no_antialias.gif,center,zoom,160x120);
** フォントの作り方 [#ofde738b]
デフォルトでlib内にfont0.dat～font5.datが用意されている。font3.datが特殊文字。しかし、同位体屋に必須である「パーミル（‰）」が存在しない。よって自作を試みる。
われながら暇人だ。

- 現在使えるフォントの確認
 font number
で一覧が表示される。デフォルトでは０－５。チルダ（～）がどのフォントでもまったく同じ文字に当てられている。

- フォントファイルの解読

http://iges.org/grads/gadoc/font.htmlを参照。

fontX.datには、ASCIIコードの32（空白）～126番目（チルダ）の95のキーを入力したときに呼ばれる記号が、それぞれの行に記述されてある。つまり、あるキー（例えばチルダ）に対応する記号を変更したければ、その行（チルダの場合は95行目）を変更すればよい。

続いてフォントの中身を解読する。例えばスラッシュ（／：16行目）は
   3G][BIb
とある。まずはじめの3は、3つの座標ペアが続く、という意味。座標というのは、ASCIIコードの順番でちょうど真ん中の「R」をゼロとして、若いキーがマイナス、後のキーがプラス、としている。よってGは-11、]は+11が当てられる。

はじめの数字の直後は、記号の幅を示す。つまりG]は-11～+11を占めるということ。

これ以降は、一筆書きなら座標を並べていき、間をおくなら「空白＋R」でつなげる。よって例の場合は、[B（+9，-16）からIb（-9，+16）まで直線を引く、ということになる。
*** GrADSでパーミルを使う [#qe8ff029]
font3.datの最後の行（チルダ用）を以下のように変更する。
  14F^IUISJPLONOOPPQTTVUXUZT[Q[O
   ↓
  44D`WFE[ RJFLHLJKLIMGMEKEIFGHFJFLGOHRHUGWF RSTQUPWPYR[T[VZWXWVUTST R[TYUXWXYZ[\[^Z_X_V]T[T

GrADS内では、「｀３～」で呼び出す。例えば、
 ga> draw ylab `3d`0`a18`nO [`3~`0]
で、『δ18O [‰]』（18は上付き）が表示される。

*** Using Permil (part per thousand) in GrADS [#wc5250f9]
Modify the last line of '''gradsdir'''/lib/font3.dat as follows: 
  14F^IUISJPLONOOPPQTTVUXUZT[Q[O
   ↓
  44D`WFE[ RJFLHLJKLIMGMEKEIFGHFJFLGOHRHUGWF RSTQUPWPYR[T[VZWXWVUTST R[TYUXWXYZ[\[^Z_X_V]T[T

The permil symbol is written by "`3~", such as
 ga> draw ylab `3d`0`a18`nO [`3~`0]
This writes "d18O [permil]" in Y-axis.
** lats4dを使ったフォーマット変換 [#t434eeb9]
別項目へ[[IT memo/GrADS memo/lats4d]]。


** 注意点アレコレ [#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の結果を欠損値を除いた個数で割る。

- 欠損値の取り扱い２
 変数同士の演算ですべて欠損値を返すため、
 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以上、ってことにしているから。