いちからはじめる GMT その53 - grdimage -I, grdgradient

2011年 8月 1日 掲載

今回は，grdimage の出力に陰影をつけて， メリハリの効いた立体感のある図を描く方法を解説する．

まずは，サンプルデータを作る． 次のような AWKスクリプト crater.awk を書き， 直径 5，高さ 1 のクレーター状の地形データ crater.txt を作ってみた．

$ cat crater.awk
BEGIN{
  for(x = -5; x <= 5; x += 0.03125){
    for(y = -5; y <= 5; y += 0.03125){
      r = sqrt(x^2 + y^2)

      if(r < 1.5){
        z = 0
      }else if(r < 2.5){
        z = r - 1.5
      }else if(r < 3.5){
        z = 3.5 - r
      }else{
        z = 0
      }

      print x, y, z
    }
  }
}
$ awk -f crater.awk > crater.txt

この地形の断面は二等辺三角形になっている．

$ awk '$1==0{print $2, $3}' crater.txt > profile.txt
$ psxy profile.txt -Jx1 -R-5/5/-2/2 -Ba2 > profile.eps

断面図

このデータを使って，grdimage で地形図を描いてみる． カラーパレットは，範囲が 0〜1 の copper を使った．

$ xyz2grd crater.txt -Gcrater.grd -R-5/5/-5/5 -I0.03125
$ makecpt -Ccopper -Z > c.cpt
$ grdimage crater.grd -Jx1 -Cc.cpt -Ba2 > crater.eps

クレーター

で，できあがった図だが…，リング状の地形があるのはわかるが， なんかもっさりしていて，断面図に見られるような三角っぽさがわからない．

そこで，この図に陰影をつけてみたいのだが， そのためには，標高データのほかに各グリッドの明暗のデータが必要となる． 各グリッドの明暗を決める方法であるが， 光源方向から見て上り坂で明るく，下り坂で暗くなるので， 光源方向から見た傾きを計算すればよい（下図）． 傾きが正のグリッドを明るく，負のグリッドを暗くするのである．

概念図

各グリッド面が，ある方角から見てどの程度傾いているかを計算するコマンドが，grdgradient である． grdgradient では，最低限，元データの grdファイル名， 計算結果の grdファイル名，どちらから傾きを計算するかの方角を指定する必要がある． 計算結果のファイル名は -Gオプションで指定し，方角は -Aオプションで指定する． 方角の値は，北（y軸正方向）が 0で，時計回りに度単位であたえる．

では，北から見た傾きを計算してみよう．

$ grdgradient crater.grd -Ggradient.grd -A0

計算結果 gradient.grd の分布は，以下のようになっている．

$ makecpt -Cpolar -Z > p.cpt
$ grdimage gradient.grd -Jx1 -Cp.cpt -Ba2 > gradient.eps

傾斜分布

赤が正，青が負である． 北半分の外側と，南半分の内側が正（上り坂）， 北半分の内側と，南半分の外側が負（下り坂）になっているのがわかる．

この傾斜の grdファイルを，grdimage に -Iオプションであたえると，陰影のある図ができあがる．

$ grdimage crater.grd -Jx1 -Cc.cpt -Igradient.grd -Ba2 > crater.eps

クレーター

山裾の立ち上がりや稜線がシャープで，斜面がツルツルであることがよくわかる図になった．

光源の位置を変えてみよう． -A90 で，東（X軸正方向）から光が当たった図になる．

$ grdgradient crater.grd -Ggradient.grd -A90
$ grdimage crater.grd -Jx1 -Cc.cpt -Igradient.grd -Ba2 > crater.eps

クレーター

-A180 で，光源が南になる．

$ grdgradient crater.grd -Ggradient.grd -A180
$ grdimage crater.grd -Jx1 -Cc.cpt -Igradient.grd -Ba2 > crater.eps

クレーター

-A270 で，光源が西になる．

$ grdgradient crater.grd -Ggradient.grd -A270
$ grdimage crater.grd -Jx1 -Cc.cpt -Igradient.grd -Ba2 > crater.eps

クレーター

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

(C) 2011 KDO