2012年 6月 1日 掲載
今回は，GMT に用意されている地図データを使って， 陸域または海域をマスクする grdlandmask というコマンドを紹介しよう．
まずはサンプルデータの用意である． その54 で作った mountains.awk を使って， 地形データの grdファイル mountains.grd を作る．
$ awk -f mountains.awk > mountains.txt $ xyz2grd mountains.txt -Gmountains.grd -R138.5/139.2/34.5/35.5 -I0.002
このデータには，海域にも標高データ（0m）が入っているので， 描画の際には，pscoastコマンドの -Sオプションで海域を塗りつぶして対処していた．
$ makecpt -Cglobe > g.cpt $ grdimage mountains.grd -Cg.cpt -Jm10 -K > mountains.eps $ pscoast -J -R -Sblue -Df -O >> mountains.eps
grdlandmask と grdmath を使えば，元のデータから海域部分を除去することができる．
まずは grdlandmask でマスクファイルを作ってみる． 最小限のオプションは grdmask と同じで， -G（マスクファイル名），-R（範囲），-I（グリッド間隔） である．
$ grdlandmask -Gmask.grd -R138.5/139.2/34.5/35.5 -I0.002
grdmask と違って，マスクに使うのが GMT に用意されている地図データなので， データファイルを与える必要はない．
できあがったマスクファイルを grdimage で描画してみよう．
$ makecpt -Cpolar -Z > pl.cpt $ grdimage mask.grd -Jm10 -Cpl.cpt > mask.eps
一応，伊豆半島だとはわかるが，かなり荒っぽい地図だ． 地図の解像度を変更するには，pscoast同様 -Dオプションを使う． -Dオプションで指定できる解像度は，高い方から順に，
である．最高の解像度 -Df を指定して作りなおしてみよう．
$ grdlandmask -Gmask.grd -R138.5/139.2/34.5/35.5 -I0.002 -Df $ grdimage mask.grd -Jm10 -Cpl.cpt > mask.eps
うん，いい感じだ．
このマスクファイルでは，陸域が 1（赤），海域が 0（白）になっている． 海陸の出力値の変更は，grdmask同様 -Nオプションで行う． -N には，海域の値，陸域の値を「/」で区切って指定する． 今回は海域をマスクしたいので，-NNaN/0 としてみた．
$ grdlandmask -Gmask.grd -R138.5/139.2/34.5/35.5 -I0.002 -Df -NNaN/0 $ grdimage mask.grd -Jm10 -Cpl.cpt > mask.eps
海域が NaN（グレー），陸域が 0（白）のマスクファイルができあがった． このマスクファイル mask.grd を使って，地形データ mountains.grd の海域をマスクする．
$ grdmath mountains.grd mask.grd OR = land.grd $ grdimage land.grd -Cg.cpt -Jm10 > land.eps
海域のデータの抜けた，きれいな地形データができあがった． このデータを使うと，下のような図も作れる．
$ grdview land.grd -Jm10 -Jz0.001 -E135/45 -Qs -Cg.cpt -Ba10m > land.eps
この例では，-Bオプションに「a10m」と書いて， 目盛りの間隔を「10分」に指定している．「m」が「分」の単位の指定である． また，できあがった図では，目盛りの数値が「度分」で描画されている． このあたりについては，そのうちゆっくり解説することにしよう．
[ KDOホームにもどる] [ いちからはじめる GMT インデックスにもどる ]
(C) 2012 KDO