2012年 4月 1日 掲載
これまで，2次元のグリッドファイル 「grdファイル」をあつかった描画をいくつか紹介してきたが， 実際のデータを使っていると， 不都合な部分（意図せずに補間された部分など）をかくしたい， という場面に出会うこともときどきある． そういうときに使うのが grdmask と grdmath で， これらのコマンドを用いると， grdファイルから特定の範囲のデータを取り除くことができる． もう少しくわしく手順を書くと， grdmask で「マスクする」範囲を書きこんだ grdファイルを作り， grdmath でその範囲のデータを取り除くのである．
では，grdmask から解説していくことにしよう． grdmask でマスクする範囲を指定する方法はいくつかあるが， 今回は多角形の頂点を読み込む方法を紹介する． まずは，多角形のデータを作ろう． AWKスクリプトで，星形の頂点の座標を star.txt に書き出した．
$ awk 'BEGIN{for(i=0;i<10;i++){d=2*3.1415/10*i;r=1+i%2*1.5;print r*cos(d), r*sin(d)}}' > star.txt
$ cat star.txt
1 0
2.02257 1.46943
0.309052 0.951045
-0.77241 2.37768
-0.808973 0.587845
-2.5 0.000231634
-0.809082 -0.587695
-0.772851 -2.37754
0.308876 -0.951102
2.0223 -1.4698
このデータを psxy -G で描画すると，以下のようになる．
$ psxy star.txt -Jx1 -R-6/6/-6/6 -Ggray > star.eps
この星形の範囲をマスクする grdファイル（以後，マスクファイルと呼ぶ）を grdmask で作るのだが， 最初は最小限のオプションでためしてみよう．
$ grdmask star.txt -Gstar.grd -R-6/6/-6/6 -I0.125
-Gオプションでマスクファイルのファイル名， -Rオプションでマスクファイルのデータの範囲， -Iオプションでグリッド間隔を指定する． できあがったマスクファイル star.grd を ncdumpコマンドで見てみると， 0 と 1 がズラッと並んでいる．
$ ncdump star.grd
netcdf star {
dimensions:
x = 97 ;
y = 97 ;
variables:
double x(x) ;
x:long_name = "x" ;
x:actual_range = -6., 6. ;
double y(y) ;
y:long_name = "y" ;
y:actual_range = -6., 6. ;
float z(y, x) ;
z:long_name = "z" ;
z:_FillValue = NaNf ;
z:actual_range = 0., 1. ;
// global attributes:
:Conventions = "COARDS/CF-1.0" ;
:title = "star.grd" ;
:history = "grdmask -R-6/6/-6/6 -I0.125 star.txt -Gstar.grd" ;
:GMT_version = "4.5.7 [64-bit]" ;
data:
x = -6, -5.875, -5.75, -5.625, -5.5, -5.375, -5.25, -5.125, -5, -4.875,
-4.75, -4.625, -4.5, -4.375, -4.25, -4.125, -4, -3.875, -3.75, -3.625,
-3.5, -3.375, -3.25, -3.125, -3, -2.875, -2.75, -2.625, -2.5, -2.375,
-2.25, -2.125, -2, -1.875, -1.75, -1.625, -1.5, -1.375, -1.25, -1.125,
-1, -0.875, -0.75, -0.625, -0.5, -0.375, -0.25, -0.125, 0, 0.125, 0.25,
0.375, 0.5, 0.625, 0.75, 0.875, 1, 1.125, 1.25, 1.375, 1.5, 1.625, 1.75,
1.875, 2, 2.125, 2.25, 2.375, 2.5, 2.625, 2.75, 2.875, 3, 3.125, 3.25,
3.375, 3.5, 3.625, 3.75, 3.875, 4, 4.125, 4.25, 4.375, 4.5, 4.625, 4.75,
4.875, 5, 5.125, 5.25, 5.375, 5.5, 5.625, 5.75, 5.875, 6 ;
y = -6, -5.875, -5.75, -5.625, -5.5, -5.375, -5.25, -5.125, -5, -4.875,
-4.75, -4.625, -4.5, -4.375, -4.25, -4.125, -4, -3.875, -3.75, -3.625,
-3.5, -3.375, -3.25, -3.125, -3, -2.875, -2.75, -2.625, -2.5, -2.375,
-2.25, -2.125, -2, -1.875, -1.75, -1.625, -1.5, -1.375, -1.25, -1.125,
-1, -0.875, -0.75, -0.625, -0.5, -0.375, -0.25, -0.125, 0, 0.125, 0.25,
0.375, 0.5, 0.625, 0.75, 0.875, 1, 1.125, 1.25, 1.375, 1.5, 1.625, 1.75,
1.875, 2, 2.125, 2.25, 2.375, 2.5, 2.625, 2.75, 2.875, 3, 3.125, 3.25,
3.375, 3.5, 3.625, 3.75, 3.875, 4, 4.125, 4.25, 4.375, 4.5, 4.625, 4.75,
4.875, 5, 5.125, 5.25, 5.375, 5.5, 5.625, 5.75, 5.875, 6 ;
z =
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
:
:
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1,
:
:
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
}
これを grdimage で見てみると，
$ makecpt -Cpolar -Z > pl.cpt $ grdimage star.grd -Jx1 -Cpl.cpt > star.eps
星形の内側が 1（赤），外側が 0（白）になっているのがわかる．
さて，マスクファイルは grdmath を使って， データファイルからデータを取り除くのに使うわけだが， 0 と 1 の grdファイルでは，ちょっと都合がわるい． そこで，内側の値が「NaN」になるように，マスクファイルを作りなおすことにする． grdファイルに書きこむ値を変更するには，-Nオプションを使う． -N では，外側の値，境界の値，内側の値を「/」で区切って指定する． デフォルトでは「0/0/1」となっているのを，「0/0/NaN」に変更しよう．
$ grdmask star.txt -Gmask_star.grd -I0.125 -R-6/6/-6/6 -N0/0/NaN
今回は，外側と境界を 0 にしたが，NaN でなければ何でもいい． できあがったマスクファイル mask_star.grd を，また grdimage で見てみよう．
$ grdimage mask_star.grd -Jx1 -Cpl.cpt > mask_star.eps
星型の内側が NaN（グレー）になっているのが確認できた．
マスクファイルの用意ができたので， 次はマスクする対象となるデータファイルを作成する． マンネリでもうしわけないが，またガウス分布の山を AWK で作る．
$ awk 'BEGIN{for(x=-6;x<=6;x+=0.125)for(y=-6;y<=6;y+=0.125) print x,y,exp(-(x^2+y^2)/16)}' > gauss.txt
$ xyz2grd gauss.txt -Gdata_gauss.grd -R-6/6/-6/6 -I0.125
$ grdimage data_gauss.grd -Jx1 -Cpl.cpt > data.eps
データファイルから，マスクファイルの NaN の範囲のデータを取り除くには， grdmath で 2つの grdファイルの OR演算をする．
$ grdmath data_gauss.grd mask_star.grd OR = masked_data.grd
grdmath は grdファイルの演算を行うコマンドで， 演算式の書き方は逆ポーランド記法である． 今回の OR演算は，「A OR B」でなく「A B OR」と書く． 「=」は結果を出力する記号で，続けて出力ファイル名を記す． grdmath の OR演算は，man grdmath には，
OR 2 1 NaN if A or B == NaN, else A.
と書いてある． A，B のどちらかあるいは両方が NaN なら結果は NaN， どちらも NaN でなければ結果は「A」である．
| - | A==NaN | A!=NaN |
| B==NaN | NaN | NaN |
| B!=NaN | NaN | A |
今回の例では，マスクファイル mask_star.grd で NaN のグリッドは NaN， それ以外のグリッドはデータファイル data_gauss.grd の値となるので， マスクファイルでデータファイルがマスクされたことになるのである．
得られたファイル masked_data.grd は以下のようになっている．
$ grdimage masked_data.grd -Jx1 -Cpl.cpt > masked_data.eps
ちゃんと星形の内側がマスクされている． grdview で描画すると，データが抜けているのが，より一層はっきりとわかる．
$ grdview masked_data.grd -Jx1 -Jz5 -E60/30 -Cpl.cpt -Qs > masked_data.eps
[ KDOホームにもどる] [ いちからはじめる GMT インデックスにもどる ]
(C) 2012 KDO