2009年 10月 1日 掲載
前回までの 3回で，カラーパレットファイルを使ったカラフルな図を作成したが， どの色がどれくらいの値なのかを示す凡例がないと， 図を見てもらう人に対して不親切である． そこで，今回はカラーパレットのスケールを描画する psscaleコマンドを紹介する．
今回は，ストーリーの流れから psscaleコマンドを取り上げることにしたが， 動作にバグっぽいところがあったりして扱いにくく，筆者はいまいち使いこなせていない． GMT のコマンドでは解決できずに，出力した EPSファイルを直接編集したりもしているが， 実は GMT のコマンド，オプションで何とかできるかもしれないし， GMT のバージョンがちがうと，こういう作業は必要ないのかもしれない． 今回使った GMT のバージョンは 4.2.1 で， PSファイルではなく EPSファイルを出力するように設定してある．
まずは，psscale を単独で実行し，usage を見てみる．
$ psscale
psscale 4.2.1 - To create grayscale or colorscale for maps
usage: psscale -D<xpos/ypos/length/width>[h] [-A] [-C<cpt_file>] [-E[b|f][<length>]] [-B<params>] [-I[<max_intens>|<low_i>/<high_i>]
[-K] [-L[i][<gap>]] [-M] [-N<dpi>] [-O] [-P] [-Q] [-S] [-U[/<dx>/<dy>][<label>]] [-V] [-X[a|c|r]<x_shift>[u]]
[-Y[a|c|r]<x_shift>[u]] [-Z<zfile>] [-c<ncopies>]
-D set mid-point position and length/width for scale.
Give negative length to reverse the scalebar
Append h for horizontal scale
OPTIONS:
-A Place the desired annotations on the other side of the colorscale instead
-B Set scale annotation interval and label. Use y-label to set unit label
If no annotation interval is set it is taken from the cpt file
-C Color palette file. If not set, stdin is read.
:
:
-Dオプションで描画位置を決め， -Cオプションでカラーパレットを指定すればよさそうである （-C は必須オプションではないが，これをセットしない場合は， かわりに標準入力でデータを与える必要がある）．
では，さっそく前回も作った双曲放物面のグラフに， カラーバーを追加してみよう．
-Dオプションでは，中心の位置と，長さ，幅をセットする必要がある． 中心の位置って，どこからの位置だ？まあ，最初は (0, 0) にして様子を見てみるか． 長さは，グラフの大きさにあわせて 8，幅はバランスを考えて 1 かな．
-Cオプションでは，polar のカラーパレットファイルを使おう．
$ awk 'BEGIN{for(x=-4;x<=4;x+=0.1)for(y=-4;y<=4;y+=0.1)print x, y, x^2-y^2}' > hp.txt
$ xyz2grd hp.txt -Ghp.grd -I0.1 -R-4/4/-4/4
$ makecpt -Cpolar -T-16/16/16 -Z > pl.cpt
$ grdimage hp.grd -Jx1 -Cpl.cpt -Ba1 -K > hp.eps
$ psscale -D0/0/8/1 -Cpl.cpt -O >> hp.eps
縦長のカラーバーが描かれ， バーの左端中央が，グラデーション図の範囲の左下隅に一致するように配置された． psscale が基準とする座標は，グラフの中の座標ではなく，描画する「紙」の座標で， 「紙」の座標の原点は最初に描画した図の左下隅なのである．
描画する座標系がわかれば，適切な位置も計算できる． 縦長のバーなので，図の右 (10, 4) に配置することにしよう．
$ grdimage hp.grd -Jx1 -Cpl.cpt -Ba1 -K > hpb.eps $ psscale -D10/4/8/1 -Cpl.cpt -O >> hpb.eps
あれ？バーの位置はいいが，目盛が切れちゃってる…． どうも，出力された EPSファイルの %%BoundingBox の値に， 目盛の分が加算されていないようだ．
$ tail hpb.eps S 0 0 M 58 F0 (16) E /PSL_dimx exch def bby /PSL_dimy exch def 945 -142 M V -90 R PSL_dimx neg PSL_dimy 2 div neg G 58 F0 (16) Z U -90 R -118.11 0 T %%Trailer %%BoundingBox: 225 0 601 383 S -1181.1 -0 T 4.16667 4.16667 scale 0 A showpage end
正常な図を作るためには， %%BoundingBox の値を変更しなくてはならないが， GMT のコマンドでどうにかなるのか， いろいろ調べたり，試したりしたが，よくわからない． 結局，出力された EPSファイルを編集して解決することにした． EPSファイルをエディタで開いて書きかえるのが，手っ取り早い解決法ではあるが （EPSファイルはアスキーファイルだからね）， 図を何度も作り直したり，大量の図を出力したりする場合は， いちいち手作業で編集するなど，やってられない． そういう作業は，awkスクリプトを一行書いてやらせるのがいい． %%BoundingBox の行の 4つ目の数値に，80 をプラスしてみる．
$ awk '{if($1=="%%BoundingBox:" && $5>0){print $1, $2,$3, $4, $5+80}else{print $0}}' < hpb.eps > hpb2.eps
何とか，のぞみ通りの図をつくることができた．
-Dオプションに「h」をつけると，バーが水平になる． グラフの下に配置してみよう． バーの上端中央の位置を -Dオプションで指定してやる．
$ grdimage hp.grd -Jx1 -Cpl.cpt -Ba1 -K > hpc.eps $ psscale -D4/-1/8/1h -Cpl.cpt -O >> hpc.eps
今度は，目盛は表示されたが，下の余白が多すぎる． これも awk のスクリプトで %%BoundingBox の値を変更する． %%BoundingBox の行の 3つ目の数値から 100 をマイナスする．
$ awk '{if($1=="%%BoundingBox:" && $4>0){print $1, $2,$3, $4-100, $5}else{print $0}}' < hpc.eps > hpc2.eps
やれやれ，やっとまともな図になった．
以上で，psscale の使い方，というか使いこなす方法の解説を終えることにする． 出力された EPSファイルをいじるなんて，スマートなやり方ではないが， 今の筆者には，これがせいいっぱいである．
実は，psbasemap で大きめの白い図を最初に作っておく， という手もあるのだが，PNGファイルなどに変換したとき，透明じゃなくなるのと， 水平バーを下に描画したときのように，図を小さくするということができないので， 今回は紹介しなかった．
[ KDOホームにもどる ] [ いちからはじめる GMT インデックスにもどる ]
(C) 2009 KDO