いちからはじめる GMT その65 - grdmask -S

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

元々の値 cos(x)cos(y)

これを観測したデータには,中央部に大きな欠測域があったとしよう.

$ 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