IT memo/GrADS memo2
[
Front page
] [
New
|
List of pages
|
Search
|
Recent changes
|
Help
]
Start:
* GrADSメモ [#bf4d626d]
個人的なメモ書きです。別の（楽な）方法があればぜひ教えて...
Pukiwikiの調子が悪く、[[IT memo/GrADS memo]]が更新できな...
#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%...
-- おお、これは知りませんでした。情報提供ありがとうござい...
-- [[Bin Guan's GrADS Script Library:http://www.atmos.umd...
- [Done] ある面的データの、球面調和関数スペクトルの算出
-- OpenGrADSを使えば解決。run power.gs ugrdprs
- [Done]テキスト出力のファイルへの保存。 -- [[kei]] &new{...
-- 自作した。[[./textout]]
-- UDFを使うともっと気の利いたものが作れそう。たとえば[[...
-- Fさんから[[情報提供>IT memo/GrADS memo#q8b92338]]。問...
- きれいな画像の保存（printimの解像度はずいぶん低いので…...
-- [[ここ:http://wind.geophys.tohoku.ac.jp/index.php?%B8%...
- [Done]鉛直断面図に、地形マスクをかける。 -- &new{2007-...
-- できたっぽい。[[ここ>IT memo/GrADS memo#zca6364e]]参照...
- [Done]気圧面データから高度毎の分布を描く。 -- [[kei]] &...
-- zinterpを使えば解決。
- bufrファイルの図化。bufrファイルを本当にGrADSで読めてい...
-- 格子情報は読める場合があることがわかったが、観測値は扱...
- station dataの散布図。及び、station地点の空間演算。イメ...
stnaave(stnvar1-stnvar2,lon=x1,lon=x2,lat=y1,lat=y2)
できたら便利ですよねえ。（注意：上のコマンドは空想のもの...
- 既に解決済みかもしれませんが、参考になれば幸いです。
den@academic-express.com
&br;テキスト出力の項で
set gxout print
set prnopts
d var
ok=write('test.txt',result,<append>)
で画面に表示されるものがtest.txtに保存されます。
（画面に表示されるものがresultという特別な変数に保存され...
--- ありがとうございます。助かります～。 -- [[kei]] &new{...
#comment2_kcaptcha
** 初級編 [#ne1495fb]
*** 入手とインストール [#q08e343d]
- http://www.iges.org/grads/downloads.html から自分にあっ...
- 好きな場所にバイナリを保存し、環境変数GADDIR, GASCRP（,...
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ファイルをデフォルトに選...
*** ある値以下（以上）をマスクして表示． [#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とすると、自動的に上記の作業...
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(...
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=2...
In GrADS,
> d prmslmsl*pow((tmpprs-0.0065*hgtsfc)/tmpprs,9.8/287/0...
誤差があることに注意。
** 中級編（scriptを使ったり、デフォルト以外の関数を使った...
*** 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でバイナリ出力したものをテキストでファイルに...
''<追加情報>''
*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.as...
いうファイルに書き込まれます。ヘッダーが1行入ってしまい...
(by 藤波さん＠HyARC)
*** （線の）凡例をつける。 [#wadc5e9d]
- cbar_l.gs/cbar_lineが便利。http://www.iges.org/grads/ga...
- デフォルトの線種をそのまま使用したい場合は、
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:/...
set t 1 10 ＃ 同じ
xanim -pause -grfill var ＃ マウスクリックでコマ送り...
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省略可...
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は、値を返してくれるが直接図化はされないため、一度d...
''追記''：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 zlo...
対策は2通りある。ひとつは散布図(set gxout scatter)を使う...
- 両対数グラフ／scatter使用時
散布図にしてXの値・Yの値にlogをかける。
set gxout scatter
set vrange X1 X2 ＃通常の線グラフならY軸に対応するが...
set vrange2 Y1 Y2
set digsiz 0.1 ＃デフォルトの点は若干小さいので大き...
d log10(lon);log10(var)
- 両対数グラフ／zlog使用時
描きたい図の横軸をZ軸にしたデータとコントロールファイルを...
スペクトルの表示に便利。
#ref(KE3_18_w.gif,/IT memo/GrADS memo,right,around,zoom,1...
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を対数表示。
*** 平均日変化グラフの描画（平均季節変化グラフも同様） [#...
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+...
＃ ミソはaveの設定。はじめに設定したt=1,4のそれぞ...
＃ で平均を取る。
あーラクチンだねー。
''追記''
上記は要はClimatologyを計算していることに他ならない。それ...
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というフ...
80行までの小さなASCIIファイルを読み込める。この機能を使う...
関連：[[テキストファイルの出力:http://meteora.ucsd.edu/~k...
*** 外部コマンドの呼び出し [#h0424e4e]
スクリプトやコマンドラインにおいて、びっくりマークから始...
tt=1
'! awk "NR=='tt+1'{printf(\"%4.1f\n\",\$3}" ./slpanm_eva...
** 上級編（UDF使ったり、station data使ったり） [#o4686eb8]
*** UDF設定方法 [#n8d866dd]
- UDFTファイルの場所を環境変数GAUDFTに記述
- 関数が増えるたびにUDFTファイルに追記していく。
*** 空間内挿（lterp） [#yed267f9]
- 便利な空間内挿関数。時間を動かせないのでそうしたい場合...
- Googleで探してもソースが出てこないので、ここに貼り付け...
- 使い方
d lterp(var1.1,var2.2) ＃ var1.1: 内挿前 var2.2: 内挿先
d lterp(var1.1,var2.2)-var2.2 ＃ などという演算が可能...
- 注意点
-- 鉛直内挿ではない。
-- 鉛直レベルが同じだと、鉛直方向にも同時作業してくれる。
-- 解像度が異なるデータで、ある点の時系列を同時に描画した...
*** 鉛直内挿（zinterp） [#if5146f2]
- 便利な鉛直内挿関数。GASCRP内になさければ、ftp://grads.i...
- 気圧面データから高度面の値を出したり、高度面データから...
- 使い方（スクリプトそのものに丁寧に書かれています）
d zinterp(field, zgrid, zlev)
例えばNICAMのモデル面（地表からの高さ）での出力データから...
open ms_pres.ctl
open topog.ctl
d zinterp(lev,ms_pres,50000)+topog.2
となる。850hPaの東西風は
open ms_pres.ctl
open ms_u.ctl
d zingerp(ms_u.2,ms_pres,85000)
- 注意点
-- zgridとzlevは一緒の単位を持つ。
-- 線形内挿のため、高度面データから気圧面を作る場合、（気...
*** 任意の断面図 [#ea4a6ac9]
- gr2stn, collect, coll2grを使用。が、鉛直分布「しか」描...
-スクリプト例
'reinit'
'open /home/kyoshimura/IsoRSM/cal10km/200702/r_pgb_prs.c...
'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,160x...
- 参考図（→右の図）
- やり方
set z 1 10
define tmp=pressfc ＃ 全層に表面気圧を入れる
などとやってから、 数多くの正値の太いコンターでd lev-tmp/...
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をうまく使って調整する。
- 例：&ref(CZC-ALT.gs);参照
*** Gradient/Laplacianの計算 [#cdaf01c9]
hdivgとhcurlをつかって、スカラー量qのdq/dx、dq/dy及びラプ...
- &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と...
+あとは、適宜それらのコントロールファイルを開き、図化すれ...
*** 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 (meridion...
- 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*co...
'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 3...
'rgbset2'
'd msf/1e9'
'set gxout contour'
'set clevs -300 -250 -200 -150 -120 -90 -60 -30 -10 10 3...
'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...
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 に基づいて相対湿...
''インストール''
+&ref(sp2rh.tar); をudfディレクトリ内で解凍し、
+適当なコンパイラでsp2rh.fをコンパイル（エンディアンには...
+udftにsp2rh.udftをコピペする（パス等は要修正）。
''使い方''
d sp2rh(temp,pressure,spec_hum)
# where temp[K], pressure[Pa], spec_hum[kg/kg] and RH[%]
# if you want to use constant value,
# use dummy variable as "pressfc-pressfc+100000". (simil...
***比湿の計算（自作UDF） [#e9dff544]
上記のsp2rhの逆。
''インストール''
+&ref(rh2sp.tar); をudfディレクトリ内で解凍し、
+適当なコンパイラでrh2sp.fをコンパイル（エンディアンには...
+udftにrh2sp.udftをコピペする（パス等は要修正）。
''使い方''
d rh2sp(temp,pressure,relative_hum)
# where temp[K], pressure[Pa], relative_hum[%] and spec_...
# if you want to use constant value,
# use dummy variable as "pressfc-pressfc+100000". (simil...
** OpenGrADS [#v81436b6]
- 以下は古い情報（ver1.9）
経緯は忘れたがOpenGrADSのページからダウンロード・インスト...
> q udct
> q udft
でずらずらっと出てきていたら、インストール済み。
-新しい情報（ver2.0以降）
--http://sourceforge.net/projects/opengrads/files/
のgrads2から自分の環境のものをダウンロードする。（僕の場...
--中にあるContents/gradsを実行する。
$ ./Contents/grads
環境変数を変えていないため、基本的に、昔のLibraryやUDFTも...
*** 流線関数・速度ポテンシャルの描画 [#ab437288]
要OpenGrADS。
- 流線関数
> d fish(hcurl(ugrdprs,vgrdprs))
or
> d fish_psi(ugrdprs,vgrdprs)
- 速度ポテンシャル
> d fish(hdivg(ugrdprs,vgrdprs))
or
> d fish_chi(ugrdprs,vgrdprs)
*** 球面調和関数のパワースペクトル [#s504239c]
要OpenGrADS
> run power.gs ugrdprs
とやるだけで、波数に対するKinetic Energyが出る。また、
> d sh_filt(ugrdprs,10)
で、波数10でのフィルタリング。
長年、これがほしかった。
*** 持ち上げ凝結高度（LCL） [#bda257a0]
要OpenGrADS
> rh2m=100*pressfc*spfh2m/(0.622+0.378*spfh2m)/satvap(tm...
> plcl=plcl(tmp2m,rh2m,pressfc) ; calculate LCL pressure
> d zinterp(hgtprs,lev,plcl/100) ; calculate height for ...
OpenGrADSを立ち上げると、これまでの自作のUDFが使えなくな...
** NetCDFファイルの読み込み [#k92a27aa]
*** COARDS規約のファイル [#if3f099e]
sdfopen XXX.nc
するだけ。
*** 読めそうだけどsdfopenが効かない場合 [#d30c4d73]
簡素版コントロールファイルをつくってxdfopenする。
xdfopen test.ctl
[[SWINGデータ:http://atoc.colorado.edu/~dcn/SWING/databas...
- templateオプションが効くが、たまにSegmentation Faultし...
-xdfopen用ctlファイルの作り方
DSET ^foo.nc # netcdfファイルを指定
UNDEF 1e20 # undef値を指定
XDEF xc # netcdfファイルのX軸変数を指定（lon, long...
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ファイルを開く。（デフォルト設定でWindowsの画...
-- 最後に、tiffファイルとして保存する。
- ただ、かなりむだに容量が大きくなるのが困り者。
[2009/5/16 更新]
-Powerpointで作った図のEPS化には、WMF2EPSというソフトが便...
--http://www.wmf2eps.de.vu/
--http://www.iplab.is.tsukuba.ac.jp/~zono/tips/wmf2eps_on...
--http://www.bunmeisha.co.jp/LaTeX2e/latex2e_eps.html
-難点
--拡張メタファイルではトリミングが無視される。そのため、...
** 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...
-http://www.atmos.ucla.edu/~munnich/grads/gxgif.html
で触れられているように、gxepsやgxpsで作った図には変な横線...
対策としては、
- shadedを使わない。contourかgrfillを使う。（grfillだと横...
- EPS/PSを使わない。gifなどで論文にしてしまう。GrADSのEPS...
-- この記述はあいまいで、printimを使って出力されたスクリ...
くらいかなあ。
''<追加情報>''
gimpでセーブしなおすと横線は消えるらしい。（by T.Yamada@G...
''<追加情報2>''
-http://www.kamiguchi.org/kenji/2006/03/grads-eps.html
の情報によると、アンチエイリアシングなるものの効果による...
-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が用意されている。...
われながら暇人だ。
- 現在使えるフォントの確認
font number
で一覧が表示される。デフォルトでは０－５。チルダ（～）が...
- フォントファイルの解読
http://iges.org/grads/gadoc/font.htmlを参照。
fontX.datには、ASCIIコードの32（空白）～126番目（チルダ）...
続いてフォントの中身を解読する。例えばスラッシュ（／：16...
3G][BIb
とある。まずはじめの3は、3つの座標ペアが続く、という意味...
はじめの数字の直後は、記号の幅を示す。つまりG]は-11～+11...
これ以降は、一筆書きなら座標を並べていき、間をおくなら「...
*** GrADSでパーミルを使う [#qe8ff029]
font3.datの最後の行（チルダ用）を以下のように変更する。
14F^IUISJPLONOOPPQTTVUXUZT[Q[O
↓
44D`WFE[ RJFLHLJKLIMGMEKEIFGHFJFLGOHRHUGWF RSTQUPWPYR[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 f...
14F^IUISJPLONOOPPQTTVUXUZT[Q[O
↓
44D`WFE[ RJFLHLJKLIMGMEKEIFGHFJFLGOHRHUGWF RSTQUPWPYR[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の...
- 一方、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-var...
は異なることがある。
ave(var1-var2,t=T1,t=T2) と等しくなるのは
ave(maskout(var1,var2+1e30,t=T1,t=T2)-ave(maskout(var2,v...
である。1e30というのは、var1のとりうる範囲が-1e30以上、...
End:
* GrADSメモ [#bf4d626d]
個人的なメモ書きです。別の（楽な）方法があればぜひ教えて...
Pukiwikiの調子が悪く、[[IT memo/GrADS memo]]が更新できな...
#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%...
-- おお、これは知りませんでした。情報提供ありがとうござい...
-- [[Bin Guan's GrADS Script Library:http://www.atmos.umd...
- [Done] ある面的データの、球面調和関数スペクトルの算出
-- OpenGrADSを使えば解決。run power.gs ugrdprs
- [Done]テキスト出力のファイルへの保存。 -- [[kei]] &new{...
-- 自作した。[[./textout]]
-- UDFを使うともっと気の利いたものが作れそう。たとえば[[...
-- Fさんから[[情報提供>IT memo/GrADS memo#q8b92338]]。問...
- きれいな画像の保存（printimの解像度はずいぶん低いので…...
-- [[ここ:http://wind.geophys.tohoku.ac.jp/index.php?%B8%...
- [Done]鉛直断面図に、地形マスクをかける。 -- &new{2007-...
-- できたっぽい。[[ここ>IT memo/GrADS memo#zca6364e]]参照...
- [Done]気圧面データから高度毎の分布を描く。 -- [[kei]] &...
-- zinterpを使えば解決。
- bufrファイルの図化。bufrファイルを本当にGrADSで読めてい...
-- 格子情報は読める場合があることがわかったが、観測値は扱...
- station dataの散布図。及び、station地点の空間演算。イメ...
stnaave(stnvar1-stnvar2,lon=x1,lon=x2,lat=y1,lat=y2)
できたら便利ですよねえ。（注意：上のコマンドは空想のもの...
- 既に解決済みかもしれませんが、参考になれば幸いです。
den@academic-express.com
&br;テキスト出力の項で
set gxout print
set prnopts
d var
ok=write('test.txt',result,<append>)
で画面に表示されるものがtest.txtに保存されます。
（画面に表示されるものがresultという特別な変数に保存され...
--- ありがとうございます。助かります～。 -- [[kei]] &new{...
#comment2_kcaptcha
** 初級編 [#ne1495fb]
*** 入手とインストール [#q08e343d]
- http://www.iges.org/grads/downloads.html から自分にあっ...
- 好きな場所にバイナリを保存し、環境変数GADDIR, GASCRP（,...
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ファイルをデフォルトに選...
*** ある値以下（以上）をマスクして表示． [#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とすると、自動的に上記の作業...
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(...
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=2...
In GrADS,
> d prmslmsl*pow((tmpprs-0.0065*hgtsfc)/tmpprs,9.8/287/0...
誤差があることに注意。
** 中級編（scriptを使ったり、デフォルト以外の関数を使った...
*** 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でバイナリ出力したものをテキストでファイルに...
''<追加情報>''
*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.as...
いうファイルに書き込まれます。ヘッダーが1行入ってしまい...
(by 藤波さん＠HyARC)
*** （線の）凡例をつける。 [#wadc5e9d]
- cbar_l.gs/cbar_lineが便利。http://www.iges.org/grads/ga...
- デフォルトの線種をそのまま使用したい場合は、
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:/...
set t 1 10 ＃ 同じ
xanim -pause -grfill var ＃ マウスクリックでコマ送り...
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省略可...
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は、値を返してくれるが直接図化はされないため、一度d...
''追記''：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 zlo...
対策は2通りある。ひとつは散布図(set gxout scatter)を使う...
- 両対数グラフ／scatter使用時
散布図にしてXの値・Yの値にlogをかける。
set gxout scatter
set vrange X1 X2 ＃通常の線グラフならY軸に対応するが...
set vrange2 Y1 Y2
set digsiz 0.1 ＃デフォルトの点は若干小さいので大き...
d log10(lon);log10(var)
- 両対数グラフ／zlog使用時
描きたい図の横軸をZ軸にしたデータとコントロールファイルを...
スペクトルの表示に便利。
#ref(KE3_18_w.gif,/IT memo/GrADS memo,right,around,zoom,1...
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を対数表示。
*** 平均日変化グラフの描画（平均季節変化グラフも同様） [#...
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+...
＃ ミソはaveの設定。はじめに設定したt=1,4のそれぞ...
＃ で平均を取る。
あーラクチンだねー。
''追記''
上記は要はClimatologyを計算していることに他ならない。それ...
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というフ...
80行までの小さなASCIIファイルを読み込める。この機能を使う...
関連：[[テキストファイルの出力:http://meteora.ucsd.edu/~k...
*** 外部コマンドの呼び出し [#h0424e4e]
スクリプトやコマンドラインにおいて、びっくりマークから始...
tt=1
'! awk "NR=='tt+1'{printf(\"%4.1f\n\",\$3}" ./slpanm_eva...
** 上級編（UDF使ったり、station data使ったり） [#o4686eb8]
*** UDF設定方法 [#n8d866dd]
- UDFTファイルの場所を環境変数GAUDFTに記述
- 関数が増えるたびにUDFTファイルに追記していく。
*** 空間内挿（lterp） [#yed267f9]
- 便利な空間内挿関数。時間を動かせないのでそうしたい場合...
- Googleで探してもソースが出てこないので、ここに貼り付け...
- 使い方
d lterp(var1.1,var2.2) ＃ var1.1: 内挿前 var2.2: 内挿先
d lterp(var1.1,var2.2)-var2.2 ＃ などという演算が可能...
- 注意点
-- 鉛直内挿ではない。
-- 鉛直レベルが同じだと、鉛直方向にも同時作業してくれる。
-- 解像度が異なるデータで、ある点の時系列を同時に描画した...
*** 鉛直内挿（zinterp） [#if5146f2]
- 便利な鉛直内挿関数。GASCRP内になさければ、ftp://grads.i...
- 気圧面データから高度面の値を出したり、高度面データから...
- 使い方（スクリプトそのものに丁寧に書かれています）
d zinterp(field, zgrid, zlev)
例えばNICAMのモデル面（地表からの高さ）での出力データから...
open ms_pres.ctl
open topog.ctl
d zinterp(lev,ms_pres,50000)+topog.2
となる。850hPaの東西風は
open ms_pres.ctl
open ms_u.ctl
d zingerp(ms_u.2,ms_pres,85000)
- 注意点
-- zgridとzlevは一緒の単位を持つ。
-- 線形内挿のため、高度面データから気圧面を作る場合、（気...
*** 任意の断面図 [#ea4a6ac9]
- gr2stn, collect, coll2grを使用。が、鉛直分布「しか」描...
-スクリプト例
'reinit'
'open /home/kyoshimura/IsoRSM/cal10km/200702/r_pgb_prs.c...
'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,160x...
- 参考図（→右の図）
- やり方
set z 1 10
define tmp=pressfc ＃ 全層に表面気圧を入れる
などとやってから、 数多くの正値の太いコンターでd lev-tmp/...
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をうまく使って調整する。
- 例：&ref(CZC-ALT.gs);参照
*** Gradient/Laplacianの計算 [#cdaf01c9]
hdivgとhcurlをつかって、スカラー量qのdq/dx、dq/dy及びラプ...
- &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と...
+あとは、適宜それらのコントロールファイルを開き、図化すれ...
*** 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 (meridion...
- 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*co...
'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 3...
'rgbset2'
'd msf/1e9'
'set gxout contour'
'set clevs -300 -250 -200 -150 -120 -90 -60 -30 -10 10 3...
'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...
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 に基づいて相対湿...
''インストール''
+&ref(sp2rh.tar); をudfディレクトリ内で解凍し、
+適当なコンパイラでsp2rh.fをコンパイル（エンディアンには...
+udftにsp2rh.udftをコピペする（パス等は要修正）。
''使い方''
d sp2rh(temp,pressure,spec_hum)
# where temp[K], pressure[Pa], spec_hum[kg/kg] and RH[%]
# if you want to use constant value,
# use dummy variable as "pressfc-pressfc+100000". (simil...
***比湿の計算（自作UDF） [#e9dff544]
上記のsp2rhの逆。
''インストール''
+&ref(rh2sp.tar); をudfディレクトリ内で解凍し、
+適当なコンパイラでrh2sp.fをコンパイル（エンディアンには...
+udftにrh2sp.udftをコピペする（パス等は要修正）。
''使い方''
d rh2sp(temp,pressure,relative_hum)
# where temp[K], pressure[Pa], relative_hum[%] and spec_...
# if you want to use constant value,
# use dummy variable as "pressfc-pressfc+100000". (simil...
** OpenGrADS [#v81436b6]
- 以下は古い情報（ver1.9）
経緯は忘れたがOpenGrADSのページからダウンロード・インスト...
> q udct
> q udft
でずらずらっと出てきていたら、インストール済み。
-新しい情報（ver2.0以降）
--http://sourceforge.net/projects/opengrads/files/
のgrads2から自分の環境のものをダウンロードする。（僕の場...
--中にあるContents/gradsを実行する。
$ ./Contents/grads
環境変数を変えていないため、基本的に、昔のLibraryやUDFTも...
*** 流線関数・速度ポテンシャルの描画 [#ab437288]
要OpenGrADS。
- 流線関数
> d fish(hcurl(ugrdprs,vgrdprs))
or
> d fish_psi(ugrdprs,vgrdprs)
- 速度ポテンシャル
> d fish(hdivg(ugrdprs,vgrdprs))
or
> d fish_chi(ugrdprs,vgrdprs)
*** 球面調和関数のパワースペクトル [#s504239c]
要OpenGrADS
> run power.gs ugrdprs
とやるだけで、波数に対するKinetic Energyが出る。また、
> d sh_filt(ugrdprs,10)
で、波数10でのフィルタリング。
長年、これがほしかった。
*** 持ち上げ凝結高度（LCL） [#bda257a0]
要OpenGrADS
> rh2m=100*pressfc*spfh2m/(0.622+0.378*spfh2m)/satvap(tm...
> plcl=plcl(tmp2m,rh2m,pressfc) ; calculate LCL pressure
> d zinterp(hgtprs,lev,plcl/100) ; calculate height for ...
OpenGrADSを立ち上げると、これまでの自作のUDFが使えなくな...
** NetCDFファイルの読み込み [#k92a27aa]
*** COARDS規約のファイル [#if3f099e]
sdfopen XXX.nc
するだけ。
*** 読めそうだけどsdfopenが効かない場合 [#d30c4d73]
簡素版コントロールファイルをつくってxdfopenする。
xdfopen test.ctl
[[SWINGデータ:http://atoc.colorado.edu/~dcn/SWING/databas...
- templateオプションが効くが、たまにSegmentation Faultし...
-xdfopen用ctlファイルの作り方
DSET ^foo.nc # netcdfファイルを指定
UNDEF 1e20 # undef値を指定
XDEF xc # netcdfファイルのX軸変数を指定（lon, long...
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ファイルを開く。（デフォルト設定でWindowsの画...
-- 最後に、tiffファイルとして保存する。
- ただ、かなりむだに容量が大きくなるのが困り者。
[2009/5/16 更新]
-Powerpointで作った図のEPS化には、WMF2EPSというソフトが便...
--http://www.wmf2eps.de.vu/
--http://www.iplab.is.tsukuba.ac.jp/~zono/tips/wmf2eps_on...
--http://www.bunmeisha.co.jp/LaTeX2e/latex2e_eps.html
-難点
--拡張メタファイルではトリミングが無視される。そのため、...
** 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...
-http://www.atmos.ucla.edu/~munnich/grads/gxgif.html
で触れられているように、gxepsやgxpsで作った図には変な横線...
対策としては、
- shadedを使わない。contourかgrfillを使う。（grfillだと横...
- EPS/PSを使わない。gifなどで論文にしてしまう。GrADSのEPS...
-- この記述はあいまいで、printimを使って出力されたスクリ...
くらいかなあ。
''<追加情報>''
gimpでセーブしなおすと横線は消えるらしい。（by T.Yamada@G...
''<追加情報2>''
-http://www.kamiguchi.org/kenji/2006/03/grads-eps.html
の情報によると、アンチエイリアシングなるものの効果による...
-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が用意されている。...
われながら暇人だ。
- 現在使えるフォントの確認
font number
で一覧が表示される。デフォルトでは０－５。チルダ（～）が...
- フォントファイルの解読
http://iges.org/grads/gadoc/font.htmlを参照。
fontX.datには、ASCIIコードの32（空白）～126番目（チルダ）...
続いてフォントの中身を解読する。例えばスラッシュ（／：16...
3G][BIb
とある。まずはじめの3は、3つの座標ペアが続く、という意味...
はじめの数字の直後は、記号の幅を示す。つまりG]は-11～+11...
これ以降は、一筆書きなら座標を並べていき、間をおくなら「...
*** GrADSでパーミルを使う [#qe8ff029]
font3.datの最後の行（チルダ用）を以下のように変更する。
14F^IUISJPLONOOPPQTTVUXUZT[Q[O
↓
44D`WFE[ RJFLHLJKLIMGMEKEIFGHFJFLGOHRHUGWF RSTQUPWPYR[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 f...
14F^IUISJPLONOOPPQTTVUXUZT[Q[O
↓
44D`WFE[ RJFLHLJKLIMGMEKEIFGHFJFLGOHRHUGWF RSTQUPWPYR[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の...
- 一方、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-var...
は異なることがある。
ave(var1-var2,t=T1,t=T2) と等しくなるのは
ave(maskout(var1,var2+1e30,t=T1,t=T2)-ave(maskout(var2,v...
である。1e30というのは、var1のとりうる範囲が-1e30以上、...
Page:
