GrADSのTips
okanishi- お知らせ:
- 初めてGrADSを使う人向けの内容ではありませんので、初学者は他の方が作ったページを参照してください。
- スクリプトについてはGrADSスクリプトのTipsをご覧下さい。
- このTipsは本講座メンバーをはじめ様々な方々の協力・助言・激励によって作成されたものです。
- 以前のページはこちら
目次
- 凡例
- 起動と終了
- ファイル操作
- 状態表示（querry）
- 描画一般
- 図を描画する
- 表示されている図を消去する
- ファイル番号を指定して描画する
- デフォルトのファイル番号を変更する
- 各次元をステップ数で設定する
- 各次元を座標値（world coordinate）で設定する
- 様々な次元でアニメーションを見る
- 次元が低い場合でも強制的にアニメーションを見る / 解除する
- 任意の断面を切り出す
- アニメーションをゆっくり見る
- タイムステップを指定して描画する
- GrADSのロゴ（左下）と日付（右下）を表示する／しない
- 日付（右下）を表示する／しない
- 描画の度に毎回実行するスクリプトを指定する／指定を解除する。
- 背景色を変更する
- Gradsウインドウに他のウインドウを重ねても消えないようにする
- 描画位置を設定する
- 領域全体が未定義値の場合、「Entire Grid Undefined」と表示する / しない
- 図の外側に枠を表示する / しない / 丸型にする
- 地図と投影法
- 座標軸とラベル・補助線
- 等値線（contour）
- ベクトル（vector）
- 折線グラフ (line)
- 折れ線の塗りつぶし（linefill）
- バー（bar）
- 散布図（scatter）
- 文字
- 図形
- コントロールファイルの文法
- 緯度経度座標以外のデータの読み込み（PDEF）
- 変数
- 関数
- ユーザ定義関数
- その他
- 番号一覧
凡例
- 特に断りがない限り、以下のTipsはLinuxを想定しています。
- "ga->" はGrADS中のコマンドライン端末（"ga->" を省略する場合も多い）
- "[mydir@hosts]" はLinuxのコマンドライン端末
- コマンド説明の便宜上、"//～(説明文)～" と表記する場合があります。コマンド入力の際は、//から文末までは無視してください。
- "/"は又は、"[]"は省略可能。
起動と終了
GrADSを起動する
[mydir@hosts] grads Grid Analysis and Display System (GrADS) Version 1.8SL11 Copyright (c) 1988-2001 by Brian Doty Center for Ocean-Land-Atmosphere Studies Institute for Global Environment and Society All Rights Reserved > Config: v1.8SL11 32-bit little-endian readline printim Issue 'q config' command for more information. Landscape mode? (no for portrait): GX Package Initialization: Size = 11 8.5 ga->
"Landscape mode?"の部分で何も入力せずにEnterキーを押すと、landscape(横長)になります。 ここでnoと入力してEnterキーを押すと、portrait(縦長)になります。landscape、portraitを指定して起動する場合は、
[mydir@hosts] grads -l // landscape [mydir@hosts] grads -p // portrait
でもOKです。
なお、Version 1.9 以下のGrADSでは、gradsの代わりにgradscと入力する必要があります。これが面倒でしたら、予めリンクやエイリアスを設定するとよいでしょう。
非X環境で起動する
[mydir@hosts] grads -b
引数-bを付けて起動することで、グラフィックス表示以外は通常通り使用することができます。windowsからTeraTermで接続してGrADSを使う場合に役立ちます。描画処理がない分、処理が高速になることも期待できます。
起動直後に実行するコマンドを指定する
grads -c コマンド
- (例)GrADS起動直後にhelpを実行する。
[mydir@hosts] grads -c help
GrADSを起動したまま親シェルのコマンドを実行する
!コマンド
- (例)GrADSを起動したままemacsでu.ctlを編集する
ga-> !emacs u.ctl & ga-> (emacs実行中もGrADSはそのまま使える)
GrADSを終了する
quit
ファイル操作
コントロールファイルをオープンする
open ファイル名
- (例)test.ctlというコントロールファイル(データの情報が入ったファイル)を開く。
ga-> open test.ctl Scanning description file: test.ctl Data file test.dr is open as file 1 LON set to 0 360 LAT set to -87.864 87.864 LEV set to 1000 1000 Time values set: 1998:8:1:0 1998:8:1:0 ga->
コントロールファイルをオープンすると番号が振られます(上の例では1)。複数のファイルを同時に開いた場合は2,3,...と付番され、複数ファイルの区別に使われます。蛇足ですが、末尾の.ctlは省略できます(ただし、testというファイルが別に存在すればエラーになる)。
- (例)test.ctlを開く(その2)。
ga-> open test Scanning description file: test Data file test.dr is open as file 1 LON set to 0 360 LAT set to -87.864 87.864 LEV set to 1000 1000 Time values set: 1998:8:1:0 1998:8:1:0 ga->
コントロールファイルをクローズする
close ファイル番号
ファイル番号はオープン時に付番されます。
- (例1)ファイル番号1のファイルを閉じる
ga-> close 1
複数ファイルを開いている場合、ファイル番号の一番大きなファイルしか閉じることができません。
- (例2:誤)2つのファイル(ファイル番号1,2)を閉じる
ga-> close 1 Close Error: Only last file may be closed
- (例2:正)2つのファイル(ファイル番号1,2)を閉じる
ga-> close 2 File 2 has been closed ga-> close 1 File 1 has been closed
netcdf(*.nc)をオープンする
sdfopen ファイル名
- (例)slp.1995.nc(NCEPの海面較正気圧)を開く
ga-> sdfopen slp.1995.nc
Version 1.9以下のGrADSを使用している場合、GrADS起動時のコマンドは、gradscではなくgradsncとする必要があります。Version 2.0系ではgradsでOKです。
蛇足ですが、netcdfのファイルの情報を見たい場合、Linuxのコマンドラインで
[mydir@hosts] ncdump -c ファイル名
とすると見られます。
epsとして保存する
enable print ファイル名(.gmf) print disable print !gxeps -c -i ファイル名(.gmf) -o ファイル名(.eps)
1～3行目ででgmfファイル(GrADS Metacode Format)を作り、4行目でepsに変換しています。
- (例)現在表示されている図をtest.epsとして保存
ga-> enable print test.gmf ga-> print ga-> disable print ga-> !gxeps -c -i test.gmf -o test.eps
netCDFファイルを開くときの注意
64-bit機のGrADSの場合、「sdfopen」で開くと、まれに「289.615」という値を持つ箇所が、欠損値または「-32767.0」として認識される。
32-bit機のGrADSで開くか、Ctlファイルを作ってnetCDFを開くようにすると、解決できる。
状態表示（querry）
ファイルの情報を表示する
q file ファイル番号
データファイル名や次元、変数などの情報が表示されます。
- (例)ファイル番号1番のファイルの情報を表示する
ga-> q file 1 File 1 : TEST DATA Descriptor: test.ctl Binary: test.dr Type = Gridded Xsize = 128 Ysize = 64 Zsize = 16 Tsize = 3100 Number of Variables = 3 u 16 0 ZONAL VELOCITY v 16 0 MERIDIONAL VELOCITY t 16 0 TEMPERATURE
なお、ファイル番号が1の場合は番号を省略できます(q file のみでOK)。
コントロールファイルの中身を表示する
q ctlinfo ファイル番号
「q file」では整形されて表示されますが、「q ctlinfo」は整形されずにそのまま全部表示されるようです。
現在の次元を表示する
q dims
- (例)
ga-> q dims Default file number is: 1 X is fixed Lon = 0 X = 1 Y is varying Lat = -90 to 90 Y = 0.228044 to 64.772 Z is varying Lev = 1000 to 10 Z = 1 to 19 T is fixed Time = 00Z01AUG1998 T = 1
現在のグラフィック環境の情報を表示する
q gxinfo
直前に描いた図の種類やページのサイズなどを取得できます。
- (例)
ga-> q gxinfo Last Graphic = Contour Page Size = 11 by 8.5 X Limits = 4.4 to 7.4 Y Limits = 1 to 4 Xaxis = Lat Yaxis = Lev Mproj = 2
シェードの色を取得する
q shades
- (例)
ga-> q shades Number of levels = 6 16 < -10 17 -10 -5 18 -5 0 19 0 5 20 5 10 21 10 >
左側から、色番号、下端値、上端値を表しています。
座標の相互変換を行う
q xy2w v1 v2 // (X,Y)系 → world系 q xy2gr v1 v2 // (X,Y)系 → (gX, gY, gZ, gT)系 q w2xy v1 v2 // world系 → (X,Y)系 q w2gr v1 v2 // world系 → (gX, gY, gZ, gT)系 q gr2w v1 v2 // (gX, gY, gZ, gT)系 → world系 q gr2xy v1 v2 // (gX, gY, gZ, gT)系 → (X,Y)系
- (X,Y)系 : display上の座標。portraitの場合、0<=X<=11、0<=Y<=8.5。
- world系 : データ自身の座標。例えばlat-lon断面を描いている場合、(lat,lon)系。
- (gX,gY,gZ,gT)系 : 格子番号で表した座標。
現在定義されている変数名を得る
q define
現在の未定義値を得る (ver. 2.0.a6以降)
q undef
描画一般
図を描画する
display 変数名 又は d 変数名
dはdisplayの省略形です。描画のされ方(2D,3D,アニメーション etc)は、次元やgxoutの設定に依存します。
- (例)変数uを描画する
ga-> d u
表示されている図を消去する
clear 又は c
cはclearの省略形です。
clearを行うと、パラメータの多く(例えば色)も同時に初期化されます。これを回避するには、以下のようにnorsetオプションを使います。
c norset
ファイル番号を指定して描画する
d 変数名.ファイル番号
open時に割り振られたファイル番号を使うと、複数のファイルの変数を同時に扱うことができます。
- (例)ファイル番号1の変数Uと、ファイル番号2の変数Vを同時に図示する
ga-> d U.1 ga-> d V.2
- (例)ファイル番号1の変数Uと、ファイル番号2の変数Uの差を図示する
ga-> d U.2-U.1
デフォルトのファイル番号を変更する
set dfile ファイル番号
ファイル番号を指定せずにdrawした場合のデフォルトを指定します。
- (例)ファイル番号1の変数Uと、ファイル番号2の変数Vを同時に図示する
ga-> set dfile 1 // デフォルトのファイルを1に指定 ga-> d u // Uを描画 ga-> set dfile 2 // デフォルトのファイルを2に指定 ga-> d v // Vを描画
各次元をステップ数で設定する
set x 開始ステップ [終了ステップ] // X（経度）次元 set y 開始ステップ [終了ステップ] // Y（緯度）次元 set z 開始ステップ [終了ステップ] // Z（高度）次元 set t 開始ステップ [終了ステップ] // T（時間）次元 set e 開始ステップ [終了ステップ] // E（アンサンブル）次元
終了ステップは省略可能です。省略すると固定（fixed）次元、省略せず「開始≠終了」の場合は変動（varying）次元となります。
- (例1) Z=3におけるUの緯度経度断面(1<=x<=15, 1<=y<=20)を作成する
set x 1 15 set y 1 20 set z 3 d U
- (例2) X=20におけるTの緯度高度断面(1<=y<=50, 1<=z<=30)を作成する
set x 20 set y 1 50 set z 1 30 d T
- (例3) MSLの緯度経度断面(10<=x<=25, 10<=y<=25)のアニメーション(1<=t<=10)を見る
set x 10 25 set y 10 25 set t 1 10 d MSL
各次元を座標値（world coordinate）で設定する
set lon 開始座標値 [終了座標値] // 経度（X）次元 set lat 開始座標値 [終了座標値] // 緯度（Y）次元 set lev 開始座標値 [終了座標値] // 高度（Z）次元 set time 開始座標値 [終了座標値] // 時間（T）次元 set ens 開始座標値 [終了座標値] // アンサンブル（ENS）次元
使い方はほとんど「set x ～」等と同様です。但し、座標値はコントロールファイルの書式に従います。
- (例1) 100hPaにおけるU*Vの緯度経度断面(0<=lon<=360, -90<=lat<=90)を作成する
set lon 0 360 set lat -90 90 set lev 1000 d U*V
- (例2) 2000年における北緯40度、東経135度のSSTの移り変わりを見る
set lon 135 set lat 40 set time 01jan2000 31dec2000 d SST
様々な次元でアニメーションを見る
set loopdim x/y/z/t
- (例)x-y図でz=1からz=20までの各層をアニメーションで見る
set x 1 100 set y 1 100 set z 1 20 set t 1 set loopdim z d u
- (例)x-t図でz=1からz=20までの各層をアニメーションで見る
set x 1 100 set y 1 set z 1 20 set t 1 124 set loopdim z d u
次元が低い場合でも強制的にアニメーションを見る / 解除する
set looping on/off
折れ線グラフのアニメーションを見たい場合等に有効です。
任意の断面を切り出す
Gradsの本家にスクリプトが載っています。http://www.iges.org/grads/gadoc/usingstationdata.html#xsection
ただし、undefを含むデータを描く場合、Version2.0以降ではundef領域が塗りつぶされるバグがあります。注意してください。
アニメーションをゆっくり見る
ftp://grads.iges.org/grads/scripts/xanim.gs
以上のスクリプトを用いる。
まず、時間の指定。
set t 1 10
あとは、3通りほど書き方がある。
xanim -pause -grfill var ＃ マウスクリックでコマ送り（かつgrfillで表示） xainm -sec 1 var ＃ １コマ１秒ずつ表示 xanim -skip 4 -repeat 5 var ＃ ４コマ目ごと、5回繰り返し
早すぎて見えないときに便利。
参考:http://meteora.ucsd.edu/~kyoshimura/?IT%20memo%2FGrADS%20memo
タイムステップを指定して描画する
d 変数名(t=タイムステップ)
「set t」を使わずにタイムステップを指定して描画することができます。
- (例)変数Uのt=1～2の変化を描画する
ga-> d U(t=2)-U(t=1)
GrADSのロゴ（左下）と日付（右下）を表示する／しない
set grads on/off
日付（右下）を表示する／しない
set timelab on/off
描画の度に毎回実行するスクリプトを指定する／指定を解除する。
set imprun スクリプト名/off
背景色を変更する
set background 色番号
- (例)背景色を赤に変更する(そんな人いない？)
ga-> set background 2
Gradsウインドウに他のウインドウを重ねても消えないようにする
- grads内で設定するのではないのだが。
- 「/etc/X11/xorg.conf」の「Section "Screen"」に以下を追加する(Vine4.0 or Fedora5) Option "BackingStore" "Yes"
- 「/etc/X11/xorg.conf」の「Section "Device"」に以下を追加する(Vine3.2) Option "BACKING STORE" "Yes"
- 「/etc/X11/xorg.cnf」がない場合 (CentOS 6.2) init 3 //これでテキストログイン Xorg -configure // /root/xorg.conf.new ができる mv xorg.conf.new /etc/X11/xorg.cnf (「/etc/X11/xorg.cnf」がある場合と同様に編集) init 5 //これでGUIに戻る。もしかしたら再起動が必要かも。
描画位置を設定する
set parea xmin xmax ymin ymax set vpage xmin xmax ymin ymax
pareaは図の描画位置を設定します。ラベルやタイトルなどは含まれません。縦横比が決まっている場合（例：水平図）はほとんどの場合、設定した領域全体に描画が行われないので注意してください。
vpageは全体の描画位置を設定します。pareaと違いラベルやタイトルも含みます。
本家のマニュアルには、複数の図を同時に出力する場合、pareaよりもvpageの方が適切であると記述されています（一応どちらでもうまくいくと思いますが）。
領域全体が未定義値の場合、「Entire Grid Undefined」と表示する / しない
set datawarn on/off
図の外側に枠を表示する / しない / 丸型にする
set frame on/off/circle
circleは北極中心図などのときに有効です。
地図と投影法
以下の描画例はNCEP/NCAR再解析 を使用しています。
背景に地図を描く／描かない
set mpdraw on/off
地図のみを表示する
draw map
地図の色・スタイル・太さを設定する／自動設定する
set map 色番号 線種番号 太さ set map auto
地図の解像度を変える
set mpdset lowres/mres/hires
地図として政治境界を描く／描かない
set poli on/off
地図の解像度がmres/hiresの場合のみ有効です。デフォルトはon。
地図の様々な種類の線を個別に設定する
set mpt 種類番号 色番号 線種番号 太さ // 個別に指定 set mpt 種類番号 -1 // デフォルトの色・線種・太さ set mpt 種類番号 off // 表示しない
種類番号は海岸線や国境、州境などを表す番号です(何が何に対応するかはmapファイル依存)。種類番号に*を指定すると「全て」の意味になります。
- (例)海岸線のみ表示する
set mpt * off set mpt 0 -1 draw map
メルカトル図法で描く
set mproj latlon
デフォルトの図法です。
- (例) 変数Tをメルカトル図法で描く
set mproj latlon set lat -90 90 set lon 0 360 d T
画面全体に直交座標で描く（縦横比を無視する）
set mproj scaled set mproj off
scaledは地図あり、offは地図なし（緯度経度のラベルもなし）です。
- (例) 変数Tを縦横比を無視した直交座標で描く
set mproj scaled set lat -90 90 set lon 0 360 d T
ロビンソン図法で描く
set mproj robinson
lat, lon はそれぞれ[-90:90], [-180:180]のみ対応しています。
- (例) 変数Tをロビンソン図法で描く
set mproj robinson set lat -90 90 set lon -180 180 d T
モルワイデ図法で描く
set mproj mollweide
- (例) 変数Tをモルワイデ図法で描く
set mproj mollweide set lat -90 90 set lon 0 360 d T
正投影図法で描く
set mproj orthogr
「set lat -90 90」およびlonの範囲を180度にする必要があります。
- (例) 正投影図法で日本中心のTを描く
set mproj orthogr set lat -90 90 set lon 45 225 d T
- (応用例)地球はまわる
- 重い (1.6MB) ので要注意: orthogr_anim.gif
- GrADSスクリプト、およびconvertコマンドを使って作成。 orthogr_anim.pdf
北極(南極)中心の図を描く
- 北極中心
set mproj nps
- 南極中心
set mproj sps
npsの場合、緯度の設定は0～90Nの範囲で設定するのが無難でしょう。あまり範囲を広げすぎると、GrADSが怒ってデフォルトの図法で描いてしまうことがあります。
- 地図の回転 北極を中心にポーラーステレオで描く場合、日本を下側にするか上側にするかを変えたい時があります。それには、
set lon 0 360 set lon -180 180
などとすると、中心の経度を動かすことが出来ます。
- (例1) 変数Tを北極中心に描く
set mproj nps set lat 0 90 set lon 0 360 d T
- (例2) 変数Tを南極中心に描く
set mproj sps set lat -90 0 set lon 0 360 d T
- (例3) 日本中心で±90度のTを描く
set mproj nps set lat 0 90 set lon 45 225 d T
緯度(or 経度)鉛直断面で緯度経度を表すN, S, E, Wを消す
set mproj off
座標軸とラベル・補助線
鉛直軸を対数軸にする／解除する
set zlog on/off
座標軸のラベル・補助線を自由な位置に設定する
set xlevs レベル1 レベル2 ... set ylevs レベル1 レベル2 ...
- (例)y軸を設定する
ga-> set ylevs 1000 500 300 100 50 30 10 5 3 1
(注意)指定可能なレベルの数は49個までのようです。
ラベル・補助線の間隔を設定する
set xlint 間隔 set ylint 間隔
等間隔にラベル・補助線をつけることができます。
ラベルとして任意の文字列を設定する
set xlabs ラベル1 | ラベル2 | ... set ylabs ラベル1 | ラベル2 | ...
ラベル・補助線は等間隔に設定されます。
- (例)風速の頻度分布を描く
ga-> set xlabs 0-10[m/s] | 10-20[m/s] | 20-30[m/s] ga-> d u_freq
縦軸と横軸を入れ換える／元に戻す
set xyrev on/off
座標軸の向きを逆にする
set xflip on/off set yflip on/off
ふつうは右 or 上に行くほど座標の値が大きくなりますが、これをonにすることで右 or 上に行くほど座標の値を小さくすることができます。たとえばホフメラー図を描く場合、時間経過を「上→下」「下→上」どっちで表すか、これで指定できます。
ラベル・補助線を表示する／しない
set xlab on/off set ylab on/off
- (例)X軸にラベルをつけるが、Y軸にはラベルをつけない
ga-> set xlab on ga-> set ylab off
ラベルの書式指定する
set xlab format set ylab format
formatはC言語のprintfとほぼ同じ形式です（以下の例を参照）。
- (例)Y軸のラベルを小数指数形式にする
ga-> set ylab %f
- (例)Y軸のラベルを小数点以下2桁までの指数形式にする
ga-> set ylab %.2e
- (例)Y軸のラベルを指数部3桁、小数点以下2桁までの指数形式にする
ga-> set ylab %3.2e
- (例)X軸のラベルを対数グラフ的にする
ga-> set xlab 10`a%.1f
ラベルのフォントを設定する
set xlopts 色番号 [太さ [大きさ]] set ylopts 色番号 [太さ [大きさ]]
ラベルに用いるフォントの色、太さ、大きさを設定します。デフォルト値は色番号=1、太さ=4、大きさ=0.12です。xloptsはX軸(横軸)、yloptsはY軸(縦軸)のラベルを設定します。
- (例)x軸のラベル：色1、太さ2、大きさ0.2、を設定する
ga-> set xlopts 1 2 0.2
- (例)y軸のラベル：色1、太さ1、大きさ0.2、を設定する
ga-> set ylopts 1 1 0.2
x軸のラベルを描く場所を指定する
set xlpos 位置 側
「位置」にはx軸からのずれ、「側」にはt(上側)、b(下側)のどちらかを指定します。
- (例)x軸のラベルを、下側の通常より0.5下に描画する
ga-> set xlpos -0.5 b
y軸のラベルを描く場所を指定する
set ylpos 位置 側
「位置」にはy軸からのずれ、「側」にはr(右側)、l(左側)のどちらかを指定します。
- (例)y軸のラベルを、左側の通常より0.5左に描画する
ga-> set ylpos -0.5 l
補助線のみ表示し、ラベルを表示しない
経験則ですが、
ga-> set xlopts 1 0 0 // x軸 ga-> set ylopts 1 0 0 // y軸
で、補助線のみを表示することができます。
線種・色を指定して補助線を表示する／補助線を表示しない
set grid on 線種 色 // 水平、鉛直 set grid horizontal 線種 色 // 水平のみ set grid vertical 線種 色 // 鉛直のみ set grid off // 補助線を表示しない
特定の座標軸(緯線・経線)の書き方
- (簡便な方法) 日本付近の画像で北緯35度線と東経135度線を描く。
ga-> set lat 20 50 ga-> set lon 120 150 ga-> set cint 35.0 ga-> d lat ga-> set cint 135.0 ga-> d lon)
- (より完全な方法) 北緯35度線と東経135度線を描く。
ga-> set cint 35.0 ga-> d maskout(lat,(34.0-lat)*(lat-36.0)) ga-> set cint 135.0 ga-> d maskout(lon,(134.0-lon)*(lon-136.0))
時間軸のラベルに年／月日を入れない
set tlsupp year set tlsupp month
yearを指定すると年、monthを指定すると年と月日の表示が消えます。yearは気候値を描く時に便利です。この設定は"c"を行なうと解除されます。
等値線（contour）
等値線を引く
set gxout contour d 変数
- (例)変数Uの等値線を引く
ga-> set gxout contour ga-> d U
等値線を引く値を指定する
set clevs レベル1 レベル2 ...
- (例)複数の値を指定して変数Uの等値線を引く
ga-> set clevs -5 -2 0 2 5 ga-> d U
- (例)0線を引く
ga-> set clevs 0 ga-> d U
等値線の最大値、最小値、間隔を指定する
set cmax 最大値 set cmin 最小値 set cint 間隔
- (例)変数Uについて、0～20の範囲で1毎に等値線を引く。
ga-> set cmin 0 ga-> set cmax 20 ga-> set cint 1 ga-> d U
等値線を描かない範囲を指定する
set black 最小値 最大値
- (例)変数Uについて、0～20の範囲で1毎に等値線を引く。但し-5～5の範囲は描かない。
ga-> set cmin 0 ga-> set cmax 20 ga-> set cint 1 ga-> set black -5 5 ga-> d U
等値線の色(一色)を設定する
set ccolor 色番号
- (例)変数Uの等値線を黒で描く
ga-> set ccolor 1 ga-> d U
等値線の色(複数)を設定する
set ccols 色番号1 色番号2 ...
普通はset clevsで等値線の値を設定してから使います。
- (例)等値線を複数の色で描く
ga-> set clevs 0 5 10 15 20 ga-> set ccols 1 2 3 4 5 ga-> d U
- 設定した色の情報は、一回のdに対してのみ有効です。また cコマンドを実行すると、設定した内容も消えてしまいます。
- 色番号一覧
等値線の太さを変える
set cthick 太さ
太さは1～20まで設定可能です。
- (例)等値線の太さを4に設定する
ga-> set cthick 4
等値線の線種を設定する
set cstyle 線種番号
等値線にラベルをつける／つけない／強制的につける
set clab on/off/forced
onとforcedの違いは、onは全ての等値線にラベルがつくとは限らないのに対して、forcedは全ての等値線に強制的にラベルを付けます。
等値線のラベルのフォントを設定する
set clopts 色番号 [太さ [大きさ]]
等値線のラベルに用いるフォントの色、太さ、大きさを設定します。デフォルト値は色番号=1、太さ=1、大きさ=0.09です。
- (例)等値線のラベル：色1、太さ2、大きさ0.2、を設定する
ga-> set clopts 1 2 0.2
場所によって等値線の間隔を変える
場所によって値が急激に変化する場合に有効です。たとえば、ある経度における物理量Aの鉛直断面を書く場合、
ga-> set lon 120 //東経120°で断面を切る ga-> set lat -90 90 ga-> set t 1 ga-> set lev 1000 100 ga-> A1=A //1000～100hPaのデータを一時的に格納 ga-> set lev 100 10 ga-> A2=A // 100～ 10hPaのデータを一時的に格納 ga-> set lev 1000 10 ga-> set cmin 1 //1000～100hPaの等値線の設定 ga-> set cint 1 ga-> set cmax 10 ga-> d A1 ga-> set cmin 0.1 // 100～ 10hPaの等値線の設定 ga-> set cint 0.1 ga-> set cmax 1 ga-> d A2
上の例では、1000hPa～100hPaと100hPa～10hPaについて、等値線を描き分けています。
maskoutを用いた方法も示しておきます。こちらの方がスマートかもしれません。
ga-> set cmin 1 ga-> set cint 1 ga-> set cmax 10 ga-> d maskout(A, -(lev-1000)*(lev-100)) ga-> set cmin 0.1 ga-> set cint 0.1 ga-> set cmax 1 ga-> d maskout(A, -(lev-100)*(lev-10))
- (例)30S～30N以外を描画する
ga-> d maskout(A, (lat+30)*(lat-30))
ベクトル（vector）
ベクトルを表示する
set gxout vector d 変数1;変数2
変数1が横軸方向、変数2が縦軸方向の成分になります。
- (例)変数U,Vをベクトルで表示する
ga-> d U;V
ベクトルの色を変える
set ccolor 色番号
- (例)ベクトルの色を色番号2に設定する
ga-> set ccolor 2
ベクトルの太さを変える
set cthick 太さ
太さは1～20まで設定可能です。
- (例)ベクトルの太さを6に設定する
ga-> set cthick 6
ベクトルの長さスケールを変える
set arrscl 基準長 値
- 基準長：画面上での基準の長さ。画面右下の凡例に反映される。
- 値：ベクトルの長さが基準長であるときの、各成分の値の大きさ
- (例)画面上の長さが1の矢印(基準)が、大きさ200を表すようなベクトルを設定する
ga-> set arrscl 1 200
ベクトルの表示間隔を変える
d skip(変数1,Xステップ,Yステップ);変数2
- (例)ベクトル(U,V)をX=2,Y=3毎に表示する。
ga-> d skip(U,2,3);V
ベクトルの長さスケール(凡例)を表示する／しない
set arrlab on/off
好きな位置にスケールを表示できればいいのですが・・・。
バグ：ベクトルの向きがずれる
ベクトルを描く前に不正な処理（エラーに限らない）を行うと、ベクトルの向きがずれることがあります。これはreinitしても解消せず、GrADSの再起動が必要です。(気象研究所の方からの情報)
抜本的な対処法はGrADSのインストールを参照。
折線グラフ (line)
折線グラフを引く
set gxout line d 変数
次元が1次元の場合のデフォルトの図法です。
グラフの色を指定する
set ccolor 色番号
色番号はデフォルトの色一覧を参照してください。
グラフの太さを変える
set cthick 太さ
太さは1～20まで設定可能です。
- (例)グラフの太さを6に設定する
ga-> set cthick 6
マークの大きさを設定する
set digsiz サイズ
本家のマニュアルではサイズを0.10～0.15にするのがおすすめだそうです。
グラフのマーク(プロット点)の種類を指定する
set cmark マーク番号
グラフの線種を指定する
set cstyle 線種番号
y軸の範囲を指定する
set vrange 最小値 最大値
- (例)-3～3の値の範囲を指定
ga-> set vrange -3 3
緯度軸をcosでスケーリングする／しない
set coslat on/off
高緯度に行くほど単位緯度あたりの軸の長さが短くなります。全球平均値の緯度毎の寄与を見る際に有効かも。
未定義値を線で結ぶ / 元に戻す
set missconn on/off
折れ線の塗りつぶし（linefill）
あらかじめ次元を1次元にしておく必要があります。
2本の折れ線の間を塗りつぶす
set gxout linefill d 変数1;変数2
- (例)変数t_obsと変数t_modelの間を塗りつぶす
set gxout linefill d t_obs;t_model
- (例)変数uについて、値が10を越える領域を塗りつぶす
set gxout linefill d u;u-u+10
塗りつぶす色を指定する
set lfcols 色1 色2
「一番目の折れ線 > 二番目の折れ線」 のとき色1で、「一番目の折れ線 < 二番目の折れ線」 のとき色2で塗りつぶされます。
バー（bar）
あらかじめ次元を1次元にしておく必要があります。
バーを描く
set gxout bar d 変数
バーを塗りつぶす／塗りつぶさない
set baropts filled/outlined
バー同士の間隔を設定する
set bargap 間隔[%]
デフォルトの間隔は0です。間隔を100にするとバーではなく線で描かれます。
バーを描く基準値を設定する
set barbase 基準値 set barbase bottom/top
どの値を基準にしてバーを描くか設定します。値ではなくbottom(top)を指定すると、最小値よりも小さい値(最大値よりも大きい値)が基準になります。つまり、バーは常に上(下)に向かって伸びることになります。
- (例)t=300を基準にバーを描く
ga-> set barbase 300 ga-> d t
散布図（scatter）
散布図を描く
set gxout scatter d 変数1;変数2
- (例)横軸U、縦軸Vとして、t=1～100の範囲で散布図を描く。
ga-> set gxout scatter ga-> set t 1 100 ga-> d U;V
横軸、縦軸の範囲を設定する
set vrange 最小値 最大値 set vrange2 最小値 最大値
vrangeは横軸、vrange2は縦軸の値の範囲を設定します。
- (例)横軸U(0<=U<=100)、縦軸V(-10<=V<=10)として、散布図を描く。
ga-> set gxout scatter ga-> set vrange 0 100 ga-> set vrange2 -10 10 ga-> d U;V
マークの大きさを設定する
set digsiz サイズ
本家のマニュアルではサイズを0.10～0.15にするのがおすすめだそうです。
マークの種類を設定する
set cmark マーク番号
文字
大きさと幅を指定して文字列を描画する
set strsiz 幅 [高さ]
以下に、座標(3,3)に幅0.1、高さ0.15の大きさの文字で「hello」と描画する例を示します。
ga-> set strsiz 0.1 0.15 ga-> draw string 3 3 hello
色、座標の基準、太さ、回転角を指定して文字列を描画する
set string 色番号 [座標の基準 [太さ [回転角]]]
座標の基準とは、drawの際に指定する座標が、文字列のどの位置を基準にするかを表します。以下に指定のための文字列を示します。(GrADs本家より引用)
tl tc tr tl - top left
+-------------+--------------+ tc - center top
| | tr - right top
l + + c + r etc.
| |
+-------------+--------------+
bl bc brたとえば"c"の場合、drawで指定した座標を中心に文字列が配置されます。回転角(度)を指定した場合は、座標の基準を中心に時計回りに回転します。
- (例)色番号2(赤)を指定
ga-> set string 2
- (例)色番号2(赤)、座標の基準center、太さ0.1を指定
ga-> set string 2 c 0.1
- (例)色番号2(赤)、座標の基準center、太さ0.1、回転角90度を指定
ga-> set string 2 c 0.1 90
タイトルを書く
draw title This is Title.
すると「This is Title.」が結果画像の上に記入される。(事前に「d var」など、何か描いておく必要がある。)
別のやり方として、
set strsiz 0.2 0.25 set string 1 c 6 draw string 5.5 8.05 This is Title.
でも、似たように書くことが出来る。(事前に「d var」など、何か描いておく必要はない。)
文字中に「\」を使うことで、改行することもできる。
draw title This is\Title.
とすれば、以下のように出る。
This is Title.
X軸、Y軸のタイトル(ラベル)を書く
draw xlab X_LABEL draw xlab Y_LABEL
するとX軸には「X_LABEL」、Y軸には「Y_LABEL」と記入される。
上付き文字を表示する
`a文字列
なお、上付き文字を解除したい場合は`nを使います。
- (例)x^2+y^2=4 (^2は2乗)を上付き文字を使って表示
ga-> draw string 1.0 2.0 x`a2`n+y`a2`n=4
下付き文字を表示する
`b文字列
なお、下付き文字を解除したい場合は`nを使います。
- (例)HCl+NaOH=H2O+NaClを上付き文字を使って表示
ga-> draw string 1.0 2.0 HCl+NaOH=H`b2`nO+NaCl
様々な記号(e.g., ギリシャ文字)を表示する
`n文字列
nはフォント番号を表し、0から9まで指定可能です。デフォルトでは0～4のフォント番号が使えるはずです。`nの後の文字列には、表示したい記号に対応する文字を並べます。これは例えば
ga-> font 3
等で確認できます(左が入力する文字、右が実際に表示される文字)なお、デフォルトに戻る場合は`0を使います。
- (例)α≦β+γを表示する
ga-> draw string 1.0 2.0 `3a<b`0+`3c
図形
線を描く
draw line x1 y1 x2 y2
(x1,y1)から(x2,y2)へ直線を引きます。
- (例)座標(1,2)～(4,5)を結ぶ直線を引く
ga-> draw line 1 2 4 5
線の色、線種、太さを設定する
set line 色番号 [線種番号 [太さ]]
四角形を描く
draw rec xmin ymin xmax ymax
- (例)座標(1,2)～(4,5)に青(色番号4)で四角形を描く。
ga-> set line 4 ga-> draw rec 1 2 4 5
四角形(塗りつぶし)を描く
draw recf xmin ymin xmax ymax
- (例)座標(1,2)～(4,5)に青(色番号4)で内側を塗りつぶした四角形を描く。
ga-> set line 4 ga-> draw recf 1 2 4 5
多角形(塗りつぶし)を描く
draw polyf x1 y1 x2 y2 x3 y3 ... xn yn
- (例)青(色番号4)で内側を塗りつぶした5角形を描く
ga-> set line 4 ga-> draw polyf 1 1 2 1 3 2 2 2 1 3
マークを描く
draw mark マーク番号 x座標 y座標 マークの大きさ
マークの色を変える場合は
set line 色番号
天気記号を描く
draw wxsym 記号番号 x座標 y座標 大きさ [色 [太さ]]
記号番号は、
ga-> wxsym
と入力すると表示されます。
コントロールファイルの文法
データファイル名を指定する
バイナリデータのファイル名data-filenameを指定する場合、
DSET data-filename
ここでdata-filenameは、以下のような記述方法があります。
DSET ^test.dr //コントロールファイルから見た相対パス DSET /home/hoge/test.dr //絶対パス
うるう年を考慮しない
コントロールファイル内で、
OPTIONS 365_DAY_CALENDAR
と記述します。こうすると、1980年などの閏年でも一年が365日として扱われます。
複数ファイルを1つのコントロールファイルで読み込む
長い期間のデータを扱う場合、一つのファイルにまとめるとファイルサイズが大きくなることがあります。このようなとき、ファイルを分割すると扱いが楽になります。例えばデータを以下のように、月毎に分けることを考えます。
u199009.dr u199010.dr u199011.dr u199012.dr u199101.dr u199102.dr u199103.dr
このとき、コントロールファイルをひとつにまとめることができます。
DSET u%y4%m2.dr OPTIONS TEMPLATE // これが必要 .....(省略)..... TDEF 212 LINEAR 1SEP1990 1DY .....(省略).....
%y4には年が4桁で、%m2には月が2桁で置き換えられます。この他にも以下のような表記があります (よく分からないものには？が付いています…)。
| %y2 | 2桁の年 |
| %y4 | 4桁の年 |
| %m1 | 1桁又は2桁の月 |
| %m2 | 2桁の月(必要に応じて0が入る) |
| %mc | 3文字の月(JANなど？) |
| %d1 | 1桁又は2桁の日 |
| %d2 | 2桁の日(必要に応じて0が入る) |
| %h1 | 1桁又は2桁の時 |
| %h2 | 2桁の時 |
| %h3 | 3桁の時(例えば120又は012？) |
| %f2 | 2桁又は3桁の予報時間？ |
| %f3 | 3桁の予報時間？ |
| %n2 | 2桁の分(必要に応じて0が入る) |
| %n2 | 2桁の分(必要に応じて0が入る) |
| %ch | 任意の文字列（後述するCHSUBで指定） |
また、Version2.0aより、アンサンブルメンバーに関してもテンプレートで読み込めることが可能となりました。
| %e | 任意の文字列（tdefなどと同じように「edef 7 names 00 01 02 03 04 05 06」で指定可能） |
なお、
/mnt/sky/data/u1990_1991.dr /mnt/air/data/u1992_1993.dr
のように、任意の名前のファイルを連続して読み込むことも可能です。以下のように、タイムステップの始点・終点とファイル名を指定します。
DSET /mnt/%ch.dr OPTIONS TEMPLATE CHSUB 1 731 sky/data/u1990_1991 // はじめの731日分はこっち CHSUB 732 1461 air/data/u1992_1993 // 次の730日分はこっち .....(省略)..... TDEF 1461 LINEAR 1JAN1990 1DY .....(省略).....
エンディアンを指定して読み込む
OPTIONS little_endian //リトルエンディアンとして読み込む OPTIONS big_endian //ビッグエンディアンとして読み込む
他にもcray_32bit_ieeeというのがあるそうですが、誰か使っている人いますか？。
y方向またはz方向が逆順のバイナリファイルを使って描きたい場合の設定
コントロールファイルの「OPTIONS」に「yrev」または「zrev」を加える。
X次元の設定
二通りの設定方法があります。
XDEF 格子数 LINEAR 開始座標 増分 XDEF 格子数 LEVELS 座標1 座標2 ...
LINEARは等間隔格子、LEVELSは不等間隔格子の指定に適しています。基本的に単位は「度」ですが、kmなど他の単位で記述することも可能です。その場合、描画の前に
ga-> set mproj off
と指定するとよいでしょう。ただし地図を描きたい場合は、自前で地図データを用意して重ね描きする必要があります。
Y次元の設定
XDEFと基本的に同じです。
YDEF 格子数 LINEAR 開始座標 増分 YDEF 格子数 LEVELS 座標1 座標2 ...
X次元・Y次元・Z次元の設定についての注意点
異なるコントロールファイルのデータを比較したい場合(差・和を取りたいなど)、平面の場合ならX次元・Y次元の設定を、鉛直の場合ならX次元・Z次元、または、Y次元・Z次元の設定を、同じにしておく必要があるので注意!!
T次元の設定
TDEF タイムステップ数 LINEAR 開始時刻 時間増分
X次元などと違い、T次元は常にLINEAR（等間隔タイムステップ）です。開始時刻は
hh:mmZddmmmyyyy
もしくは簡単に、
ddmmmyyyy
等と指定します。ただし
| hh | 時(2桁)、省略時 00 |
| mm | 分(2桁)、省略時 00 |
| dd | 日(1 or 2桁)、省略時 1 |
| mmm | 月(3文字), 'JAN','OCT' など |
| yyyy | 年(2 or 4桁)、2桁の場合、1950-2049年の下2桁。 |
時間増分は、'1yr', '2mn' などのように、「2桁以内の整数 + 単位」と表現します。単位は以下が使えます。
| mn | 分 (minute) |
| hr | 時 (hour) |
| dy | 日 (day) |
| mo | 月 (month) |
| yr | 年 (year) |
4バイト浮動小数以外のバイナリデータを読み込む
変数を列記する部分
変数名 層の数 unit 説明
のunitフィールドに指定します。指定方法は以下の通り（空白を入れないように注意！）。
| unit | データ |
|---|---|
| -1,40,1 | 1バイト符号なし整数 |
| -1,40,2 | 2バイト符号なし整数 |
| -1,40,2,-1 | 2バイト符号あり整数 |
| -1,40,4 | 4バイト整数 |
- (例) 1バイト整数のデータ(変数名sst)
VAR 1 sst 0 -1,40,1 SST ENDVARS
コントロールファイルの例(書式なしファイルの場合)
DSET ^%y4%m2.dr UNDEF 9.999E+20 TITLE JRA-25 1.25 OPTIONS YREV LITTLE_ENDIAN TEMPLATE XDEF 288 LINEAR 0.000000 1.2500000 YDEF 145 LINEAR -90.000000 1.25 ZDEF 1 LEVELS 1013 TDEF 120000 LINEAR 00Z01JAN1979 6HR VARS 1 prmsl 0 33,100,0 ** pressure reduced to mean sea level [Pa] ENDVARS
NetCDF用のコントロールファイルを作成する
sdfopenでNetCDFファイルを開けない場合、コントロールファイルを作成するとうまくいく場合があります。これには、
- 完全なコントロールファイルを作成する方法(open使用)
- NetCDFヘッダの一部を上書きする方法(xdfopen使用)
の二通りがあります。
前者は基本的に通常のコントロールファイルと同様です。ここに例があります。
後者は xdfopen というコマンドを使用します。文法は通常のコントロールファイルと良く似ていますが、最低限
DSET NetCDFファイル名
があればOKです。
次元を設定するには、
XDEF NetCDF中の次元名 XDEF NetCDF中の次元名 サイズ XDEF NetCDF中の次元名 サイズ [linear/levels] 値・・・
などとします。「NetCDF中の次元名」とは、ncdump -h などで出てくる次元の名前です。３番目の使用方法のサイズ以降は、ふつうのコントロールファイルの設定方法と同様です。たとえば、
XDEF lon YDEF lat TDEF time 1200 LINEAR 00:00Z15JAN2000 1mo
とすると、時間次元のみ上書きされます。
以下、完全なファイル(a.ctl)の例。
DSET ^a.nc XDEF lon YDEF lat TDEF time 1200 LINEAR 00:00Z15JAN2000 1mo
開くときは、
ga-> xdfopen a.nc
ただし、GrADS Version 1.9以下では gradsnc を使用してください。
コントロールファイルの例(NetCDFファイル、openを用いる場合)
dset ^%y4%m2%d2.nc options template yrev undef 9.96921e+36 unpack scale_factor add_offset dtype netcdf xdef 240 linear 120.0000 0.125 ydef 252 linear 22.4000 0.100 zdef 1 levels 1013 tdef 40000 linear 00z02jul2002 1hr vars 9 psea=>psea 0 t,y,x sea level pressure u=>u 0 t,y,x eastward component of wind v=>v 0 t,y,x northward component of wind temp=>temp 0 t,y,x temperature rh=>rh 0 t,y,x relative humidity r1h=>r1h 0 t,y,x rainfall in 1 hour ncld_upper=>ncld_upper 0 t,y,x upper-level cloudiness ncld_mid=>ncld_mid 0 t,y,x mid-level cloudiness ncld_low=>ncld_low 0 t,y,x low-level cloudiness endvars
緯度経度座標以外のデータの読み込み（PDEF）
GrADSでは緯度経度座標だけでなく任意の座標系のデータを読み込むことができます。特にランベルト座標やポーラーステレオ座標のようによく使われる座標系については、コントロールファイルに記述するだけで簡単に読み込むことができます。以下ではPDEFの使い方を簡単に説明してきます。詳細はGrADS本家 をご覧ください。
ランベルト座標
NHMから直接出力されたデータやMANALなどを読み込む場合に便利です。
PDEF XNUM YNUM LCCR LAT LON X Y SLAT1 SLAT2 SLON DX DY PDEF XNUM YNUM LCC LAT LON X Y SLAT1 SLAT2 SLON DX DY
| XNUM, YNUM | X, Y方向の格子数。 |
| LAT, LON | データ上のある1点の緯度,経度。 |
| X, Y | 上記で指定したある1点の緯度,経度の格子番目。GrADSは「options yrev」以外は南西が開始点なので、南および西から数えた点。 |
| SLAT1, SLAT2, SLON | 投影基準点。 |
| DX, DY | 格子間隔 [m] |
前者（LCCR）は風が図法上の成分で書いてある場合、後者（LCC）は風が緯度経度座標の成分にすでに変換されている場合に用います。
以下、MANALを読み込む場合の例。
PDEF 361 289 LCCR 30 140 245 85 30 60 140 10000 10000
ただし、Version 1.9 b の場合、PDEFを定義した場合に lat、lon 変数がうまく描かれません。Version 2.0 a では修正されています。
一般的な座標
GrADSで定義されていない座標系を用いる場合、以下のように BILIN キーワードを用いてコントロールファイルに記述します。ただし、後で示すように座標変換を行うためのマッピングファイルを作成する必要があります。
PDEF isize jsize BILIN format byteorder fname
| isize | 変換前座標のi格子数 |
| jsize | 変換前座標のj格子数 |
| format | マッピングファイルが direct access の場合 STREAM、sequential access の場合 SEQENTIAL を指定。 |
| byteorder | マッピングファイルのエンディアン。BINARY-BIG または BINARY-LITTLE。 |
| fname | マッピングファイル名。 |
マッピングファイルを作るのに必要な最低限の情報は、
変換後の格子（緯度経度座標）における、変換前格子のグリッド番号
です。簡単のため、以下のような座標で考えてみましょう。
| | | (1,4)-------(2,4)-------(3,4)--- | | | | <3,8> |<4,8> | | | | (1,3)-------(2,3)-------(3,3)--- | | | | <3,7>| <4,7> | | | (1,2)-------(2,2)-------(3,2)--- | | | | |<3,6> | | | | <4,6> (1,1)-------(2,1)-------(3,1)---
ここで(*,*)は変換後の緯度経度座標のグリッド番号、<*,*>は変換前の一般座標のグリッド番号を表しているものとします。このとき、たとえば(2,3)の格子点は、変換前の格子点<3,8>, <4,8>, <3,7>, <4,7> に囲まれています。大体<3.8,7.5>といったところでしょうか。この「(2,3)は変換前格子の<3.8,7.5>に相当する」という情報が必要となります。GrADSはこの情報をもとに線形内挿を行うのです。
ここでは触れませんが、「風速がどの座標系に基づいているか」によっては単純な線形内挿では不十分な場合があります。風速(u,v)が緯度経度座標を向いていない場合、どれくらい回転することで緯度経度座標に直せるかという情報も必要になります。
マッピングファイルには、XDEF*YDEF個の4バイト浮動小数点データが3組格納されています。
- 各緯度経度座標におけるi座標値
- 各緯度経度座標におけるj座標値
- 各緯度経度座標における風速の回転の大きさ i座標値として -999 を入れると未定義値扱いとなります。また、回転の大きさを-999にすると回転なしになりますので、とりあえず-999にするとよいでしょう。
変数
予約済みの変数
lat lon lev
これらはコントロールファイルで開いた変数と同様に使用できます。
- (例)緯度線を表示する
d lat
- (応用例)場所によって等値線の間隔を変える
関数
注意事項
「http://meteora.ucsd.edu/~kyoshimura/?IT%20memo%2FGrADS%20memo#ea5481a4」を参考に、関数使用時の注意点を以下にあげる。
以下を踏まえて、注意して扱っていくこと。
- 欠損値の扱い(ave,sum)
ave(var,t=1,t=10)*10 sum(var,t=1,t=10)
以上二つは同じではない！（欠損値の取り扱いが原因）
sum → 欠損値をゼロとして（欠損値を飛ばして）計算する。 ave → sumの結果を欠損値を除いた個数で割る。
欠損値を含んでいる場合などは注意して扱わなければいけない。
- 平均計算の順序(時間平均,空間平均)
平均をするときは、時間平均→空間平均の順にする。
× 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)
- aveとaaveの違い
aveとaaveはundefのデータを含んでいなければ
ave(ave(var,x=1,x=10),y=1,y=20) aave(var,x=1,x=10,y=1,y=20)
は同じになる。しかし、undefを含んでいるような以下のようなデータ(Uはundef)の場合を考えてみる。
3 4 6 7 10 10 10 10 15 U U U
aveで求めると、Undefを除いて計算するので、各行の平均は、1行目が5、2行目が10、3行目が15となり、全体の平均は(5+10+15)/3=10となる。
一方、aaveで求めると、(3+4+6+7+10+10+10+10+15)/9=75/9=8.33となり、両者は異なる。
Undefを含んだデータの平均の場合はaaveを使用する必要がある。
- aveとmean(aaveとamean)
aveは緯度による重み付けをした平均に対して、meanは重み付けをしない平均になる。欠損値の扱いは、どちらも同じ。
aave,ameanも同じ関係。 - asumについて
asumは計算がうまく行っていない気がするので(バグか？)、
○ sum(sum(var,x=1,x=10),y=1,y=20) × asum(var,x=1,x=10,y=1,y=20)
平均を計算する
ave(変数, 開始, 終了 [,増分] [,-b]) mean(変数, 開始, 終了 [,増分] [,-b])
平均(時間 or 空間)を計算します。平均の開始と終了の指定方法は例を御覧下さい。増分は時間平均の際の時間増分を指定することができます。また、-bは境界を正確に扱うためのフラグです。
格子がlinearでない場合、格子間隔に応じて重みをつけて平均されます。また、緯度平均を取る場合、aveを用いると緯度に応じた重みがかかりますが、meanでは重みがかかりません。
以下の例では、空間3次元と時間の次元をもつ変数Uを取り上げて説明します。全ての例は以下の次元の状態からスタートするものとします。
ga-> q dims //とりあえず現在の次元を表示 Default file number is: 1 X is varying Lon = 0 to 360 X = 1 to 129 Y is varying Lat = -87.864 to 87.864 Y = 1 to 64 Z is fixed Lev = 1000 Z = 1 T is fixed Time = 00Z01AUG1998 T = 1
- (例)t=1～t=5で時間平均を行う(高度1000hPa面のまま)。
ga-> a=ave(u, t=1, t=5) Averaging. dim = 3, start = 1, end = 5 Define memory allocation size = 33024 bytes ga-> d a
- (例)01aug1998～10aug1998で時間平均を行う(高度1000hPa面のまま)。
ga-> b=ave(u,time=01aug1998,time=10aug1998) Averaging. dim = 3, start = 1, end = 37 Define memory allocation size = 33024 bytes ga-> d b
- (例)東西平均(帯状平均)を行う。1000hPaにおける東西平均値を折れ線グラフで表示する。
ga-> set lon 0 ga-> a=ave(u,lon=0,lon=360,-b) Averaging. dim = 0, start = 1, end = 129 Define memory allocation size = 256 bytes ga-> d a
最初のset lonは、便宜上指定したものです。-bをつけると、周期的な平均で値がより正確になります。
- (例)東西平均(帯状平均)を行う。1000hPa～1hPaにおける東西平均値を等値線で表示する。
ga-> set lon 0 ga-> set lev 1000 1 LEV set to 1000 1 ga-> b=ave(u,lon=0,lon=360,-b) Averaging. dim = 0, start = 1, end = 129 Averaging. dim = 0, start = 1, end = 129 Averaging. dim = 0, start = 1, end = 129 Averaging. dim = 0, start = 1, end = 129 Averaging. dim = 0, start = 1, end = 129 Averaging. dim = 0, start = 1, end = 129 Averaging. dim = 0, start = 1, end = 129 Averaging. dim = 0, start = 1, end = 129 Averaging. dim = 0, start = 1, end = 129 Averaging. dim = 0, start = 1, end = 129 Averaging. dim = 0, start = 1, end = 129 Averaging. dim = 0, start = 1, end = 129 Averaging. dim = 0, start = 1, end = 129 Averaging. dim = 0, start = 1, end = 129 Averaging. dim = 0, start = 1, end = 129 Averaging. dim = 0, start = 1, end = 129 Averaging. dim = 0, start = 1, end = 129 Define memory allocation size = 4096 bytes ga-> d b
- (例)気候値を計算する。ここでは10年分の月平均データを用いた例。
ga-> set t 1 12 ga-> a=ave(u,t+0,t=120,12mn)
最初に指定したtの範囲で、aveのt+0にtが代入されてそれぞれ平均されます。
上記で「12mn」の部分には「set t」で指定した期間（上記では1年間）を記入します。
（単位は、コントロールファイルの「tdef」の単位と同じのにする必要がありそうです。）
- (例)気候値からの偏差を計算する。ここでは10年分の月平均データを用いた例。
ga-> set t 1 12 ga-> clim=ave(u,t+0,t=120,12mn) ga-> modify clim seasonal ga-> set t 1 120 ga-> anomaly=u-a
modifyは年に依存しない変数を指定するコマンドです。つまり、何年を指定しても、同じ値が出ます。
ga-> set time jan1960 ga-> d clim ga-> set time jan1970 ga-> d clim
以上2つの値は、同値。
他に、日に依存しない変数を指定することも可能です。
modify 変数名 seasonal modify 変数名 diurnal
- (例)移動平均を計算する。時間ステップが100のuデータに対して、両側3ステップを平均する場合、
ga-> set t 4 97 ga-> a=ave(u,t-3,t+3)
- (例)移動平均を計算する。2000年1月のslpデータに対して、両側1日を平均する場合。
ga-> set time 02jan2000 30jan2000 ga-> a=ave(slp,time-1dy,time+1dy)
領域平均を計算する
最初の注意事項に書いた通り、1次元方向でなく、2次元の領域で平均する場合にはaaveを用いる(欠損値の扱いの都合上)。
d aave(var,x=1,x=10,y=1,y=20)
ただし、このままでは時間方向に動かすことができない(「set t 1 10」の後に上記コマンドを実行しても、エラーメッセージが出るだけ)。
「tloop」関数を使用することで、時間方向に動かすことができ、時系列の変化をグラフにして出すことが可能になる。
d tloop(aave(var,x=1,x=10,y=1,y=20))
最大値／最小値を計算する
max(変数, 開始, 終了 [,時間増分]) min(変数, 開始, 終了 [,時間増分])
- (例)40Nの東西風uの最大値を計算する。
ga-> set lat 40 ga-> set lon 0 ga-> d max(u,lon=0,lon=360)
最大値／最小値が出現する位置を求める
maxloc(変数, 開始, 終了 [,時間増分]) minloc(変数, 開始, 終了 [,時間増分])
使い方は基本的にmax()、min()と同様です。
- (例)40Nの東西風uが最大となる位置を求める。
ga-> set lat 40 ga-> set lon 0 ga-> d maxloc(u,lon=0,lon=360)
鉛直積分を計算する
vint(開始気圧, 物理量, 終了気圧)
開始気圧(hPa)から終了気圧(hPa)までの物理量の質量重みつき鉛直積算を行います。開始気圧と物理量は次元を持った変数、終了気圧は定数で指定します。
- (例)エネルギーEをpsから10hPaまで鉛直積算して表示する
d vint(ps, E, 10)
- (例)エネルギーEを1000hPaから10hPaまで鉛直積算して表示する
d vint(E(z=1)-E(z=1)+1000, E, 10)
開始気圧は次元をもつ必要があるため、E(z=1)-E(z=1)という領域全体で0の変数を用意し、それに必要な気圧を加えます。
中心差分を計算する
cdiff(変数, 次元)
次元にはX, Y, Z, T どれかを指定します。指定した次元はvaryingである必要があります。
- (例)Uの緯度方向の微分（du/dlat）を行う（1000hPa面）
ga-> set lon 0 360 ga-> set lat -90 90 ga-> set lev 1000 ga-> du=cdiff(u,y) ga-> dy=cdiff(lat,y) ga-> dudy=du/dy
dudyに微分が格納されます。計算は中心差分で行っており、関数cdiffはすべてのjについて u(j+1)-u(j-1) という計算を行っています。
- バグ？：次元にZを指定すると、「d cdiff(u,z)」ならうまく表示されますが、「dz=cdiff(u,z)」とすると「Data Request Warning～」もしくは「Specifued Dimension non varying」となって失敗してしまいます。以下に代替案を示します（情報提供：気象庁の上杉さん）
- (例)Uの鉛直微分（du/dlev）を行う（lon=0面）。
- cdiffと同様、中心差分を行っています（上流差分にも応用できます）。
- u(z+1), u(z-1)を使うので、levの指定には上下1層以上の余裕が必要です。
ga-> set lon 0 ga-> set lat -90 90 ga-> set lev 925 150 ga-> du=u(z+1)-u(z-1) ga-> dz=lev(z+1)-lev(z-1) ga-> dudz=du/dz
水平渦度(curlの鉛直成分)を計算する
hcurl(東西風, 南北風)
東西風と南北風の変数から、水平渦度を計算します。
- (例)水平渦度を計算して表示
ga-> a=hcurl(u,v) ga-> d a
水平発散(div)を計算する
hdivg(東西風, 南北風)
東西風と南北風の変数から、水平発散を計算します。
- (例)水平発散を計算して表示
ga-> a=hdivg(u,v) ga-> d a
データにマスクをかける
maskout(expr, mask)
exprは描画したい変数、maskはマスクの判定を行うための変数です。ある地点におけるmaskの値が0以上の場合、exprをそのまま返します。maskの値が0未満の場合は、未定義値として扱われます。
- (例)陸上ならmask=1、海上ならmask=-1という値が入ったマスクがある場合、
ga-> d maskout(u,mask) //陸上のみuを描画 ga-> d maskout(u,-mask) //海上のみuを描画 ga-> d aave(maskout(u,mask),lon=0,lon=360,lat=0,lat=90) //陸上のみ平均
未定義値／未定義値以外／全値を強制的に定数にする
const(expr, value, -u) // exprの未定義値 → value const(expr, value) // exprの未定義値以外 → value const(expr, value, -a) // exprの全ての値 → value
- (例) 山岳（未定義値）を塗りつぶした気温ta(>0)の等値線図を描く
set gxout shaded set clevs 0 set ccols 1 0 d const(ta,-100,-u) set gxout contour d ta
時間についての相関係数を計算する
tcorr(変数1, 変数2, t=開始時間ステップ, t=終了時間ステップ) or tcorr(変数1, 変数2, time=開始時刻, time=終了時刻)
時間方向の相関係数を計算します。変数1は時間の次元しか持てないことに注意してください(変数2は時間以外を持つことも可)。
- (例)U.1とU.2の相関を求める
ga-> d tcorr(U.1, U.2, t=1, t=10)
線形回帰係数を計算する
tregr(変数1, 変数2, t=開始時間ステップ, t=終了時間ステップ) or tregr(変数1, 変数2, time=開始時刻, time=終了時刻)
変数1、2はともに時間に依存する必要があります。変数2は時間に加えて空間(X,Y)に依存することもできます。横軸に変数1、縦軸に変数2をとった散布図を書いたとき、最小自乗法でフィッティングした直線の傾きを求めていることになります。
- (例)変数aoに回帰したzを求める
tregr(ao, z, t=1, t=30)
べき乗を計算する
pow(変数1, 変数2)
変数1の「変数2」乗を計算します。
- (例)Uから運動エネルギーを求める
ga-> d 0.5 * pow(u, 2)
スムージングを行う
smth9(変数)
自身＋周囲8格子の重みつき平均を行います。重みは、自身が1、上下左右が0.5、斜め上・下が0.3です。但し、1次元データの場合は自身＋周囲2格子のデータを使用します。
水平・時間内挿を行う(Ver2.0以降)
lterp(変数1, 変数2)
変数1を変数2の格子に内挿する計算をします。格子数・時間間隔の異なるデータ間の計算を、この関数を用いることで出来るようになります。(鉛直方向の内挿は出来ない)
ただし、以下の注意が必要。
関数を使用する場合は、変数1・変数2が未定義値を含んでてはダメ。未定義値を含むデータを内挿すると計算が落ちます(セグメンテーションエラー)。
2つの領域が異なる場合、小さい方の領域に範囲をセット(set x, set y)してから計算を行う必要がある。
鉛直座標が圧力・高度の場合、地表面付近は未定義値が入っているので、注意。
- (例)JRA125の気温データをNCEPの格子に内挿し、JRA125-NCEPの差を見てみる
ga-> open NCEP/t.ctl ga-> open JRA125/t.ctl ga-> set time 01jan2000 ga-> d lterp(t.2,t.1)-t.1
ユーザ定義関数
GrADSではave()やpow()のような組み込み関数だけでなく、独自の関数を追加することができます。これをユーザ定義関数（User Defined Function, 以下UDF）と呼びます。
UDFは自分の好きなコンピュータ言語(たとえばFortran、C言語)を用いて作成することができます。以下では実際に関数を作成しながらUDFの解説をします。
注意！：2009.01.18現在、2.0系ではユーザ定義関数が使用できません。
ga-> q udf Warning: User Defined Functions have been disabled in this version of GrADS
となります。GrADSソースコードのgafunc.cに
JMA remainder of subroutine commented out pending implementation of DLLs
とあるので、そのうち復活するのではと期待しています・・・。
UDFT (ユーザ定義関数テーブル)
UDFTとは関数の名前やプログラムの場所、引数の渡し方などを定義したテキストファイルです。任意の場所に置くことができます (ここでは/usr/local/grads/にudftというファイルを作成します)。以下にUDFTの例を示します。
test 1 3 expr value value sequential /usr/local/grads/test/test /home/kodama/grads/test/test.out /home/kodama/grads/test/test.in
各行の簡単な説明は以下の通りです。
(関数名) (引数の数の最小値) (引数の数の最大値) (引数の種類1) (引数の種類2) ... (関数データ交換ファイルの読み書き方法) (実行ファイル名) (関数データ交換ファイル名) (関数結果ファイル名、実行ファイル->GrADS)
| (関数名) | 関数の名前。1文字以上8文字以内。1文字目は必ずアルファベット。大文字と小文字の区別はされないので注意。 |
| (引数の数の最小値) | 関数に与える引数の数の最小値。 |
| (引数の数の最大値) | 関数に与える引数の数の最大値。ちなみに引数の数は8個以内。 |
| (引数の種類 n) | 引数の種類。引数の数（の最大値）だけ列挙する。expr:グリッドデータ等、value:数値、char:文字列。 |
| (関数データ交換ファイルの読み書き方法) | sequential or direct。 |
| (実行ファイル名) | 実行ファイルの名前（e.g. a.out） |
| (関数データ交換ファイル名) | GrADSからプログラムにデータを渡す際の中間ファイル名。 |
| (関数結果ファイル名 | プログラムの結果をGrADSに渡す際の中間ファイル名。 |
次に、環境変数GAUDFTにUDFTのファイル名を設定します。 bashなら
[mydir@hosts] export GAUDFT="/usr/local/grads/udft"
cshなら
[mydir@hosts] setenv GAUDFT /usr/local/grads/udft
ここまでうまくいったかどうかを確認してみましょう。GrADSを立ち上げて、
ga-> q udft
によって、現在登録されているUDFを確認することができます。 複数のUDFをUDFTに記述する場合のUDFTは、
test 1 1 expr sequential /usr/local/grads/test/test /home/kodama/grads/test/test.out /home/kodama/grads/test/test.in test2 1 1 expr sequential /usr/local/grads/test2/test2 /home/kodama/grads/test2/test2.out /home/kodama/grads/test2/test2.in
のように前のテーブルの後ろにつなげて記述します。決して空行を空けたりしてはいけません。
関数データ交換ファイル
GrADSからプログラムへデータを渡すには、 UDFT4行目に指定した関数データ交換ファイル(function data transfer file)を介して行います。ファイルはヘッダレコードと引数レコードに分かれています。
ヘッダレコードは浮動小数点(4bytes)20個で構成されています。各浮動小数点の意味は以下の通り。
- 関数が呼ばれたときの引数の数
- 常に0。関数データ交換ファイルのフォーマットを表す。(0以外の場合はエラーらしい)
- これ以降は現在未定義(～20)
引数部には、UDFで指定した引数の順番にデータが入っています。引数の種類（expr, value, char）によって書式が違います。
exprを指定した場合、
- 第1レコード(4bytes浮動小数点×20個)
- UNDEFの値
- i次元の種類を表す値(idim)
- -1:なし
- 0:X次元(経度)
- 1:Y次元(緯度)
- 2:Z次元(高度)
- 3:T次元(時間)
- j次元の種類を表す値(jdim: 値の意味はitypeと同じ)
- i次元の要素の数(isiz)
- j次元の要素の数(jsiz)
- i次元が線形か非線形か(ilinear: 0なら非線形)
- j次元が線形か非線形か(jlinear: 0なら非線形)
- i次元の開始値(istrt)、ただしi次元が線形の場合
- i次元の増分(iincr)、ただしi次元が線形の場合
- j次元の開始値(jstrt)、ただしj次元が線形の場合
- j次元の増分(jincr)、ただしj次元が線形の場合
- 時間次元の開始時刻(年)、ただしどちらかの次元が時間の場合
- 時間次元の開始時刻(月)、ただしどちらかの次元が時間の場合
- 時間次元の開始時刻(日)、ただしどちらかの次元が時間の場合
- 時間次元の開始時刻(時)、ただしどちらかの次元が時間の場合
- 時間次元の開始時刻(分)、ただしどちらかの次元が時間の場合
- 時間次元の増分(分)、ただしどちらかの次元が時間の場合
- 時間次元の増分(月)、ただしどちらかの次元が時間の場合
- これ以降は現在未定義(～20)
- 第2レコード
- グリッドデータ(isiz×jsiz個の4bytes浮動小数点)
- 第3レコード
- i次元の座標値(isiz個の4bytes浮動小数点)
- i次元が存在しない(idim=-1)場合、このレコードは存在しない。
- i次元が存在する(idimが-1ではない)場合、線形格子(ilinear=1)であってもこのレコードは存在する。
- i次元の座標値(isiz個の4bytes浮動小数点)
- 第4レコード
- j次元の座標値(jsiz個の4bytes浮動小数点)
- j次元が存在しない(jdim=-1)場合、このレコードは存在しない。
- j次元が存在する(jdimが-1ではない)場合、線形格子(jlinear=1)であってもこのレコードは存在する。
- j次元の座標値(jsiz個の4bytes浮動小数点)
valueを指定した場合、4bytes浮動小数点×1個。これがそのまま引数の値です。
charを指定した場合、80bytesの文字列。80bytesに満たない部分には空白で埋められ、80bytesを超える場合は捨てられます。
関数結果ファイル
プログラムで処理した結果をGrADSへ戻すには、UDFT5行目に指定した関数結果ファイル(function result file)を介して行います。ファイルはヘッダレコードとグリッドレコードに分かれています。
ヘッダレコードは浮動小数点(4bytes)20個で構成されています。各浮動小数点の意味は以下の通り。
- プログラムの戻り値。0で正常終了、0以外で異常終了。異常終了の場合はこれ以上のデータ読み込みは行われない。
- 常に0。関数結果ファイルのフォーマットを表す。
- これ以降は現在未定義(～20)
グリッドレコードは基本的に、関数データ交換ファイルの引数部と同じ形式で記述します。ただし第3レコードと第4レコードは扱いが少し異なるので注意が必要です。
- 第3レコード
- i次元の座標値(isiz個の4bytes浮動小数点)
- i次元が存在しない(idim=-1)場合、又は線形格子(ilinear=1)の場合、このレコードは存在しない。
- i次元が存在し、かつ非線形格子である場合に限り、このレコードは存在する。
- i次元の座標値(isiz個の4bytes浮動小数点)
- 第4レコード
- j次元の座標値(jsiz個の4bytes浮動小数点)
- j次元が存在しない(jdim=-1)場合、又は線形格子(jlinear=1)の場合、このレコードは存在しない。
- j次元が存在し、かつ非線形格子である場合に限り、このレコードは存在する。
- j次元の座標値(jsiz個の4bytes浮動小数点)
簡単な例
2次元の気温場T[K]を与えると、単位を[℃]に変換する関数kelvin()を作成します。関数とデータ交換用のファイルは /home/hoge/gudf/kelvin 以下に置くことにします。予め環境変数GAUDFTを設定するのも忘れずに。
- udf
kelvin 1 1 expr direct /home/hoge/gudf/kelvin/a.out /home/hoge/gudf/kelvin/kelvin.out /home/hoge/gudf/kelvin/kelvin.in
- kelvin.f90
program kelvin
implicit none
real(4) :: header(20), argh(20)
integer :: def(2)
real(4),allocatable :: gdata(:,:)
integer :: i, j, pos
write(*,*) 'This is kelvin.f90'
open(8, file='/home/hoge/gudf/kelvin/kelvin.out', &
& status='old', form='unformatted', access='direct', recl=4)
open(10,file='/home/hoge/gudf/kelvin/kelvin.in', &
& status='new', form='unformatted', access='direct', recl=4)
!***** Input from GrADS output file (kelvin.out) *****!
! Header
pos = 1
do i=1, 20
read(8,rec=pos) header(i)
pos = pos + 1
end do
!*** 1: expr data ***!
! 1st Record
do i=1, 20
read(8,rec=pos) argh(i)
pos = pos + 1
end do
def(1) = int( argh(4) )
def(2) = int( argh(5) )
allocate( gdata( def(1), def(2) ) )
write(*,*) 'Dimension : ', def(1), def(2)
! 2nd Record
do j=1, def(2)
do i=1, def(1)
read(8,rec=pos) gdata(i,j)
pos = pos + 1
end do
end do
!***** Analysis *****!
gdata(:,:)= gdata(:,:) - 273.15
!***** Output to GrADS input file (test.in) *****!
! 1st Record
pos = 1
do i=1, 20
write(10,rec=pos) 0.0
pos = pos + 1
end do
! 1st Record of arguement header
do i=1, 20
write(10,rec=pos) argh(i)
pos = pos + 1
end do
! 2nd Record of arguement header
do j=1, def(2)
do i=1, def(1)
write(10,rec=pos) gdata(i,j)
pos = pos + 1
end do
end do
deallocate( gdata )
stop
end program kelvin- 注意
- バイナリデータの読み書きにはdirectを使っています。が、sequentialを使った方が楽＆速いかも。
kelvin.f90コンパイルしてa.outを作成すれば、準備は完了です。GrADSを立ち上げて、
d kelvin(T)
とすれば実行できます。
その他
色を定義する
set rgb 色番号 red green blue
RGB値を色番号に結びつけます。色番号は0～15が定義済み(デフォルトの色一覧)なので、16～99が使えます。red, green, blueは0～255の整数です。
シェルのコマンドを実行する
!コマンド名
!を先頭に付けると、シェルのコマンドと解釈されます。 gradsを実行しながらlsを実行したい場合は、
ga-> !ls index.html index.html~
とします。
ga-> !emacs test.ctl &
みたいなこともできます(GrADSを閉じずにコントロールファイルを編集できる)。
カラーバーを表示する
cbar
「set gxout shaded」などの設定がされている場合、
ga-> cbar
とすると、適切な位置にカラーバーが表示されます。
ga-> cbarc ga-> cbarn
などのバリエーションもあります。これらのカラーバーは、/usr/local/lib/grads/ といったフォルダに、スクリプトとして置いてあります。 任意の場所にカラーバーを描画したい場合はスクリプト集を参考にして下さい。
cbar*がシステムに存在しない場合は、
ftp://grads.iges.org/grads/scripts/
からcbar*.gsをダウンロードして任意の場所(/usr/local/lib/grads/scripts/ など)に置き、 環境変数GASCRPにダウンロードしてきたファイルを置いたディレクトリを指定すると使えるようになります。
バイナリデータの格納順
たとえば以下のようなコントロールファイル(test.ctl)があったとします。
DSET test.dr UNDEF 9.99E+33 TITLE UVT XDEF 3 LEVELS 0 120 240 YDEF 4 LEVELS -90 -30 30 90 ZDEF 5 LEVELS 1000 800 600 400 200 TDEF 10 LINEAR 1JAN2000 1DY VARS 3 U 5 99 horizontal velocity V 5 99 horizontal velocity T 5 99 temperature ENDVARS
このとき、対応するバイナリファイルtest.drの中身は(図はbinary.png)
U(x=1, y=1, z=1, t=1) U(x=2, y=1, z=1, t=1) U(x=3, y=1, z=1, t=1) U(x=1, y=2, z=1, t=1) U(x=2, y=2, z=1, t=1) U(x=3, y=2, z=1, t=1) U(x=1, y=3, z=1, t=1) …………(中略)………… U(x=3, y=4, z=1, t=1) U(x=1, y=1, z=2, t=1) U(x=2, y=1, z=2, t=1) U(x=3, y=1, z=2, t=1) U(x=1, y=2, z=2, t=1) …………(中略)………… U(x=3, y=4, z=5, t=1) V(x=1, y=1, z=1, t=1) V(x=2, y=1, z=1, t=1) …………(中略)………… T(x=2, y=4, z=2, t=1) T(x=3, y=4, z=2, t=1) ………(以上をt=2,3,...についても繰り返す)……… U(x=1, y=1, z=1, t=2) U(x=2, y=1, z=1, t=2) U(x=3, y=1, z=1, t=2) U(x=1, y=2, z=1, t=2) ……………………………
Fortran 90 でデータを吐き出す例(test.f90)を以下に示しました。(注意：open文のreclの指定方法は、コンパイラやオプションによって異なります。ifortで特にオプションがない場合は4*を付けない等)
program test
implicit none
integer,parameter :: xdef=3, ydef=4, zdef=5
real(4) :: U(xdef,ydef,zdef), V(xdef,ydef,zdef), T(xdef,ydef,zdef)
integer :: i, r
open(10,file='test.dr',status='unknown',form='unformatted', &
& access='direct',recl=4*xdef*ydef*zdef)
r = 1 !書きこみ位置
do i=1,10
! ここで時刻t=iのU,V,Tを格納
…………(中略)…………
! 書きこみ
write(10,rec=r) U
r = r + 1
write(10,rec=r) V
r = r + 1
write(10,rec=r) T
r = r + 1
end do
close(10)
end program test計算したデータなどをGradsファイルで出力する
「set gxout fwrite」を使う。以下はnetCDFからGrADSファイルを作成する場合の例。しっかりと格子情報を指定してからにしてください。(うまくいかないときがあるらしい・・・)
'sdfopen netcdf.nc' 'set gxout fwrite' 'set undef dfile' ;* <- GrADS 2.0.a6以降の場合あった方がいいかも 'set fwrite test.grd' t=1 while(t<=tmax) 'set t 't z=1 while(z<=zmax) 'set z 'z 'set x 1 96' 'set y 1 96' 'd var' z=z+1 endwhile t=t+1 endwhile 'disable fwrite'
出力時にはオプションがつけられる。「set fwrite (-le,-be) (-sq,-st) (-ap,-cl) test.grd」。
「-le」-> littel endian 出力(default)、「-be」-> big endian 出力。
「-sq」-> sequential 出力、「-st」-> stream(direct) 出力(default)。
「-ap」-> (ファイルが存在した場合)指定先のファイルに書き足す、「-cl」-> (ファイルが存在した場合)指定先のファイルを書き換える(default)。
'set fwrite'を複数回実行する場合、'disable fwrite'でファイルを閉じる必要があります。
GrADS 2.0.a6以降の場合、未定義値は自動的に-9.99e+8になってしまいます。未定義値を元のファイルに合わせるには、'set undef dfile' とするとよいでしょう。
出力する未定義値を変更する (Ver.2.0.a6以降)
set undef 未定義値 set undef file ファイル番号 // 指定したファイル番号の未定義値を使用 set undef dfile // 現在のファイル番号の未定義値を使用
ここで指定した未定義値は、'set gxout fwrite'でファイルに出力する際などに利用されます。Ver.2.0.a6以降で 'set undef' を使わないと、Ver.2.0.a6より前のGrADSとは挙動が異なるので要注意です。
計算したデータなどをnetCDFファイルで出力する(Ver.2.0a3-)
「set gxout fwrite」を使うときとは若干違う。先に全ての時間・格子情報を定義してから、出力する。
'open grads.ctl' 'set x 1 96' 'set y 1 96' 'set z 1 38' 'set t 1 10' 'define temp = t' 'set sdfwrite grads.nc' 'sdfwrite temp'
一度「define」しないといけないので(netCDFの定義の都合上)、注意。デフォルトは、1ファイル1変数しかできない。「set sdfattr」で設定できるようです。
できたファイルは「ncdump -h grads.nc」などで、簡易的に内容をチェックすると良い。
格子サイズが小さい場合（e.g.200m格子）に、格子点でのデータが表示されない問題
- 問題: 領域全体の描画は可能だが、各格子点での値や水平での平均値を表示しようとすると、次のようなエラーが出て、表示できない場合があります。
'open grads.ctl'
'set x 1 96'
Data Request Error: Invalid grid coordinates
World coordinates convert to non-integer grid coordinates
Variable = u Dimension = 0
Error ocurred at column 1
DISPLAY error: Invalid expression
Expression = u- 原因: 格子サイズが小さい場合、位置情報（緯度・経度の値）が桁落ちによって、隣の格子の位置と区別できなくなるため。
- 回避策
- 位置情報の有効桁数を小数点以下に多く確保するため、基準経度・緯度をゼロに近い値にする。
- （不便を承知ならば）緯度・経度座標を捨て、距離をkmで表示し、軸の数字を後から緯度・経度に書き換える。
バグ or 仕様：描画回数の限界
- v1.9b4で発生確認、v1.8SL11では発生せず。
- dコマンドを連続して57112回以上実行すると、セグメンテーション違反でGrADSが落ちる。
- ただし、clearやreinitを行うと回避できる。
- 描画の度にset ccolorで色を指定しても回避可能。
- おそらく色の設定の不具合がバグの原因。
- 以下のスクリプト参照
- 落ちる例（実用的ではありませんが…）
i = 1 while( i < 100000 ) 'd u' i = i + 1 endwhile
- 落ちない例
i = 1 while( i < 100000 ) 'set ccolor 1' 'd u' i = i + 1 endwhile
"Out of buffer space"が出る場合
描画に時間がかかるような図（高解像度の図など）を描くときに表示される場合があります。これが出ても描画する分には問題ありませんが、ファイルに保存をする場合は図が切れてしまいます。このような場合、図を分けて保存するなどの対処が必要でしょう。
ソースコードインストールができる人限定になりますが、より根本的な解決方法として、ソースコードのバッファ値を増やす方法があります。grads-2.0.a7.1 の場合、gxmeta.c の 249 or 250 という数字を大きくすることでエラーが出なくなります。たとえば 249 を 499 に、250 を 500 にするのが一案でしょう。なお、メモリを大量に食いつぶす恐れがありますので、物理メモリが小さなマシンでは注意が必要でしょう。
番号一覧
色
| 番号 | 名前 | 色 |
|---|---|---|
| 0 | background |  |
| 1 | foreground |  |
| 2 | red |  |
| 3 | green |  |
| 4 | dark blue |  |
| 5 | light blue |  |
| 6 | magenta |  |
| 7 | yellow |  |
| 8 | orange |  |
| 9 | purple |  |
| 10 | yellow/green |  |
| 11 | medium blue |  |
| 12 | dark yellow |  |
| 13 | aqua |  |
| 14 | dark purple |  |
| 15 | gray |  |
マークの種類
| 番号 | マーク |
|---|---|
| 0 | (なし) |
| 1 | ＋ |
| 2 | ○ |
| 3 | ● |
| 4 | □ |
| 5 | ■ |
| 6 | × |
| 7 | ◇ |
| 8 | △ |
| 9 | ▲ |
| 10 | ○｜ |
| 11 | ●｜ |
線種
| 線種番号 | 線種 |
|---|---|
| 0 | (なし) |
| 1 | 実線 |
| 2 | 長破線 |
| 3 | 短破線 |
| 4 | 長破線+短破線 |
| 5 | 点線 |
| 6 | 一点鎖線 |
| 7 | 二点鎖線 |
