2011年 9月 1日 掲載
陰影のある地形図を描くつづきである． 前回あつかった地形データは，xyz の値が「km」単位だったが， ふつう地図といったら，xy が「度」，z が「m」だろう． 今回は，このふつうの単位のデータを使って地図を描いてみる．
まずは，もとの地形データの用意である． たいていの GMT の解説記事では， ここで公開されている地形データの紹介があって， ダウンロードしたデータの加工の仕方を説明するところだが， この連載では AWK で作ってしまう． 今回使う地形データを作る AWKスクリプトが，以下の mountains.awk である．
BEGIN{
for(lon = 138.4; lon <= 139.2; lon += 0.002){
for(lat = 34.5; lat <= 35.6; lat += 0.002){
h = height(lon, lat, 138.73, 35.36, 3776)
h = max(h, height(lon, lat, 139.02, 35.23, 1438))
h = max(h, height(lon, lat, 139.00, 34.86, 1405))
print lon, lat, h
}
}
}
function height(lon, lat, lon0, lat0, h0){
dx = (lon - lon0) * cos(lat0 * 3.1415 / 180) * 1852 * 60
dy = (lat - lat0) * 1852 * 60
r = sqrt(dx ^ 2 + dy ^ 2)
z = h0 - r * 0.2
if(z < 0) z = 0
return z
}
function max(a, b){
if(a > b){
return a
}else{
return b
}
}
このスクリプトを実行すると， 富士山，箱根山，天城山を円錐で表した地形データができあがる．
$ awk -f mountains.awk > mountains.txt $ xyz2grd mountains.txt -Gmountains.grd -R138.5/139.2/34.5/35.5 -I0.002
このデータには海がないので， 描画の際には，pscoast -S でマスクをかけなくてはならない．
$ makecpt -Cglobe > g.cpt $ grdimage mountains.grd -Cg.cpt -Jm10 -K > mountains.eps $ pscoast -J -R -Sblue -Df -O >> mountains.eps
で，できた図がこれ↓である．
しょぼっ． まあ，GMTコマンドの解説に使うだけだから…．
さて，さっそくこの図に陰影をつけてみる． 前回と同じようにしてやってみよう．
$ grdgradient mountains.grd -Gintensity.grd -A90 $ grdimage mountains.grd -Cg.cpt -Jm10 -Iintensity.grd -K > mountains.eps $ pscoast -J -R -Sblue -Df -O >> mountains.eps
これはひどい．日向は真っ白，日陰は真っ黒けである． 陰影強調のためのファイル intensity.grd がどうなっているか， grdinfoコマンドでみてみよう．
$ grdinfo intensity.grd intensity.grd: Title: Directional derivative(s) intensity.grd: Command: grdgradient mountains.grd -Gintensity.grd -A90 intensity.grd: Remark: intensity.grd: Gridline node registration used intensity.grd: Grid file format: nf (# 18) GMT netCDF format (float) (COARDS-compliant) [DEFAULT] intensity.grd: x_min: 138.5 x_max: 139.2 x_inc: 0.002 name: x nx: 351 intensity.grd: y_min: 34.5 y_max: 35.5 y_inc: 0.002 name: y ny: 501 intensity.grd: z_min: -18237.5175781 z_max: 18237.5175781 name: z intensity.grd: scale_factor: 1 add_offset: 0
z がとんでもない値になっている． こんなことになるのは，「m」単位の標高差を「度」単位の水平距離で割っているからである． 傾斜をちゃんと計算するには，水平座標と標高の単位をあわせなくてはならない． grdgradientコマンドには，この計算をやってくれる -Mオプションが用意されている． -Mオプションをつけると，度単位の xy座標から，m単位の距離を計算して，傾斜の計算に使ってくれるのである． では，さっそく．
$ grdgradient mountains.grd -Gintensity.grd -A90 -M $ grdimage mountains.grd -Cg.cpt -Jm10 -Iintensity.grd -K > mountains.eps $ pscoast -J -R -Sblue -Df -O >> mountains.eps
おお，今度はまともな図になった． けど，ちょっと陰影がうすいかな． また，intensitiy.grdファイルを grdinfo で見てみよう．
$ grdinfo intensity.grd intensity.grd: Title: Directional derivative(s) intensity.grd: Command: grdgradient mountains.grd -Gintensity.grd -A90 -M intensity.grd: Remark: intensity.grd: Gridline node registration used intensity.grd: Grid file format: nf (# 18) GMT netCDF format (float) (COARDS-compliant) [DEFAULT] intensity.grd: x_min: 138.5 x_max: 139.2 x_inc: 0.002 name: x nx: 351 intensity.grd: y_min: 34.5 y_max: 35.5 y_inc: 0.002 name: y ny: 501 intensity.grd: z_min: -0.199882388115 z_max: 0.199882388115 name: z intensity.grd: scale_factor: 1 add_offset: 0
z の大きさが ±0.2 くらいしかない． そりゃそうだ． この地形データは，傾き 0.2 の円錐形の山で作ってあるんだから，しょうがない．
grdgradientコマンドには，傾斜の値を正規化する -Nオプションがある． -Nオプションには，いろいろパラメータを指定できるのだが， とりあえず，今回は無指定で使ってみよう．
$ grdgradient mountains.grd -Gintensity.grd -A90 -M -N $ grdinfo intensity.grd intensity.grd: Title: Normalized directional derivative(s) intensity.grd: Command: grdgradient mountains.grd -Gintensity.grd -A90 -M -N intensity.grd: Remark: intensity.grd: Gridline node registration used intensity.grd: Grid file format: nf (# 18) GMT netCDF format (float) (COARDS-compliant) [DEFAULT] intensity.grd: x_min: 138.5 x_max: 139.2 x_inc: 0.002 name: x nx: 351 intensity.grd: y_min: 34.5 y_max: 35.5 y_inc: 0.002 name: y ny: 501 intensity.grd: z_min: -1 z_max: 1 name: z intensity.grd: scale_factor: 1 add_offset: 0
z の値が -1 から 1 の範囲に正規化された． このファイルを使うと，地形図もだいぶよくなる．
$ grdimage mountains.grd -Cg.cpt -Jm10 -Iintensity.grd -K > mountains.eps $ pscoast -J -R -Sblue -Df -O >> mountains.eps
-N にあたえるパラメータと作用については，次回説明することにしたいが， とりあえず，上でやった無指定の場合にどういうことをしているかくらいは， 今回説明しておかなきゃならないかな？
-Nオプションをパラメータ無指定で使うと， もとの傾斜の値 g は，以下の式で正規化される．
平均値との差を，平均値との差の大きさの最大値で割っているので， -1 から 1 の範囲の値になる，というわけだね．
[ KDOホームにもどる] [ いちからはじめる GMT インデックスにもどる ]
(C) 2011 KDO