GrADSメモ

個人的なメモ書きです。別の(楽な)方法があればぜひ教えてください。誤りなどの指摘もよろしくお願いします。

できたらいいな。

ご存知のことがあれば、コメントしてくださいm(_ _)m

  • ある2つの面的な時系列データの、同じグリッド同士の時系列データの相関係数分布の表示。
  • ある面的データの、球面調和関数スペクトルの算出
  • [Done]テキスト出力のファイルへの保存。 -- kei 2007-08-01 (Wed) 13:33:06
    • 自作した。./textout
    • UDFを使うともっと気の利いたものが作れそう。たとえばここのような。
  • きれいな画像の保存(printimの解像度はずいぶん低いので…) -- kei 2007-08-01 (Wed) 13:34:39
    • ここで書かれているように、epsとして保存するのがベター。 -- 2007-11-07 (Wed) 10:50:21
  • [Done]鉛直断面図に、地形マスクをかける。 -- 2007-08-15 (Wed) 14:55:58
    • できたっぽい。ここ参照 -- 2007-11-07 (Wed) 02:17:21
  • 気圧面データから高度毎の分布を描く。 -- kei 2007-11-07 (Wed) 02:37:56

#comment2_kcaptcha

初級編

入手とインストール

  • 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

基本コマンド

  • open controlfile
    • netCDFの場合は sdfopen ncfile
  • d expr
  • set x (y, z, t) num1 num2
  • set lon (lat, lev) val1 val2

出力形式

set gxout outputtype
  • よく使うのは、contour, grfill, shaded, vector, scatter 等
  • ベクトル表示
    d u;v
    d skip(u,2,2);skip(v,2,2) # 間引き

演算

  • 変数同士の四則計算
    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ファイルからデフォルトを選ぶ

set dfile num

ある値以下(以上)をマスクして表示.

d maskout(var,var-10) # 10以上の値を表示

背景を白に

set display color white

中間の値に色を塗らない

set black -0.2 0.2 # -0.2から0.2は塗らない
  • 背景色(色番号1)を塗っているので透過はしないことに注意。

カラーバー表示

run cbar # fg,bgなし
run cbarn # fg,bgあり
  • オプション
    run cbarn 大きさ 水平0/鉛直1 x位置 y位置

タイトル、軸ラベルをつける

draw title xxxxx
draw xlab (ylab) xxxxx

プリントアウト

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

バイナリで保存。(データの切り出しなど。)

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

中級編(scriptを使ったり、デフォルト以外の関数を使ったり)

scriptで引数を使う

function main(args)をスクリプトの頭で定義する。

(example)
function main(args)
var1=subwrd(args,1)
var2=subwrd(args,2) 
...

テキスト出力

set gxout print
set prnopts %10.3e 1 1 # 小数点以下3桁、一列に表示、ブランクは一文字分

画面に出たものをファイルに保存するやり方がわからん!コピペするしかないのか??スクリプトなら何とかなりそうだけど。。。

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

(線の)凡例をつける。

  • 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標準は下記の通り。
    NumColorLineMark
    1112
    2313
    3714
    4215
    5611
    6912
    71013
    81114
    91215
    101511

アニメーション

普通に

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回繰り返し

など。オプションは変数より前にないとダメらしい。

画面内の分割

page.gs (ftp://ftp.cpc.ncep.noaa.gov/wd51we/grads/page.gs)が使える。

page q1 # 右上
d var1
page q2 # 左上
d var2

という具合。元に戻すには

page reset

移動平均の図化

set t 1 100               # 例えば。
define a=ave(var,t-0,t+9) # 前方10個平均。define省略可。第二引数はtだけだとダメ。
set missconn on           # 欠損は無視してつなぐ。
d skip(a,10)

領域平均値の時系列グラフ

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)

統計情報の表示

set gxout stat
d var
 欠損値の数や和、標準偏差、分散などが表示される。
 sublin, subwrdなどで取り出して使用することが多い。

相関係数の算出(時間変動)

ある時間軸一次元のデータの、領域全体の時系列変化に対する相関係数を計算する。

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の変動に対しての、全グリッドの時間変動の相関を表示

相関係数の算出(空間変動)

同領域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)

両対数グラフの作成

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軸にしたりする。) スペクトルの表示に便利。
    KE3_18_w.gif
    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+120,4))
   # ミソはaveの設定。はじめに設定したt=1,4のそれぞれから4つ飛び
   # で平均を取る。

あーラクチンだねー。

グラフの中にグラフを書く

vpageを使う。ただ、位置などの微調整が難しい。

KE_9_2.gif
# 現在の座標情報を得る(あくまでも目安)
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以前の設定は持ち越されないので注意。

上級編(UDF使ったり、station data使ったり)

UDF設定方法

  • UDFTファイルの場所を環境変数GAUDFTに記述
  • 関数が増えるたびにUDFTファイルに追記していく。

空間内挿(lterp)

  • 便利な空間内挿関数。時間を動かせないのでそうしたい場合は少し面倒。もともとetaについてくる関数っぽい?
  • Googleで探してもソースが出てこないので、ここに貼り付けてしまう。オープンソースでしょう。たぶん。&attachref(lterp.tar);
  • 使い方
    d lterp(var1.1,var2.2) # var1.1: 内挿前 var2.2: 内挿先
    d lterp(var1.1,var2.2)-var2.2 # などという演算が可能となる。楽。
  • 注意点
    • 鉛直内挿ではない。
    • 鉛直レベルが同じだと、鉛直方向にも同時作業してくれる。
    • 解像度が異なるデータで、ある点の時系列を同時に描画したいときには、前もってdefineする必要がある。

任意の断面図

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

鉛直分布図などに地形マスクをかける

x-sec_06Z13FEB2007_38.571.gif
  • 参考図(→右の図)
  • やり方
    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 でマスクしようとしても、「透過色」を指定できないのでうまくいかない。

鉛直分布図と非鉛直分布図を同時に描く。

  • 参考図(→上の図)
  • やり方:vpageをうまく使って調整する。
  • 例:&attachref(CZC-ALT.gs);参照

Gradient/Laplacianの計算 NEW!

hdivgとhcurlをつかって、スカラー量qのdq/dx、dq/dy及びラプラシアン(d2q/dx2+d2q/dy2)を計算できる。簡単にできるスクリプトはこちら。

  • &attachref(lap.gs);
  • &attachref(gradx.gs);
  • &attachref(grady.gs);

NetCDFファイルの読み込み

COARDS規約のファイル

sdfopen XXX.nc

するだけ。

読めそうだけどsdfopenが効かない場合

簡素版コントロールファイルをつくってxdfopenする。

xdfopen test.ctl

SWINGデータ用のサンプル→&attachref(s1a.echam4.ctl);

  • templateオプションが効くが、たまにSegmentation Faultしてしまう。(特に、演算による時系列アニメで、ファイルをまたぐ場合)

Office使いのための論文用画質維持について

  • 複数の図を一枚のパネルにする際、パワポで切り貼りしてワードに貼り付けるという手法は、論文だけでなくそのまま発表にも使えて便利。
  • が、論文最終段階で、そのパネルを画像(Tiff or eps)として提出する場合、画質がかなり落ちる。(画面上の画質と同じになる(?))
  • 対処法としては、
    • パワポで、複数の図からなるパネル(背景含む)を全部選択し、右クリックで「画像として保存」を選択、「拡張メタファイル(.emf)」で保存。
    • 次に、emfファイルを開く。(デフォルト設定でWindowsの画像Viewerが起動する)
    • 最後に、tiffファイルとして保存する。
  • ただ、かなりむだに容量が大きくなるのが困り者。

EPS/PSで保存した場合の横線の謎 NEW!

tmp.gif

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

対策としては、

  • shadedを使わない。contourかgrfillを使う。(grfillだと横線だけでなく格子線になるが、比較的目立たない。)
  • EPS/PSを使わない。gifなどで論文にしてしまう。GrADSのEPSファイルって、結局ほぼラスターデータだから(文字やshadeが相当小さな単位でのベクトルデータ)、容量や画質はそんなに変わらん。

くらいかなあ。

フォントの作り方

デフォルトでlib内にfont0.dat~font5.datが用意されている。font3.datが特殊文字。しかし、同位体屋に必須である「パーミル(‰)」が存在しない。よって自作を試みる。 われながら暇人だ。

  • 現在使えるフォントの確認
    font number
    で一覧が表示される。デフォルトでは0-5。チルダ(~)がどのフォントでもまったく同じ文字に当てられている。
  • フォントファイルの解読

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でパーミルを使う

font3.datの最後の行(チルダ用)を以下のように変更する。

 14F^IUISJPLONOOPPQTTVUXUZT[Q[O
  ↓
 44D`WFE[ RJFLHLJKLIMGMEKEIFGHFJFLGOHRHUGWF RSTQUPWPYR[T[VZWXWVUTST R[TYUXWXYZ[\[^Z_X_V]T[T

GrADS内では、「`3~」で呼び出す。例えば、

ga> draw ylab `3d`0`a18`nO [`3~`0]

で、『δ18O [‰]』(18は上付き)が表示される。

Using Permil (part per thousand) in GrADS

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を使ったフォーマット変換

別項目へIT memo/GrADS memo/lats4d

注意点アレコレ

  • 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以上、ってことにしているから。