2008年 11月 1日 掲載
お待たせしました．PS出力の合成をやります．
今回は，ばらつきのある観測データと，近似直線を描画してみる．
まずはデータ作りだが，適当に用意した観測データから， 近似直線を最小自乗法とかで計算するのは， awk 1行で，というわけにはいかない． そこで，近似直線は y=x と決めてしまい， 観測データは y=x に乱数でばらつきを加えて作ることにする．
次のようにして，観測データ obs.txt と，近似直線データ x.txt をつくる．
$ awk 'BEGIN{for(x=0;x<=8;x+=0.2) print x, x+(rand()-0.5)}' > obs.txt
$ awk 'BEGIN{for(x=0;x<=8;x+=8) print x, x}' > x.txt
それぞれを描画してみる．
$ psxy obs.txt -R0/8/0/8 -Jx1 -Bg1a1 -Sc0.2 > obs.eps
$ psxy x.txt -R0/8/0/8 -Jx1 -Bg1a1 > x.eps
ちなみに，直線を描画するのは，このシリーズでは，はじめてだったりする．
さて，この 2つの図を合成したい． 出力の合成は「>>」だったから，これでひとつのファイルにつないでみよう．
$ psxy obs.txt -R0/8/0/8 -Jx1 -Bg1a1 -Sc0.2 > obs_x.eps $ psxy x.txt -R0/8/0/8 -Jx1 -Bg1a1 >> obs_x.eps
すると，obs_x.eps は 1ページ目が観測データ，2ページ目が近似直線 …， という複数ページの EPSファイルになってしまった． しかも，白紙のページがあったり，前のページにもどれなかったり，何かへんである．
psxyコマンド（に限らず，GMT の PS出力コマンド）の出力を， 1枚の図に重ねて描画するためには，-Kオプションと -Oオプションを使う必要がある． 2つのオプションは man psxy では，以下のように説明されている．
PostScriptコードをさらに後ろに追加したいときには -Kオプションを， 前の図に重ねてプロットするときには -Oオプションをつけるのである．
今回の例では，obs.txt の出力に -Kオプションをつけ， 「>>」でつづける x.txt の出力に -Oオプションをつければよい． なお，-Bg1a1 は重複するので，片方だけに記述する．
$ psxy obs.txt -R0/8/0/8 -Jx1 -Bg1a1 -Sc0.2 -K > obs_x.eps $ psxy x.txt -R0/8/0/8 -Jx1 -O >> obs_x.eps
ほしい図を描画することができた．
2つの出力を合成するには，
$ psxy 〜 -K > foo.eps $ psxy 〜 -O >> foo.eps
だった．では，3つ以上の出力を合成するにはどうすればいいのかというと， 以下のようにすればよい．
$ psxy 〜 -K > bar.eps $ psxy 〜 -O -K >> bar.eps ： $ psxy 〜 -O -K >> bar.eps $ psxy 〜 -O >> bar.eps
最初に -Kオプションの出力，次から -O -K をつけた出力をつづけ， 最後に -Oオプションの出力でしめるのである．
以下に具体例を示す． obs.txt の座標に大きさのことなる円を重ねて描画して， 二重丸にしてみた（あまりいい例ではないが…）．
$ psxy obs.txt -R0/8/0/8 -Jx1 -Bg1a1 -Sc0.3 -K > obs_x2.eps $ psxy obs.txt -R0/8/0/8 -Jx1 -Sc0.15 -O -K >> obs_x2.eps $ psxy x.txt -R0/8/0/8 -Jx1 -O >> obs_x2.eps
ここで，-R，-Jオプションに同じパラメータを何度も打つのはめんどうだ，と感じたかと思う． 実は，2つ目以降の -R，-Jオプションのパラメータは省略できる． 以下のようにしても，同じ図ができあがる．
$ psxy obs.txt -R0/8/0/8 -Jx1 -Bg1a1 -Sc0.3 -K > obs_x2.eps $ psxy obs.txt -R -J -Sc0.15 -O -K >> obs_x2.eps $ psxy x.txt -R -J -O >> obs_x2.eps
1つ目の psxyコマンドが実行されたとき， カレントディレクトリに .gmtcommands4 というファイルが作られて， 指定した各パラメータが保存され，以降のコマンドのデフォルト値として参照されるのだ． 今回の例では，以下の内容の .gmtcommands4 が作られる．
# GMT common arguments shelf -Bg1a1 -Jx1 -R0/8/0/8 -jx1 EOF
パラメータを打ちなおす手間がはぶけるし，打ち間違いをする危険も回避できる．
おっと，忘れるところだった． その4 の最後に紹介した図↓の描き方を説明しなくては．
この図を描画するには，psxy で，曲線 y=exp(x) および，グリッド，左辺（W）・下辺（S）の目盛・ラベルを描き， psbasemap という，データを描画しない図の土台だけ描くコマンドで，右辺（E）・上辺（N）の目盛を追加すればよい．
$ awk 'BEGIN{for(x=-3;x<=3;x+=0.1)print x, exp(x)}' > exp.txt
$ psxy exp.txt -R-3/3/0/12 -Jx2/1 -Bg1a1:x:/g1a1:"y = exp(x)":WS -K > exp.eps
$ psbasemap -R -J -Ba1/a1EN -O >> exp.eps
[ KDOホームにもどる ] [ いちからはじめる GMT インデックスにもどる ]
(C) 2008 KDO