GMT

更新日2013年7月9日

GMT

以下では包括的な説明をするつもりはなく、GrADSの代わりにGMTを使うとどんなことができるか、を紹介するのが目的ですので、本当に勉強したいなら
homepage1.nifty.com/~kdo/gmt00_index.html(いちからはじめるGMT)
をみて勉強するほうが得策です。
中級者向けとして、これは、GMTのコマンドをmanすると出てくる文章をただ日本語訳したまとめサイトですが、
http://gmt.shin-gen.jp/GMT4.1.4/gmt_man.html(GMT Online Man Page)
も紹介しておきます。

GMTを用いる信念：

論文のような地図を書く場合、投影法の選択が重要である
解析ツールを使う代わりに、手で計算することを経験することが重要である
どんな便利な解析ツールも最初の一人は手作業を経験している

インストール

インストール自体は略。
というのも、ubuntuだと何も考えなくても

$ sudo apt-get install gmt

でインストールはできてしまうので。各自調べてみてください。

パスの通し方は、 .bashrcでは、

export PATH=$PATH:/usr/local/bin（←もし最初からあるなら書き加える必要ない）
export GMTHOME=/usr/lib/gmt
export PATH=$GMTHOME/bin:$PATH
export MANPATH=$GMTHOME/man:$MANPATH

とします。

グラフを書く

まず、データを用意します。ファイル名は"1.dat"とします。

これを、psxyをつかって、

psxy 1.dat -Jx1 -R-0.5/6.5/6/14 -K -Sc0.5 -Gblue>tmp.ps
psxy 1.dat -Jx1 -R-0.5/6.5/6/14 -K -W5,blue -O>>tmp.ps

と、丸と線を重ね書きします。

コメントを入れる

echo 0 9 8 0 0 "MC" aaa | psxy -J -R -K -O -N>>tmp.ps
（左から、緯度、経度、文字の大きさ、0、0、位置、文字列）

と、上のグラフの[0,9]の点に、文字さイズ8で"aaa"と文字がかかれます。なお日本語はこのやり方ではできない。
グラフでなく地図に入れる場合も考え方は同じで、緯度経度を指定します。

地図を書く

これがよくある使い方。地図に情報を重ね合わせます。アメダスによる2008-12年の5年間の猛暑日日数。

以下のようなデータを用意し、

#経度　#緯度　#地点名　#観測値
140.77 36.84 北茨城 0
140.65 36.60 日立 2
140.35 36.78 大子 69
140.33 36.61 常陸大宮 17
140.47 36.38 水戸 17
…

スクリプト

#!/bin/sh
basecpt='basemap.cpt'
file='moushobi.dat'
out='o.eps'
region='138/141/34.5/37.5'
proj='M16'
grid='0.01'

awk '{print $1,$2,$4}' $file>tmp
awk '{print $1,$2,14,0,0,"MC",$4}' $file>tmp2
awk '{print $1,$2,12,0,0,"MC",$3}' $file>tmp3
surface tmp -Gtest.grd -I$grid -R$region -T0.25
grdimage japan.grd -J$proj -R$region -C$basecpt -K > $out
grdcontour test.grd -J$proj -C20 -A20 -K -O>> $out
pscoast -R$region -J$proj -W -Df -S -Bg0.5a0.5 -K -O>> $out
psxy -R$region -J$proj tmp -Sc0.7 -G255/255/128 -K -O >>$out
pstext tmp2 -R$region -J$proj -K -O >> $out

で描画しています。こういうときにawkが大変便利です。

例をもう1つ。

Rignot et al.,2013,Ice-Shelf Melting Around Antarcticaで推定された棚氷底面融解量のマップ。単位はGt/yr.
赤色で示した棚氷と、○印が重ならないようにわざわざ手作業でずらしています。本当はもう少し文字を大きくしたいですが。

投影法について

メルカトル図法

pscoast -Jm1:200000000 -R0/360/-80/80 -K -W -Bg15a15/g15a15> test.eps

もっとも馴染みのある図法だと思います。
普通の書き方ですが、GrADSのlatlonと違うのは高緯度は引きのばされることです。全球図には向いていない。
pscoastで-Jxととることは許されません。かわりに-Jqを使います。

等距円筒図法

pscoast -Jq150/1:150000000 -R0/360/-90/90 -K -W -Bg30a30/g15a15> test.eps

GrADSのlatlonに相当。

pscoast -Jq240/-83/1:29000000

などとすれば縦に引き伸ばせます。

エケルト第4図法

pscoast -Jkf150/1:150000000 -R0/360/-90/90 -K -W -Bg30a60/g15a15> test.eps

エケルト第4図法。全球の分布図をかくのに適している（だろう）
気象庁の世界の平均気温（陸上における地上付近の気温と海面水温）の偏差の分布図→（http://www.data.kishou.go.jp/climate/cpdinfo/temp/map/temp_map.html）でも使われています。
（いかにもそうみえますが）これは海面水温ではありません。なんでしょう

極投影図法polar projection

pscoast -JS0/-90/6.5i -R0/360/-90/-60 -K -W -Bg15a15/g5a10> test.eps

極投影polar projection。GrADSのnps,spsに相当。北極周辺、南極周辺図に適しています。

pscoast -JS135/35/16 -R0/360/45/90 -K -W -Bg15a15/g5a10> test.eps

極以外を中心にとることもできなくもない。

ミラー図法

pscoast -Jj150/1:150000000 -R0/360/-90/90 -K -W -Bg30a30/g15a15> test.eps

ミラー図法。latlonとメルカトルの中間？

カッシーニ図法

pscoast -Jc135/35/1:40000000 -R90/15/210/50r -K -W -Bg5a10/g5a10> test.eps

-Rの最後に-rを入れることで正方形に切り取れます。-R左下隅の点の経度/左下隅の点の緯度/右上隅の点の経度/右上隅の点の緯度rです。

http://topex.ucsd.edu/marine_topo/mar_topo.html
の地形データ。AMOCの沈み込んだ水がアイスランドridge（？）の谷間を通って北大西洋に排出される。
もうすこし視点を離すと

こんな地形をしている。

ユニバーサル横メルカトル図法(UTM図法)

pscoast -Ju53/1:4000000 -R132/32/139/38r -K -G192/192/192 -X2 -Bg1a2/g1a2 -Df>$out
echo 135.518 34.682 16 0 0 "MC" Osaka | pstext -J -R -K -O >>$out
echo 133.917 34.660 16 0 0 "MC" Okayama | pstext -J -R -K -O >>$out
echo 136.965 35.167 16 0 0 "MC" Nagoya | pstext -J -R -K -O >>$out

1:25000地形図で使われる投影方法です。近畿はu53、関東はu54、九州はu52を使います。

等値線を書く

あまりにも単純ですが等値線。縦軸が水深(dbar)、横軸が塩分濃度で、コンターが海水の凝固点です。

色塗りをする

ベクトルを描画する

ベクトルは、角度（北を0度として時計回りの角度、0から360まで）と長さを指定します。

135 35 45 1

なら135E 35Nで北東からの風（南西へ向かう風）の矢印がかかれます。
ほかに、風向東西成分と南北成分の2つの.grdファイルを用意してgrdvectorするやりかたもあります。

grdvector u.grd v.grd

JMやJQなど、角が正しい図法以外の地図でgrdvectorをする場合、-Tをつけないとベクトルの向きが正しくならない。

3次元遠近法を用いた立体図を描画する

BEDMAPより。grdviewで立体図を描画できます。水平方向の角度と仰角を指定します。

海洋モデルの地形を３次元的に表現して見ようとおもったがあまりうまくはいっていない。わかりづらいが南極周辺を北から見下ろしている。

側線に沿った断面の形状を描画する

これは南極の断面図です。緑が氷床もしくは棚氷、茶色が地面、青色が海洋、赤線が海水準を示しています。データはBEDMAP2から。grdファイルを用意すれば、fortranプログラムをまわすことなく、スクリプトで始点と終点を指定することでこのような断面図を一発で書くことができます。とはいいながらこれはgrdtrackで書いたわけではなく頑張ってfortranとawkを組み合わせて書いています。

GMTでアニメーション

サンプルスクリプトでは、iで時刻管理しています。

#!/bin/sh
out=test.eps
region=225/-45/45/-45r
line=g30/g15
cpt=cpt/ai.cpt
proj=S0/-90/16

i=10
while [ $i -lt 71 ]
do
rm -rf $out
nt=`expr $nt + 253941 `
echo $nt
head -n $nt tmp | tail -n 253941 > tmp2
echo 0 -90 16 16 | psxy -J$proj -R$region -Sr -Ggray -K >$out
psxy -J$proj -R$region tmp2 -Sc0.06 -C$cpt -K -O -B$line>> $out
grdcontour -J -R shelf.grd -C0.5 -L0.3/0.7 -K -O -W2,black >>$out
psscale -D8/-0.5/16/0.5h -C$cpt -K -O -N -E -L>>$out
j=`expr $i - 9`
echo 0 -90 24 0 0 "MC" $j | pstext -J$proj -R$region -Gblack -K -O >>$out
i=`expr $i + 1`
convert $out img$i.gif
convert -rotate 90 img$i.gif img$i.gif
done

253941というのは1時刻の行数です。最後にconvert -delay 100 -loop 10 img*.gif anime.gifなどとします。
exprの前後のは`←バッククオーテーションです。'ではありません。また半角スペースもこの通りに入れないとエラーになります。

ヒストグラムの描画

pshistogram tmp -W100 -L -BSWne -Ggray -A -Jx1/-0.008 -Z1>test.eps

データは

のように1列に値がならんでいるものとします。

以下は、なれてきた人向け

gridにどんな数値が入っているかを図示

全gridでかくととても見ていられないので、数gridおきに書いています。
ちょっとしたスクリプトで、どんな値が入っているかを視覚的に確認できます。

マスクアウトなどの演算をする。

以下、ただの趣味

大西洋子午面循環の描画

陸地はマスクアウト。背景色はないほうがよいか。
まず、GTOOL形式のvo（海洋南北流速、すでに地理座標系に直したもの）があるとします。
大西洋の地形ファイルATLを用意して、fortranプログラムで各点(lat,depthで定義)の東西積分流量を求め、低層から足し合わせることで流線函数を求めます。

back