GrADSメモ

Script関係

引数を使う

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

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

注意点アレコレ

  • 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の結果を欠損値を除いた個数で割る。

気圧重みづけ鉛直積分

d vint(start pressure (hPa),var,end pressure (hPa))

複数のctlファイルを開けて,それを選ぶ

> set dfile 番号

ある値以下{以上}をカットして表示.

> d maskout(var,{-}var {値})

シェード{コンター}で表示

> set gxout shaded {contour}

平均

> ave(var,t=1,t=24)など.ネスティングも可能

演算

  • 変数同士の四則計算,ave,sum,log,pow(x,y)なんかが使用可能.

カラーバー表示

> run cbar --> fg,bgなし
> run cbarn --> fg,bgあり

ヴェクター表示。

> d u;v
間引きしたいときは
> d skip(u,2,2);skip(v,2,2)

プリントアウト

> enable print hoge.gx
> print
> disable print
> ! gxps -c -i hoge.gx -o hoge.ps
-c -rで黒背景。

追記:printimというコマンドを使えば一発。

printim tmp.gif gif white

テキスト出力

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

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

タイトルつける

> draw title xxxxxxxxxxx

netCDFを読む

% gradsnc
> sdfopen xxx.nc
> q file など

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

> 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

収束->ラクチン!

> hdivg(x方向フラックス,y方向フラックス)
鉛直積分値があればそのまま使えるし、下のvintを使って各高度の値から鉛直積分を求めても使える。

可降水量

> vint(地表気圧(mb単位), 変数(比湿等), 上がどこまで(定数))
多分、ある層からある層までを区切る、なんてことはできない!?
気圧を重力加速度で割る、という作業はすでにされているようだ。

背景を白に

> set display color white

中間の値に色を塗らない

> set black -0.2 0.2 (-0.2から0.2は塗らない)

(線の)凡例をつける。

  • 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軸にしたりする。) スペクトルの表示に便利。

#ref(): File not found: "KE3_11_w.gif" at page "IT memo/GrADS memo"

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を対数表示。