2008-5-13:標高データを用いてGMTで絵を描く - BioEoS::環境解析学研究室::

前回の授業では，標高データの説明をしました．
今回の授業では，津市西部の2次メッシュの標高データをつかってGMTで絵を描きたいと思います．

1．GMTコマンド「GRDIMAGE」の説明

オンラインのGMTのコマンドリファレンス（GMTのホームページ内⇒DOCS⇒GMT Manual Pages⇒HTML Format）以外にもGMTのフォルダ内にあるwwwというフォルダにも同じものがあります．

今回使うコマンドは「GRDIMAGE」です．
名前の通り，グリッドファイルを画像として表示するコマンドです．

リファレンスの「grdimage」の欄を見ると，次のように書いてあります．

grdimage grd_z | grd_r grd_g grd_b −Ccptfile −Jparameters [ −B[p|s]parameters ] [ −Ei|dpi ] [ −G[f|b]color ] [ −Iintensfile] [ −K ] [ −M ] [ −N ] [ −O ] [ −P ] [ −Q ][−Rwest/east/south/north[r] ] [ −S[-]b|c|l|n[/threshold] ] [ −T ] [ −U[just/dx/dy/][c|label] ] [ −V ] [ −X[a|c|r][x-shift[u]] ] [ −Y[a|c|r][y-shift[u]] ] [ −ccopies ] [ −f[i|o]colinfo ]

マイナスの後に大文字の英字で表されているもの（-Cなど）をオプションと言います．
ここで，括弧[]の付いていないものは，絶対に必要なオプションになります．
また，縦棒(|)は「または」という意味です．スラッシュ（/）と間違えないように気をつけてください．

このコマンドを使うにはnetCDF形式のグリッドデータが必要です．

そこで，グリッドデータを作成するGMTコマンド「XYZ2GRD」を使って，グリッドデータを作成し，描画させることにします．

2．XYZデータを作成する

通常のXYZデータは3列のデータ（それぞれの列にXの値，Yの値，Zの値が一列に記述されている）ですが，今回は非常に規則的なデータなので，Zの値のみを横200×縦200に記述するデータを使用します．

使用するデータは「523603.MEM」 (O:\DATA\Japan\5236\523603.MEM)です．
1行目はヘッダー（データの説明など）なので，2行目以降200行を読み込み，並び替えします．
エディターを起動させて,「EoS」にある「ruby」のフォルダに「mem2z.rb」という名前でRubyスクリプトを作成します．
Rubyのスクリプトは以下の通りです．

file = "O:/DATA/Japan2/5236/523603.MEM"
open(file) do |f|
　f.gets
　200.times do
　　line = f.gets
　　200.times do |i|
　　　p line[9+i＊5,4].to_i
　　end
　end
end

スクリプトの説明をします．以下の番号が行番号と思ってください．

1. "O:/DATA/Japan2/5236/523603.MEM"をfileに代入します．
2. fileを開いたあとをfと置く．
3. 開いたfileから1行読み取る．
4. 3行目に行ったことを200回繰り返す．
5. f.getsをlineに代入する．
6. 200回繰り返すことをiと置く．
7. fileから1行取ったものの頭から9+i＊5の数字から4つの数字を取る．

この作成したRubyファイルをコマンドプロンプトを立ち上げて作業フォルダに移動して実行してください．

ruby mem2z.rb

すると，数値が縦に40,000個出てきます．
数が多いので結果は省略します．

この実行された結果の数値をファイルとして保存します．

ruby mem2z.rb > 523603.z

すると，作業フォルダに保存されます．

保存されたファイルを見るコマンドは以下の通りです．

type 523603.z | more

3．グリッドデータを作成する

GMTコマンド「XYZ2GRD」を使います．

xyz2grd 523603.z -G523603.grd -I2.25c/1.5c -R136:22:30/136:30/34:40/34:45 -F -Z -V

ここで，「-G」は作成するグリッドデータの指定を，「-I」はデータ間隔の指定(今回の場合はx軸方向に2.25秒,y軸方向に1.5秒という意味)を表わしています.
すると523603.grdができます．
グリッドファイルの内容確認は

grdinfo 523603.grd

でできます

4．gmtdefaultsの変更

GMTの設定を変更するために新しくデフォルトのファイルを作成します．

gmtdefaults -L > .gmtdefaults4

新しくできたgmtdefaults4というファイルをエディターで開いて内容を書き換えます．
書き換える箇所は3か所あります．

• 7行目：PAPER_MEDIA=noto

を以下のように変更します．

PAPER_MEDIA=a4+

• 27行目：PLOT_DEGREE_FORMAT=+ddd:mm:ss

を以下のように変更します．

PLOT_DEGREE_FORMAT=ddd:mm:ssF

• 89行目：MEASURE_UNIT=inch

を以下のように変更します．

MEASURE_UNIT=cm

変更点は次の3点です．

1. [noto]を[a4+]に変更．
2. [+ddd:mm:ss]の[+]を取り，最後に[F]を加える．
3. [inch]を[cm]に変更．

以上です．
変更が終了したら，上書き保存してください．

5．カラーパレット（CPT）ファイルを作る

描画する際のカラーパレットを作ります．

grd2cpt 523603.grd -Ctopo > 523603.cpt

6．グリッドを描画する

グリッドファイルをポストスクリプトとして出します．

grdimage 523603.grd -C523603.cpt -JM16 -R136:22:30/136:30/34:40/34:45 -Ba5mf1m -P -V > 523603.ps

出来上がったポストスクリプトファイルをGSviewというアプリで表示した結果は以下の通りです．
クリックすれば別ウインドウに表示されます．

上に表示された画像に陰影をつける（北西方向からライトを当てる）には

grdgradient 523603.grd -G523603_i.grd -A0/270 -Ne0.6

で新しいグリッドファイル(G523603_i.grd )を作成します．
そして，上で作成したファイルを反映させてポストスクリプトファイルを作成します．

grdimage 523603.grd -C523603.cpt -JM16 -R136:22:30/136:30/34:40/34:45 -Ba5mf1m -P -V -I523603_i.grd > 523603_i.ps

出来上がったポストスクリプトファイルをGSviewというアプリで表示した結果は以下の通りです．
クリックすれば別ウインドウに表示されます．

GMTのコマンドは各自で調べてください．
不親切ですいません．
以上です．

---