*** GMT ちょっといい話し *** last updated : Nov 15, 1997 written by : M.HAYASAKI !! 注 !! 以下に書かれていることは誰にでも教えて構いませんし、この情報を どこに勝手に流そうと自由です。 ただし、ここに書いてあることに間違いがあったとしても私は一切責任を 持ちません。また、「これについて書いてくれ」という要求もお断りします。 私の時間が空いた時だけ追加・修正を行なうのみです。 なお、「こういうのを書いたから追加してくれ」といって電子化した状態 （すなわちテキストファイル）で私に渡す場合はこの限りでは ありません。適宜追加していきます。 はやさき まさみつ ------------------------------------------------------------------------- ・欠測データの処理(No.001) 測候所のデータの中には欠測値がつきものです。 そのようなデータを使ってグラフを書く時に、実際には存在しないデータを 書かせていることがありませんか？ もちろん、データを内挿すれば線を引くことは簡単ですが、本来存在しない データは書かせるべきではありません。 そんな時には以下のようにしてみるといいでしょう。 たとえば入力データ input.dat に欠測値が -999 で入っているとします。 set ifile=input.dat set ofile=out.ps awk '$2!=-999 {print $1,$2} $2==-999 {print ">"}' $ifile >! tmp uniq tmp tmp2 psxy tmp2 -R -JX -Ba100f50/a10f10g10WESn -M -W4 -K >! $ofile psxy tmp2 -R -JX -Ba100f50/a10f10g10WESn -M -Sc0.06 -L -W4 -O -K >> $ofile ここでやっていることを言葉で説明すると、 「入力データを読み込み、欠測（-999）であれば不等号記号（>）を 出力する。それ以外はもとの数値を出力。 その後、2つ以上連続する > 記号を uniq コマンドにより消去し、 psxy に -M オプションをつけることで欠測の部分は何も書かない。」 となります。 なお、この例では -R, -JX オプションに数値が設定されていませんが、 そのへんは自分の使う値を適当に放り込んで下さい。 今回はGMTというより、UNIXツールを使った「ちょっといい話し」でしたね。 references : "uniq" on-line manual awk 本いろいろ --------------------------------------------------------------------------- ・記号・パターンの出力(No.002) pstext を使って文字を書き込む際には、様々なフォント（書体）を使うことが 出来ます。 ただし出力されるかどうかはわかりません。マシン自体にフォントがなければ 別のフォントで代用されます。 フォントのサンプルとその対応する番号については font.gmt、 font2.gmt を 参照して下さい。 今回の本題はもっと細かい話しで、主に記号の出力に関して書きます。 主なサンプルは character.gmt を参照してもらうことにしましょう。 このサンプルには気象学でよく使われる 「℃」のようなもの、指数などの 文字飾りを中心に書いておきます。 今日は疲れたのでこのへんでやめときます。 references : GMT マニュアル p.15 and 61 --------------------------------------------------------------------------- ・効率の良いプログラミング・デバッギング(No.003) 「GMT は遅い！」という人がいますが、確かに遅い部分があるのは認めます。 しかし、自分自身で少しでも早く終らせるようにしようという努力をしていますか？ ちょっと努力すれば、最終出力に至るまで修正を何度も加えてやり直す場合でも 格段に早く終らせることが可能です。 例えば、以下のようなスクリプトがあったとします： （「中略」の部分にはGMTフォーマットにデータを切り出すためのFortranがある） なお、左の数字は行番号を表します。 1 #!/bin/csh -f 2 cd /data/work/hayasaki/tmp 3 #goto GMT 4 cat >! prog.f << EOF 5 PROGRAM MAIN （中略） 45 STOP 46 END 47 EOF 48 set dpath=/data/work/hayasaki 49 f77 -o DIVTEST prog.f 50 ln -s EC200velp.89jan fort.11 51 ln -s EC200wind.89jan fort.50 52 nice +10 DIVTEST # <= set priority "10" (recommended value) 53 unlink fort.11 ; unlink fort.50 54 rm prog.f DIVTEST 55 GMT: 56 set range=30/290/-40.0/60.0 # mapping range(world size) 57 # 58 set ofile=velp_test.ps # GMT output file 59 set velp=fort.21 # $3=velp 60 set divw=fort.31 # $3=divu,$4=divv 61 ### GMT mapping ### 62 pstext -R$range -JM8.0 -N -X1.3 -Y2.0 -K << END >! $ofile 63 160.0 70.0 28 0 5 2 Velocity potential (200hPa) 64 160.0 64.0 20 0 7 2 (Kinematical method) 65 END 66 psbasemap -R -JM -Ba30g30/a30g15WESn -O -K >> $ofile 67 awk '{print $1,$2,($3)/1000000}' $velp |xyz2grd -Gt.cdf -R -I2.5 68 awk '{print $1,$2,$3}' $divw | xyz2grd -Gu.cdf -R -I5/5 69 awk '{print $1,$2,$4}' $divw | xyz2grd -Gv.cdf -R -I5/5 70 pscoast -R -JM -A3000 -Dc -Ba30g30/a30g15WESn -W2 -O -K >> $ofile 71 grdcontour t.cdf -R -JM -C1 -A2f7 -L0/30 -G1.5/10 -O -K >> $ofile 72 grdcontour t.cdf -R -JM -C1 -A2f7 -L-30/0 -Wc1ta -Wa3t8_8:0 -G1.5/10 -O -K >> $ofile 73 psxy -R -JM -Sv0.02/0.03/0.04 -L -G0 -N -W1 -O -K << EOF >> $ofile 74 300 40 90 0.25 75 EOF 76 pstext -R -JM -N -O -K << EOF >> $ofile 77 303 39 20 0 4 1 5 m/s 78 EOF 79 grdvector u.cdf v.cdf -JM -R -G0 -Q0.008/0.01/0.01 -S20 -W -T -O >> $ofile 80 rm fort.21 fort.31 *.cdf 81 exit このスクリプトは以下の順序で作業がなされています。 I: データの切りだし（Fortran） II: GMT 開始（55行以降；最初は変数の設定、解説は略） i) psbasemap で枠を書く（66行） ii) awk でデータの切りだし => データをCDF形式に変換（67〜69行） iii) pscoast で地図を書く（70行） iv) grdcontour で等値線を引く（71〜72行） v) 判例を記入（73〜78行） vi) 風ベクトルの記入（79行） これで作業をしていると、見栄えのいい絵を１枚作るのに時間がかかりすぎます。 そんなアホなことをやってる人は、今すぐおやめなさい。時間のムダです。 では、どうすればいいのかといえば、 「中間ファイルを残しておく」 というのがいちばん良い方法です。 ＃ でも最終的に消すことは忘れないこと つまり、GMT用に切り出したデータは最初のうち（最終出力より前の段階） は消さない、ということです。この例でいえば、80行目をコメントアウトして 走らせましょう。 さらにいえば、データの切りだしは１度やればいいのですから、２回目からは 3行目の goto GMT のコメントアウトを外せば、すぐにGMTから始めてくれます。 また、ASCIIデータをCDF形式に変換する部分（67〜69行）ですが、 これも一度やってしまえば（データ間隔を変更しない限り）必要ないので、 コメントしておけばさらにスピードが上がります。 pscoast でも -Dc オプションや -A オプションなどを併用すればさらに 楽になります。 中間ファイルが残っていれば、もし作図に間違いがあった場合もどこで 間違えたのかがわかるので、デバッギング（バグとり）も早いでしょう。 このような努力をすれば、遅いといわれる GMT も、それほど負担に 感じません。何よりも余計な計算をさせないのですから、マシンへの 負荷も軽くなります。 みんなが実行するようになれば、使用状況も改善されるでしょう。 ＊＊余談＊＊ なお、わたしは ASCII から CDF 形式への変換には surface を使いません。 surface を使えばデータを内挿してくれますが、本来存在しないデータを 作りだしてまで絵を書くことがはたして必要でしょうか？ ＃ 第一、surface はクソ重い！ 特に、グローバルな解析をしている人は、データ間隔が2.5度程度もあれば、十分 滑らかな等値線が書かれるはずです。 客観解析データではなく不等間隔データ（ステーションデータなど）を 使っている場合でも、surface を使うよりは pscontour を使うほうが よっぽど賢いです(ただし Polar-stereographic では失敗する場合あり)。 各共用部屋にはマニュアルがあるのですから、各個人で自分が得になるような やり方を学んで下さい。 ただし、あまりハマりすぎないよーにしましょう。 論文に掲載する・発表に使用する、などでもない限り、 「自分が理解できればいい」 というのが原則です。図を書くのに時間をとられて解釈できないようでは 研究とはいえません。 ＃ これは自分に対する警告です。私がそうだったので ... --------------------------------------------------------------------------- 次回予告： 「タイトル未定」（Release 予定：未定）