2012年 5月 1日 掲載
grdmask のつづきである． 今回は，surfaceなどにより補間された grdファイルの， 元データのない範囲をマスクする方法を紹介する．
まずは，サンプルデータを作る． 元々の値が cos(x)cos(y) の分布をしていて，
$ awk 'BEGIN{for(x=-5;x<=5;x+=0.25)for(y=-5;y<=5;y+=0.25)print x, y, cos(x)*cos(y)}' > cosx_cosy.txt
$ makecpt -Cpolar -Z > pl.cpt
$ psxy cosx_cosy.txt -Jx1 -R-5/5/-5/5 -Sc0.2 -Cpl.cpt -W1 > cosx_cosy.eps
これを観測したデータには，中央部に大きな欠測域があったとしよう．
$ awk 'BEGIN{for(x=-5;x<=5;x+=0.25)for(y=-5;y<=5;y+=0.25)if(x^2+y^2>4)print x, y, cos(x)*cos(y)}' > cosx_cosy.txt
$ psxy cosx_cosy.txt -Jx1 -R-5/5/-5/5 -Sc0.2 -Cpl.cpt -W1 > cosx_cosy.eps
この観測データを surface で補間すると， 当然中央部は元々の値と異なるものとなる．
$ surface cosx_cosy.txt -Gcosx_cosy.grd -R-5/5/-5/5 -I0.125 $ grdimage cosx_cosy.grd -Jx1 -Cpl.cpt > cosx_cosy.eps
そこで，補間された grdファイルから観測データのない範囲を取りのぞくことにする．
grdmask に -Sオプションをつけると，読み込んだ座標周辺を内側， それ以外を外側としたマスクファイルが作られる． -Sオプションには，座標からどこまでの範囲を内側とするかの半径の値をあたえる必要がある． 今回は，観測データの格子間隔が 0.25 なので，半径は 0.177 にしておこう （0.25 / √2 = 0.17677…）．
$ grdmask cosx_cosy.txt -Gmask.grd -R-5/5/-5/5 -I0.125 -S0.177
できあがったマスクファイル mask.grd を grdimage で見てみる．
$ grdimage mask.grd -Jx1 -Cpl.cpt > mask.eps
観測データのある範囲が 1（赤）で， 観測データのない中央部が 0（白）になっている． 今回は，観測データの「ない」範囲を NaN にしたいので， -Nオプションで，外側の値を NaN にする． 内側の値は NaN でなければ何でもいいが，-NNaN/0/0 としておこう．
$ grdmask cosx_cosy.txt -Gmask.grd -R-5/5/-5/5 -I0.125 -S0.177 -NNaN/0/0 $ grdimage mask.grd -Jx1 -Cpl.cpt > mask.eps
このマスクファイル mask.grd を使って， 補間された grdファイル cosx_cosy.grd をマスクすると，
$ grdmath cosx_cosy.grd mask.grd OR = masked_data.grd
期待通りの図ができあがる．
$ grdimage masked_data.grd -Jx1 -Cpl.cpt > masked_data.eps
[ KDOホームにもどる] [ いちからはじめる GMT インデックスにもどる ]
(C) 2012 KDO