ホーム＞資料＞パソコン

パソコンのプチ情報

前の記事   次の記事

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

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

１）まずデータを示します。これはある場所での電気探査の結果です。
第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

２）これをコンター図にする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

前の記事

目次に戻る

 次の記事