ホーム>資料>パソコン

パソコンのプチ情報

前の記事   次の記事

b15.GIF GMT講座 その9:自作データでコンター図

さていよいよ、自作データでコンター図を描いてみましょう。
そのためには、GRDファイルを自作する必要があります。でも用意するデータは
ある場所(X,Y)でのある値(F)がひたすら並んでいるアスキーファイルです。
これをGRDファイルの形式に変換すればよいのです。

1)まずデータを示します。これはある場所での電気探査の結果です。
第1列目が観測点位置、第2列目が電極間隔(=およその探査深度)、第5列目が
得られた見掛比抵抗(探査深度までの地下の平均的な電気の通りにくさ)です。
これらを元にして地下を可視化してみましょう。
※なおウエンナー法電気探査です。
(ファイル名:ex1.dat)
280 3 0.365 1.19 57.8158648
280 6 0.206 1.15 67.53058295
290 3 0.365 1.19 57.8158648
290 6 0.184 1.18 58.78505576
290 9 0.139 1.15 68.35012886
290 13 0.122 1.11 89.77596304
300 3 0.357 1.19 56.54866776
300 6 0.178 1.16 57.84863714
300 9 0.136 1.15 66.87494623
300 13 0.124 1.13 89.63269659
300 18 0.122 1.11 124.3051796
310 3 0.364 1.2 57.1769863
310 6 0.175 1.16 56.87366011
310 9 0.14 1.15 68.84185641
310 13 0.12 1.12 87.51579535
320 3 0.381 1.19 60.35025887
320 6 0.227 1.14 75.06752972

2)これをコンター図にするGMTスクリプトを下記に示します。
(スクリプト名:contour.gmt)
------------------------------------------------------------
#!/bin/tcsh -f
# GMT script for making bathymetry map (Japan)

gmtdefaults -D >! .gmtdefaults
gmtset LABEL_FONT_SIZE 12 HEADER_FONT_SIZE 18

awk '{print $1-300, $2, $5}' < ex1.dat > tmp.dat
blockmean tmp.dat -R-30/30/0/20 -I2/1 | \
surface -R -I2/1 -Gtmp.grd
makecpt -Crainbow -Z -I -T40/150/5 >! tmp.cpt

grdimage tmp.grd -JX6.6/-3 -R-22/22/0/20 -Ctmp.cpt \
-Ba5/a5WsNe:."ex1": -Xa1.0i -Ya2.0i -P -K >! tmp.ps

grdcontour tmp.grd -JX -R -L30/150 -C10 -A20\
-Xa1.0i -Ya2.0i -P -K -O >> tmp.ps

psscale -Ctmp.cpt -D3.5/1.3/3.0/0.2h \
-B50/:"App.Resistivity(Ohm-m)": -P -K -O >> tmp.ps

# Plot observed location
awk '{print $1, $2}' < tmp.dat | \
psxy -JX -R -Sc0.1i -G0 -Xa1.0i -Ya2.0i -P -K -O >> tmp.ps

echo 0 0 18 0 1 1 VES_Survey | \
pstext -JX4.5/1 -R0/16/0/1 -Xa3.5i -Ya6.0i -P -O >> tmp.ps

gv tmp.ps
mv tmp.ps ex1.ps
rm tmp*.dat tmp*.grd tmp*.cpt *grd.shade
rm .gmt*
------------------------------------------------------------

新たに出てきた部分のみ、下記に解説します。
・blockmean: アスキーデータ(tmp.dat)に基づいて、観測データ間の平均値を
 作り出していきます。-Rで範囲指定(横=-30~30、深さ=0~20)、-Iで作り出す間隔
 (横2ずつ、深さ1ずつ)。ここでIの値は、Rの値を割り切れるように設定してください。
・surface: アスキーデータに従いながら、GRDファイルを作成します。
 -R、-Iはblockmeanと合わせておくとよいようです。-Gで出力ファイル名を指定します。
・makecptの-Zはスムーズなカラーパレット作成を意味します。-Iは数字の大小に対する
 色を反転します。
・grdimage: WSne → WsNe として、Y軸ラベルなどを図の上側に表示しました。
・awk '{print $1, $2}' < tmp.dat | \
 psxy -JX -R -Sc0.1i -G0 -Xa1.0i -Ya2.0i -P -K -O >> tmp.ps
 この行では、観測データの位置(横・縦)を円で図示しています。
 -Sc0.1iは0.1の大きさの円を描いています。-G0は黒く塗りつぶしを意味します。
・echo 0 0 18 0 1 1 VES_Survey | \
 pstext -JX4.5/1 -R0/16/0/1 -Xa3.5i -Ya6.0i -P -O >> tmp.ps
 ここではシンプルに文字だけを表示しています。複数のグラフの共通タイトルなどに
 使うとよいでしょう。

意外に簡単ですね。blockmeanやsurfaceの-Iの値を変えるとコンター図の顔つきも
変わってきます。試してみてください。

続く。
次回は未定ですが、7月中ごろかな?

[関連画像]


図1

前の記事

目次に戻る

次の記事