making grid file from ascii xyz files - アスキー水深データからグリッドをつくる
使うコマンド xyz2grd, nearneighbor, blockmedian, surface
- ランダムなアスキーデータからグリッドデータファイルを作成するにはいくつかの方法がある。下図は、xyz2grd, nearneighbor, surfaceを使ってグリッド化した時の結果とスクリプトである。
- いずれのコマンドも-fgオプションを付けて明示的にgeographic gridであることを指定するほうがいい。このオプションがなくてもregionを緯度経度で指定し,intervalの単位を分などに指定することはできるので,一見普通にグリッドファイルができてきて,作図も問題ない(なので私も長いこと-fgはつけていなかった)。ただ,その場合は,ヘッダーに[Cartesian grid]と書き込まれてしまい,一部のグリッドを利用するコマンド(例えば後述のgrd2kml)ではエラーが出る
xyz2grd
- xyz2grdでは補間は行われず、純粋に元データのあるところだけに値が入り、他はNaN(この図ではNaNは灰色にしてある)となる。
- 領
域全体にデータが均質にあるならば、適正なグリッド間隔を指定するとこの例のようなスダレ状の図にはならない。しかし、一般に海域観測データは不均質に存
在し、細かめのグリッドを作るとこのようなスカスカの図ができてします。グリッド間隔(-I)を大きく取ると、マス目が大きくなるので画面中央部はすべて
のマス目にデータが入って埋め尽くされるが、一方で元々持っていた詳細情報が平滑化されて失われる。
- 従って、xyz2grdは観測データをいきなりグリッド化する際にはお勧めではない。が、例えば既にどこかでグリッド化されたデータをアスキーで入手した、といった場合は時々あるので、そういう場合はこのコマンドが役に立つ。
nearneighbor
- 観測
データに対してはxyz2grdはあまり使えないので、通常はnearneighborかsurfaceを利用する。nearneighborは、各方位
について指定した探索半径内で最も近くにあるデータを集め、それらの加重平均でグリッド化を行う。データ量が十分にある場合は、この方法がおすすめであ
る。
- nearneighbor
の良い点としては、欠損グリッドがあった場合に周囲のデータ(neighbor)から補間値で埋めてくれるが、neighborが近くにいない場合は無理
に埋めることはしない、というところである。結果として、図のように測深機のカバーしていない場所はNaNになり、エラーデータなどで小さな欠損が生じて
いるような場所は補間値で埋められたほどよい図ができる。
surface
- surfaceはいわば離散的な実測値の上に布をかぶせて全体を補間してしまう方法である。つまり図のように領域のはじっこまで果てしなく補間作業を続けてしまい、測線のカバーしていない(観測値で拘束できていない)場所にも値をつくる。データが本来ない場所を隠すには、グリッドマスクという作業が必要になる。
- た
だし、作図以外の目的でNaNを含まない完全なグリッドデータが必要になる場合がある。典型的な例は、重磁力データの解析で2次元FFTを使う場合な
どである。このような場合は、surfaceでとりあえずNaNなしデータを作成して計算を行い、表示するときにはグリッドマスクして元データのない場所
を隠す必要がある。
bash スクリプト例 xyz2grd
# make grid file from xyz files using xyz2grd
# parameter setting
region=137:40/137:55/16:40/17 # grid region east/west/south/north
interval=0.04m
# grid interval
xyzfile=area.xyz
# input xyz file
grdfile=area.xyzgrd
# output grid file
#
gmt xyz2grd $xyzfile -R$region -I$interval -G$grdfile -fg -V
bash スクリプト例 nearneighbor
# make grid file from xyz files using nearneighbor
# parameter setting
region=137:40/137:55/16:40/17 # grid region east/west/south/north
interval=0.04m
# grid interval
radius=0.5k
# grid search radius
xyzfile=area.xyz
# input xyz file
grdfile=area.xyzgrd
# output grid file
#
gmt nearneighbor $xyzfile -R$region -I$interval -S$radius -G$grdfile -fg -V
bash スクリプト例 surface
# make grid file from xyz files using surface
# parameter setting
region=137:40/137:55/16:40/17 # grid region east/west/south/north
interval=0.04m
# grid interval
tension=0.65
# tension factor
xyzfile=area.xyz
# input xyz file
blkfile=area.blk
# block file
grdfile=area.grd
# output grid file
#
gmt blockmedian $xyzfile -R$region -I$interval -V > $blkfile
gmt surface $blkfile -R$region -I$interval -T$tension -G$grdfile -fg -V