いちからはじめる GMT その20 - surface

2009年 7月 1日 掲載


今回は,nearneighbor とならぶグリッドデータ作成コマンド surface を紹介するが….


前回,nearneighborコマンドを用いてグリッド補完を行い, grdimage でグラデーション図を作って,きれいな図ができたと書いたが, 実はコンター図を作るとそんなにきれいでもない.

$ awk 'BEGIN{for(i=0; i<=2500; i++){x=rand()*10-5; y=rand()*10-5; print x, y, (x^2-y^2)}}' > hpr.txt
$ nearneighbor hpr.txt -Ghpr2.grd -I0.2 -R-4/4/-4/4 -S1.0
$ grdcontour hpr2.grd -Jx1 -C1 -A5 -Ba1 > hpr.eps

双曲方物面

ガタガタである. nearneighborコマンドは,グリッド点に一番近いデータを補完に使うのだが, となり合うグリッドで,補完に使われるデータがちがうという場所で,不連続が生じてしまうのだ.

GMT には,nearneighbor のほかに,surface という,なめらかな grdファイルを作成するコマンドがある. どういうことをするのか man を読んでみよう. surface の man の 1行目には次のように書いてある.

surface - adjustable tension continuous curvature surface gridding algorithm

何のことかわからない….本文を読んでみよう.

surface reads  randomly-spaced (x,y,z) triples from standard input [or
xyz_file] and produces a binary grid file of gridded  values  z(x,y)  by
solving:

    (1 - T) * L (L (z)) + T * L (z) = 0

where T is a tension factor between 0 and 1, and L indicates the Lapla-
cian operator.  T = 0 gives the "minimum curvature" solution  which  is
equivalent  to  SuperMISP  and the ISM packages.  Minimum curvature can
cause undesired oscillations and false  local  maxima  or  minima  (See
Smith  and  Wessel,  1990),  and  you may wish to use T > 0 to suppress
these effects.  Experience suggests T ~ 0.25  usually  looks  good  for
potential field data and T should be larger (T ~ 0.35) for steep topog-
raphy data.  T = 1 gives a harmonic surface (no maxima  or  minima  are
possible  except  at  control data points). 

ええと…,さっぱりわからない…. わからないなりに使ってみよう. surfaceコマンドを単独で実行して,ずらっと出てくる usage を読んでみる.

$ surface
surface 4.2.1 - Adjustable tension continuous curvature surface gridding

usage: surface [xyz-file] -G<output_grdfile_name> -I<xinc>[u][=|+][/<yinc>[u][=|+]]
	-R<west>/<east>/<south>/<north>[r] [-A<aspect_ratio>] [-C<convergence_limit>] [-H[i][<nrec>]]
	[-Ll<limit>] [-Lu<limit>] [-N<n_iterations>] ] [-S<search_radius>[m|c]] [-T<tension>[i][b]]
	[-Q] [-V[l]] [-Z<over_relaxation_parameter>] [-:[i|o]] [-bi[s|S|d|D][<ncol>]] [-f[i|o]<colinfo>]

	surface will read from standard input or a single <xyz-file>.

	Required arguments to surface:
	-G sets output grid file name
	-I<xinc>[m|c|e|k|i|n|+][=][/<yinc>[m|c|e|k|i|n|+][=]]
	   Give increment and append unit (m)inute, se(c)ond, m(e)ter, (k)ilometer, m(i)les, (n)autical miles.
	   (Note: m,c,e,k,i,n only apply to geographic regions specified in degrees)
	   Append = to adjust the domain to fit the increment [Default adjusts increment to fit domain].
	   Alternatively, specify number of nodes by appending +. Then, the increments are calculated
	   from the given domain and grid-registration settings (see Appendix B for details).
	-R specifies the min/max coordinates of data region in user units.
	   Use dd:mm[:ss] format for regions given in degrees and minutes [and seconds].
	   Use [yyy[-mm[-dd]]]T[hh[:mm[:ss[.xxx]]]] format for time axes.
	   Append r if -R specifies the longitudes/latitudes of the lower left
	   and upper right corners of a rectangular area.
	   -Rg -Rd are accepted shorthands for -R0/360/-90/90 -R-180/180/-90/90

	OPTIONS:
           :
           :

最低限,-Gオプションで grdファイル名,-Iオプションでグリッド間隔, -Rオプションで範囲を指定すれば, 「adjustable tension continuous curvature surface gridding」 してくれるようである.では,さっそく,

$ surface hpr.txt -Ghpr3.grd -I0.2 -R-4/4/-4/4
$ grdcontour hpr3.grd -Jx1 -C1 -A5 -Ba1 > hpr2.eps

双曲方物面

とりあえず,きれいなコンター図ができた.

man の解説によると,tension factor T という値を変えることによって, 出来上がりがずいぶん違うようだ. T=0 で最小曲率解?,T=1 で調和表面?だそうだが,何のことやら…. 経験的に,ポテンシャル場には T〜0.25 を使うことを, 勾配の急な地形データには T〜0.35 を使うことをすすめる, と書いてあるらしい(英語が苦手なので,間違っているかもしれない).

式を見せられても,解説を読んでも,何のことかさっぱりわからないので, 実際に T を変えて作ったデータを描画して,違いを見てみることにする. T の指定は -Tオプションで行う.デフォルトでは T=0.

まずは,元データ作り.粗なデータセットを乱数で作る.

$ awk 'BEGIN{for(x=-3; x<=3; x+=2)for(y=-3; y<=3; y+=2)print x, y, rand()*20-10}' > rand.txt

できたデータをプロットしてみる.

$ echo -16 0 0 255 0 255 255 255 > cp2.cpt
$ echo 0 255 255 255 16 255 0 0 >> cp2.cpt
$ psxy rand.txt -R-4/4/-4/4 -Jx1 -Ba1 -Sc0.2 -Ccp2.cpt -W > rand.eps

乱数

これを nearneighbor で補完するとこうなる.

$ nearneighbor rand.txt -Grandnn.grd -I0.2 -R-4/4/-4/4 -S3.0
$ grdimage randnn.grd -Jx1 -Ccp2.cpt -Ba1 -K > randnn.eps
$ grdcontour randnn.grd -J -C1 -A5 -O >> randnn.eps

乱数

この図を基準にして,surface でできたグリッドデータを見てみよう.

まずは,tension factor T=0 の場合.

$ surface rand.txt -Granda.grd -I0.2 -R-4/4/-4/4
$ grdimage randa.grd -Jx1 -Ccp2.cpt -Ba1 -K > randa.eps
$ grdcontour randa.grd -J -C1 -A5 -O >> randa.eps

乱数

nearneighbor はデータに囲まれた領域の内挿しかしなかったが, surface は外挿もするようだ.

次に,おすすめの T=0.25 の場合.

$ surface rand.txt -Grandb.grd -I0.2 -R-4/4/-4/4 -T0.25
$ grdimage randb.grd -Jx1 -Ccp2.cpt -Ba1 -K > randb.eps
$ grdcontour randb.grd -J -C1 -A5 -O >> randb.eps

乱数

最後に,T=1.0 の場合を見てみる.

$ surface rand.txt -Grandc.grd -I0.2 -R-4/4/-4/4 -T1.0
$ grdimage randc.grd -Jx1 -Ccp2.cpt -Ba1 -K > randc.eps
$ grdcontour randc.grd -J -C1 -A5 -O >> randc.eps

乱数

このケースでは,T=1.0 はとがりすぎで論外である. T=0.0 はなめらかだけど,大げさすぎかな? 最適な T は,いろいろ試して決めるのがよいようである.


ちなみに,筆者は surface をあまり使わない. 結果はきれいだけど,よくわからないものを使うのはどうも…. ふだんは,もっぱら nearneighbor を使い, コンター図がガタガタになって困るときだけ使っている.


[ KDOホームにもどる ] [ いちからはじめる GMT インデックスにもどる ]


(C) 2009 KDO