Main contents

☆主なコンテンツ
1、新着論文 2、論文概説 3、コラム 4、本のレビュー 5、雑記(PC・研究関連)
6、気になった一文集(日本語English) 7、日記(日本語English) おまけTwilog

2022年6月20日月曜日

ETOPO1の標高データを使い、色付きの等深線をGMTで描く

※GMTのバージョン4.5.8で正常に動くことを確認しています。最新のバージョンではダメのようです(2018.4.19追記)

※GMT6.3.0でも動くようにスクリプトを書き換えました。(2022.6.20追記)

#!/bin/sh

range=0/50/-50/-30
#size is in cm scale
size=12
reso=l+f
xanot=a10f5
yanot=a10f5
out=S_Atlantic
xy=S_Atlantic_location

gmt grdcut ETOPO1_Bed_c_gmt4.grd -R${range} -G${out}.grd
gmt makecpt -Cocean -T-8000/0/500 -Z > ${out}.cpt

gmt grdimage ${out}.grd -R${range} -JM${size} -C${out}.cpt -E100 -P -K > ${out}.eps


gmt grdcontour ${out}.grd -JM${size} -C2000 -W0.5 -L-6000/-200 -O -K >> ${out}.eps

gmt psscale -Ba2000g1000f1000 -C${out}.cpt -D6c/-1c/12c/0.3ch -O -K >> ${out}.eps
gmt pscoast -R${range} -JM${size} -D${reso} -W1 -Ggray -B${xanot}/${yanot} -A30 -O -K >> ${out}.eps

gmt psxy ${xy}.txt -JM -R -Sc6p -Gwhite -W0.5 -O -K >> ${out}.eps
gmt pstext ${xy}.txt -JM -R -F+f10p,Helvetica,white -D0.2/-0.3 -O -K >> ${out}.eps

open ${out}.eps

これによって以下のような図が描けます。



psxyとpstextに読み込ませているテキストファイルは以下です。

14.482 -39.637 PC10/MC14
44.970 -39.030 SK200/17
13.563 -40.863 ODP1088
9.893 -39.063 ODP1089
8.900 -41.087 ODP1090

〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜

前回の記事に引き続き、GMTを使っての海底地形図に等深線を描く方法のまとめです。
今回は海底に色を併せて塗ってみます。

以下のスクリプトを実行すると、以下の図が描けます。
今回使っているのはETOPO1_Ice_g_gmt4.grdのNOAAの無料データになります。ETOPO1_Bed_g_gmt4.grdでも大丈夫ですよ。




〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜

range=130/140/30/35
#size is in cm scale
size=12
reso=f
file=test
xy=location

grdcut ETOPO1_Ice_g_gmt4.grd -R${range} -G${file}.grd
makecpt -Cocean -T-6000/0/100 -Z > ${file}.cpt
#
grdimage ${file}.grd -R${range} -JM${size} -C${file}.cpt -E100 -P -K > ${file}.eps
#
grdcontour ${file}.grd -JM${size} -C1000 -W0.5 -L-5000/-200 -A1000tf8 -O -K >> ${file}.eps
pscoast -R${range} -JM${size} -D${reso} -W1 -A -G220 -O -K -Lf135/36/36/200k >> ${file}.eps
psscale -Ba2000g1000f1000 -C${file}.cpt -D6c/-1c/6c/0.3ch -O -K >> ${file}.eps
#psxy ${xy}.txt -JM -R -W -Sd6p -Gred -O -K >> ${file}.eps
#pstext ${xy}.txt -JM -R -O -K -Gred -D0.2/0 >> ${file}.eps
psbasemap -R${range} -JM${size} -Ba1f1/a1f1WSne -O >> ${file}.eps

〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜

1、grdcutについて
オリジナルのETOPO1_Ice_g_gmt4.grdは大変重いので、必要な箇所だけをカットします。

2、makecptについて
「-T」オプションで色のセットを選択できます。今回はoceanを使っていますが、他にも色々用意されています。
「-T」オプションで色の塗り方・間隔を指定します。z値の下限/上限/間隔の順です。

3、grdimageについて
実際に地図に色を塗っているコマンドラインです。
「-E」オプションで解像度(単位はdpi)を指定します。このオプションを付けないと大変荒い図になってしまいます。

4、grdcontourについて
「-C」オプションがコンターの間隔、
「-A」オプションがアノーテーションを付ける間隔(図中に書かれるコンターの数字のこと)数字のあとにtを付けないと、白い四角が背後に映り込んでしまいます。fのあとの数字はフォントのサイズです。
「-L」オプションがコンターを描く下限/上限になります。
「-W」オプションでコンターの太さを指定しています。

5、pscoastについて
「-D」オプションで海岸線のデータの解像度を指定します。上から「f: fine」「h: high」「i: intermediate」「l: low」「c: coarse」
「-W」オプションで沿岸線の太さを指定しています。

6、psscaleについて
色分けの説明に必要な、凡例を与えます。
「-B」オプションでラベルの間隔などを指定します。a数値g目盛りf格子の順です。
「-D」オプションで凡例の大きさを指定します。最後にhを付けると水平方向(horizontal)のスケールバー、抜かすと垂直方向に描かれます。中心の位置(左端から~cm)/上端の位置(下端から~cm)/スケールバーの長さ(~cm)/スケールバーの幅(cm)

7、psxyについて
先頭に#を付けて読み込まないようにしていますが、#を外してプロットしたい点をまとめたファイル(ここではlocation.txt)を用意しておくと、指定した緯度・経度にアイコンを打つことができます。
「-S」オプションでアイコンの形とサイズ、「-G」オプションでアイコンの色を指定します。
また、pstextで個々の点に名前をつけることも可能です。
以下のテキストファイルを読み込ませて描いた図を下の方に載せておきます。

130.650318 31.585311 8 0 1 LM Sakura-jima
134.189459 33.253410 8 0 1 CT Muroto-misaki
139.257435 34.381465 8 0 1 RT Nii-jima

左から順に、経度/緯度/テキストの大きさ/テキストのフォント(0はHelvetica)/テキストの角度/テキストの位置(アイコンの横・下・左などを指定)/書きたいテキスト
それぞれの間は半角スペースで区切ります。
「-D」オプションでテキストの位置を微調整します。x軸方向/x軸方向に移動させることができます。
より詳しい情報はこちらなどのページを参照のこと。

psxyは同じファイルの先頭2列のみを読み込んでいることになります。

8、psbasemapについて
縦軸、横軸を描きます。
「-B」オプションで目盛りやラベルの間隔を指定します。x軸のラベル・目盛りの間隔/y軸のラベル・目盛りの間隔。最後のWSneは西と南には指定したラベルを打つが、北と東はラベルなしの軸にするという意味になります。

※最後に
-Oと-Kがほとんどすべての行に出てきていますが、描画用のラインの最初には「-O」が不要、最後のラインには「-K」が不要なので、注意してください。

〜〜〜

以下は編集例。

1、makecptの-Cオプションを「-Cocean」→「-Cgray」に変更。白黒になります。


2、同じ図にポイントを打ってみます。上述の「location.txt」ファイルを読み込ませます。



3、もうちょっとクローズアップしてみる
上のスクリプトの
rangeを132/134/32/34に
makecptの「-T」を-2500/0/100に
grdcontourの「-C」を100に、「-L」を-800/-100に、「-A」を200tf8に、
pscoastの「-L」をf133/34.2/34.2/50kに変更。




間違いがあったらご指摘ください〜。