GrADSのTips

[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と入力する必要があります。これが面倒でしたら、予めリンクやエイリアスを設定するとよいでしょう。

[mydir@hosts] grads -b

引数-bを付けて起動することで、グラフィックス表示以外は通常通り使用することができます。windowsからTeraTermで接続してGrADSを使う場合に役立ちます。描画処理がない分、処理が高速になることも期待できます。

grads -c コマンド

• (例)GrADS起動直後にhelpを実行する。

  [mydir@hosts] grads -c help
  

!コマンド

• (例)GrADSを起動したままemacsでu.ctlを編集する

  ga-> !emacs u.ctl &
  ga-> (emacs実行中も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
  

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 ファイル名

とすると見られます。

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
  

64-bit機のGrADSの場合、「sdfopen」で開くと、まれに「289.615」という値を持つ箇所が、欠損値または「-32767.0」として認識される。~ 32-bit機のGrADSで開くか、Ctlファイルを作ってnetCDFを開くようにすると、解決できる。

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

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
  

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

あとは、2通りほど書き方がある。

xanim -pause -grfill var    ＃ マウスクリックでコマ送り（かつgrfillで表示）
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)
  

set grads on/off

set timelab on/off

set imprun スクリプト名/off

set background 色番号

• (例)背景色を赤に変更する(そんな人いない？)

  ga-> set background 2
  

• 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の方が適切であると記述されています（一応どちらでもうまくいくと思いますが）。

set datawarn on/off

set frame on/off/circle

circleは北極中心図などのときに有効です。

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
  

  • 描画例 (拡大版: latlon.gif)

    

set mproj scaled
set mproj off

scaledは地図あり、offは地図なし（緯度経度のラベルもなし）です。

• (例) 変数Tを縦横比を無視した直交座標で描く

  set mproj scaled
  set lat -90 90
  set lon 0 360
  d T
  

  • 描画例 (拡大版: scaled.gif)

    

set mproj robinson

lat, lon はそれぞれ[-90:90], [-180:180]のみ対応しています。

• (例) 変数Tをロビンソン図法で描く

  set mproj robinson
  set lat -90 90
  set lon -180 180
  d T
  

  • 描画例 (拡大版: robinson.gif)

    

set mproj mollweide

• (例) 変数Tをモルワイデ図法で描く

  set mproj mollweide
  set lat -90 90
  set lon 0 360
  d T
  

  • 描画例 (拡大版: mollweide.gif)

    

set mproj orthogr

「set lat -90 90」およびlonの範囲を180度にする必要があります。

• (例) 正投影図法で日本中心のTを描く

  set mproj orthogr
  set lat -90 90
  set lon 45 225
  d T
  

  • 描画例 (拡大版: orthogr.gif)

    

• (応用例)地球はまわる
  • 重い (1.6MB) ので要注意: orthogr_anim.gif;
  • GrADSスクリプト?、およびconvertコマンドを使って作成。

• 北極中心

  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
  

  • 描画例 (拡大版: nps2.gif)

    

• (例2) 変数Tを南極中心に描く

  set mproj sps
  set lat -90 0
  set lon 0 360
  d T
  

  • 描画例 (拡大版: sps.gif)

    

• (例3) 日本中心で±90度のTを描く

  set mproj nps
  set lat 0 90
  set lon 45 225
  d T
  

  • 描画例 (拡大版: nps1.gif)

    

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
  

set xlpos 位置 側

「位置」にはx軸からのずれ、「側」にはt(上側)、b(下側)のどちらかを指定します。

• (例)x軸のラベルを、下側の通常より0.5下に描画する

  ga-> set xlpos -0.5 b
  

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"を行なうと解除されます。

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))
  

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のインストール?を参照。

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 線種番号

• 線種番号一覧

set vrange 最小値 最大値

• (例)-3〜3の値の範囲を指定

  ga-> set vrange -3 3
  

set coslat on/off

高緯度に行くほど単位緯度あたりの軸の長さが短くなります。全球平均値の緯度毎の寄与を見る際に有効かも。

set missconn on/off

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で塗りつぶされます。

• 色番号一覧

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
  

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で指定した座標を中心に文字列が配置されます。太さは1-12の整数で指定します（デフォルト値は経験的に3）。回転角(度)を指定した場合は、座標の基準を中心に時計回りに回転します。

• (例)色番号2(赤)を指定

  ga-> set string 2
  

• (例)色番号2(赤)、座標の基準center、太さ1を指定

  ga-> set string 2 c 1
  

• (例)色番号2(赤)、座標の基準center、太さ1、回転角90度を指定

  ga-> set string 2 c 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.

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
  

`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
  

• 色番号一覧

多角形を使ったサンプル図（作成：大気海洋変動観測センター T.T.さんより）

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日として扱われます。

長い期間のデータを扱う場合、一つのファイルにまとめるとファイルサイズが大きくなることがあります。このようなとき、ファイルを分割すると扱いが楽になります。例えばデータを以下のように、月毎に分けることを考えます。

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桁で置き換えられます。この他にも以下のような表記があります (よく分からないものには？が付いています…)。

また、Version2.0aより、アンサンブルメンバーに関してもテンプレートで読み込めることが可能となりました。

なお、

/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というのがあるそうですが、誰か使っている人いますか？。

コントロールファイルの「OPTIONS」に「yrev」または「zrev」を加える。

二通りの設定方法があります。

XDEF 格子数 LINEAR 開始座標 増分
XDEF 格子数 LEVELS 座標1 座標2 ...

LINEARは等間隔格子、LEVELSは不等間隔格子の指定に適しています。基本的に単位は「度」ですが、kmなど他の単位で記述することも可能です。その場合、描画の前に

ga-> set mproj off

と指定するとよいでしょう。ただし地図を描きたい場合は、自前で地図データを用意して重ね描きする必要があります。

XDEFと基本的に同じです。

YDEF 格子数 LINEAR 開始座標 増分
YDEF 格子数 LEVELS 座標1 座標2 ...

異なるコントロールファイルのデータを比較したい場合(差・和を取りたいなど)、平面の場合ならX次元・Y次元の設定を、鉛直の場合ならX次元・Z次元、または、Y次元・Z次元の設定を、同じにしておく必要があるので注意!

TDEF タイムステップ数 LINEAR 開始時刻 時間増分

X次元などと違い、T次元は常にLINEAR（等間隔タイムステップ）です。開始時刻は

hh:mmZddmmmyyyy

もしくは簡単に、

ddmmmyyyy

等と指定します。ただし

時間増分は、'1yr', '2mn' などのように、「2桁以内の整数 + 単位」と表現します。単位は以下が使えます。

変数を列記する部分

変数名 層の数 unit 説明

のunitフィールドに指定します。指定方法は以下の通り（空白を入れないように注意！）。

• (例) 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

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 を使用してください。

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

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

前者（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

マッピングファイルを作るのに必要な最低限の情報は、

変換後の格子（緯度経度座標）における、変換前格子のグリッド番号

です。簡単のため、以下のような座標で考えてみましょう。

  |           |           |
(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組格納されています。

1. 各緯度経度座標におけるi座標値
2. 各緯度経度座標におけるj座標値
3. 各緯度経度座標における風速の回転の大きさ

   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
    

hcurl(東西風, 南北風)

東西風と南北風の変数から、水平渦度を計算します。

• (例)水平渦度を計算して表示

  ga-> a=hcurl(u,v)
  ga-> d a
  

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格子のデータを使用します。

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
  

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)

次に、環境変数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個で構成されています。各浮動小数点の意味は以下の通り。

1. 関数が呼ばれたときの引数の数
2. 常に0。関数データ交換ファイルのフォーマットを表す。(0以外の場合はエラーらしい)
3. これ以降は現在未定義(〜20)

引数部には、UDFで指定した引数の順番にデータが入っています。引数の種類（expr, value, char）によって書式が違います。

exprを指定した場合、

• 第1レコード(4bytes浮動小数点×20個)
  1. UNDEFの値
  2. i次元の種類を表す値(idim)
     • -1:なし
     • 0:X次元(経度)
     • 1:Y次元(緯度)
     • 2:Z次元(高度)
     • 3:T次元(時間)
  3. j次元の種類を表す値(jdim: 値の意味はitypeと同じ)
  4. i次元の要素の数(isiz)
  5. j次元の要素の数(jsiz)
  6. i次元が線形か非線形か(ilinear: 0なら非線形)
  7. j次元が線形か非線形か(jlinear: 0なら非線形)
  8. i次元の開始値(istrt)、ただしi次元が線形の場合
  9. i次元の増分(iincr)、ただしi次元が線形の場合
  10. j次元の開始値(jstrt)、ただしj次元が線形の場合
  11. j次元の増分(jincr)、ただしj次元が線形の場合
  12. 時間次元の開始時刻(年)、ただしどちらかの次元が時間の場合
  13. 時間次元の開始時刻(月)、ただしどちらかの次元が時間の場合
  14. 時間次元の開始時刻(日)、ただしどちらかの次元が時間の場合
  15. 時間次元の開始時刻(時)、ただしどちらかの次元が時間の場合
  16. 時間次元の開始時刻(分)、ただしどちらかの次元が時間の場合
  17. 時間次元の増分(分)、ただしどちらかの次元が時間の場合
  18. 時間次元の増分(月)、ただしどちらかの次元が時間の場合
  19. これ以降は現在未定義(〜20)
• 第2レコード
  1. グリッドデータ(isiz×jsiz個の4bytes浮動小数点)
• 第3レコード
  1. i次元の座標値(isiz個の4bytes浮動小数点)
     • i次元が存在しない(idim=-1)場合、このレコードは存在しない。
     • i次元が存在する(idimが-1ではない)場合、線形格子(ilinear=1)であってもこのレコードは存在する。
• 第4レコード
  1. j次元の座標値(jsiz個の4bytes浮動小数点)
     • j次元が存在しない(jdim=-1)場合、このレコードは存在しない。
     • j次元が存在する(jdimが-1ではない)場合、線形格子(jlinear=1)であってもこのレコードは存在する。

valueを指定した場合、4bytes浮動小数点×1個。これがそのまま引数の値です。

charを指定した場合、80bytesの文字列。80bytesに満たない部分には空白で埋められ、80bytesを超える場合は捨てられます。

プログラムで処理した結果をGrADSへ戻すには、UDFT5行目に指定した関数結果ファイル(function result file)を介して行います。ファイルはヘッダレコードとグリッドレコードに分かれています。

ヘッダレコードは浮動小数点(4bytes)20個で構成されています。各浮動小数点の意味は以下の通り。

1. プログラムの戻り値。0で正常終了、0以外で異常終了。異常終了の場合はこれ以上のデータ読み込みは行われない。
2. 常に0。関数結果ファイルのフォーマットを表す。
3. これ以降は現在未定義(〜20)

グリッドレコードは基本的に、関数データ交換ファイルの引数部と同じ形式で記述します。ただし第3レコードと第4レコードは扱いが少し異なるので注意が必要です。

• 第3レコード
  1. i次元の座標値(isiz個の4bytes浮動小数点)
     • i次元が存在しない(idim=-1)場合、又は線形格子(ilinear=1)の場合、このレコードは存在しない。
     • i次元が存在し、かつ非線形格子である場合に限り、このレコードは存在する。
• 第4レコード
  1. j次元の座標値(jsiz個の4bytes浮動小数点)
     • j次元が存在しない(jdim=-1)場合、又は線形格子(jlinear=1)の場合、このレコードは存在しない。
     • j次元が存在し、かつ非線形格子である場合に限り、このレコードは存在する。

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の中身は(図は&ref(binary.png,noimg);)

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

「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' とするとよいでしょう。

set undef 未定義値
set undef file ファイル番号  // 指定したファイル番号の未定義値を使用
set undef dfile              // 現在のファイル番号の未定義値を使用

ここで指定した未定義値は、'set gxout fwrite'でファイルに出力する際などに利用されます。Ver.2.0.a6以降で 'set undef' を使わないと、Ver.2.0.a6より前のGrADSとは挙動が異なるので要注意です。

「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」などで、簡易的に内容をチェックすると良い。

• 問題: 領域全体の描画は可能だが、各格子点での値や水平での平均値を表示しようとすると、次のようなエラーが出て、表示できない場合があります。

  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
  

• 原因: 格子サイズが小さい場合、位置情報（緯度・経度の値）が桁落ちによって、隣の格子の位置と区別できなくなるため。
• 回避策
  1. 位置情報の有効桁数を小数点以下に多く確保するため、基準経度・緯度をゼロに近い値にする。
  2. （不便を承知ならば）緯度・経度座標を捨て、距離をkmで表示し、軸の数字を後から緯度・経度に書き換える。

• 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
    

描画に時間がかかるような図（高解像度の図など）を描くときに表示される場合があります。これが出ても描画する分には問題ありませんが、ファイルに保存をする場合は図が切れてしまいます。このような場合、図を分けて保存するなどの対処が必要でしょう。

ソースコードインストールができる人限定になりますが、より根本的な解決方法として、ソースコードのバッファ値を増やす方法があります。grads-2.0.a7.1 の場合、gxmeta.c の 249 or 250 という数字を大きくすることでエラーが出なくなります。たとえば 249 を 499 に、250 を 500 にするのが一案でしょう。なお、メモリを大量に食いつぶす恐れがありますので、物理メモリが小さなマシンでは注意が必要でしょう。

上の「"Out of buffer space"が出る場合」と同様の問題と考えていますが、GrADS 2.1.a2 ではメッセージが表示されずに図が崩れたり切れたりすることがあります。これもソースコードインストールできる人限定の対処法ですが、 gxmeta.c の "BWORKSZ" を増やす（例えば10倍など）とうまく行くことがあります。

また、ソースコードインストールができなくても、描画の前に

ga-> set antialias off

を設定するとうまくいく場合もあります。

;