National Institute of Advanced Industrial Science and Technology (AIST) This page is a page of the former research institute. We stopped updating on March 31.2001.
E-mail to webmaster (Japanese) E-mail to webmaster (English)

GMT(ver.3)の特徴と使用法

1996/04/11
地質調査所海洋地質部
棚橋 学

この文書の目的

地質調査所海洋地質部周辺のGMT,UNIXの完全な初心者が,非常に強力な地質地球物理データ処理系であるGMTを使いはじめるようになるための手がかりとなる基礎的な知識,ノウハウを示すことです.地質調査所以外の方はローカルな話題のところは適当に無視して下さい.ここでは主に海洋地球物理データのプロットの例を示していますが,GMTは一般的なランダムデータのグリッド化,可視化などにも非常に適したツールです.

GMTの特徴

GMTは(Gravity, Magnetism, Topography)データをはじめとするxyzデータを保存し,取り出し,加工し,表示し,印刷するというシステムです.現在ではGeneric Mapping Tools (GMT)と呼ばれています.その概要は製作者たちによりWessel and Smith (EOS, v.72, p.441-446, 1991)に紹介されています.GMTが開発されたいきさつなどはAppendixに簡単に書いておきました.
海洋地質部のgsjsantsにはv.3が,gsjsantにはv.3とv.2.1.4が,gsjgmtにはv.2.1が,リモートステーションのgsjrstnにはv.2.1が,白嶺丸のsun10とsparcltにはv.2.1がインストールされています.v.2ではグリッドファイルの形式がv.1とは変更され,機種依存性の低いNetCDF形式になりました.新しく使う場合は問題ないでしょうがv.1からあるコマンドの一部は使用法が変更されています.
GMTはもともとSUNワークステーションで開発されたので,SUNの上で使いやすいようになっていますが,現在ではHP,IBM/RS,DECstation,NeXT,Silicon Graphics等のUNIXワークステーションに移植され使われています.私自身もMacintosh上のUNIXシステムであるMachTenの上でインストールしてみました.速度の点で問題はありますがPowerBook Duo250の上でGMTでPostScriptファイルを作成し,GhostScriptで結果を表示することができました.
GMTにはいろいろな特徴がありますが,後で述べるCookbookによれば

  使いやすい,
  ファイル構造が単純である,
  高品位な出力が得られる,
  ビットマップドグラフィックスが得られる,
  カラー出力が可能,
  ユーザープログラムとのインタフェースが可能,
  UNIXとのインタフェースが良好,

といった長所があり,それに対して商業製品とは違って完全にはテストされていないためバグの存在する可能性があるという短所があります.しかし,v.1,v.2.1に比べると現在のバージョンはバグが少なく非常に使いやすいものになってきました.
GMTは図面の出力をほとんどすべてPostScript形式で行います.そのため,何枚分もの出力を重ね合わせる,エディタを用いて一部の編集をする,余分に何かを加えると言うことが可能であり,SUN Openwindowsを使用する場合には画面の上にプロットして結果をすぐに見ることができます.ネットワークを利用して他から画面出力を利用することは現在はできません.ポストスクリプトファイルをFTPやrcpで転送してから表示,印刷することになります.SUNのSPARCprinterや最近の高性能なPostScriptプリンタでは400 dpi-600 dpi程度の高品質なプリントが可能です.マッキントッシュのLaserWriterもPostScriptプリンタなのでDropPSのようなプログラムでプリント可能です.PostScriptの欠点として,ファイルが大きくなります.マッキントッシュでローカルトーク接続されているプリンタで出力しする場合には非常に時間がかかることがあるので注意が必要です.SONY TextronixのPhaser II,3MRainbowなどのカラーポストスクリプトプリンタでカラーの出力をすることも可能です.GMTv.2の機能でカラーのビットマップであるSunraster形式のファイルへの変換ができるので,さらに適当なファイル変換を行うことにより,さらに画像処理を適用するためにレタッチソフトを利用したり,QuickDrawなどの普通の低価格のカラープリンタにでも出力できるでしょう.また,コンターマップのような線画出力の場合には,UNIX上のGhostScriptを利用したaimakerや,マッキントッシュ上のEPSConverterによって,GMTで作成したポストスクリプトファイルをAdobe Illustrator88やv.3形式に変換することができるので,マッキントッシュ上で線画の編集や文字を追加したりすることも可能です.
GMTに関してはオンラインマニュアル用のファイルgmt.manとThe GMT-SYSTEM Version 3 Technical Reference and Cookbook by Paul Wessel and Walter H.F. Smith, 1993の2種類のドキュメントがあります.gmt.manには各コマンドの説明とGMTSYSTEMという概説の項目があります.

GMTのコマンド一覧

GMTを構成しているコマンドの説明を簡単にまとめておきます.(gmt.man中のGMTSYSTEMより)
   blockmean L2ノルムを用いた (x,y,z) データフィルタ
   blockmedian L1ノルムを用いた (x,y,z)データフィルタ
    filter1d 1次元(時系列)データセットにフィルタを適用する
   fitcircle 点のデータセットに対する最適な大円を求める
   gmtdefaults 現在のデフォールトの設定条件をリストする
   gmtset 現在の.gmtdefaultsファイルのパラメータの一部を変更する
   grdfilter 2次元データセットに対して空間領域でのフィルタを適用する
   grd2cpt グリッドファイルからGMTのカラーパレットファイルを作成する
   grd2xyz 2次元グリッドファイルからアスキーファイルへ変換する
   grdclip グリッドデータセットのzの値のレンジを制限する
   grdcontour 2次元グリッドデータセットのコンタリング
   grdcut 2次元グリッドファイルから一部の領域を切り出す
   grdedit 2次元グリッドファイルのヘッダ情報を書き換える
   grdfft 周波数領域での2次元グリッドファイルの操作
   grdgradient 2次元グリッドファイルからの方向勾配の計算
   grdhisteq 2次元グリッドファイルのヒストグラム平滑化
   grdimage 2次元グリッドファイルからのイメージの作成
   grdinfo 2次元グリッドファイルに関する情報の抽出
   grdmask 2次元グリッドファイルの特定の領域の外部の値を定数で置き換える
   grdmath 2次元グリッドファイルに数学的操作を加える
   grdpaste 複数の2次元グリッドファイルを縁に沿って連結する
   grdproject 2次元グリッドデータを新しい座標系に投影する
   grdreformat グリッドファイルを他のグリッドフォーマットに変換する
   grdsample 2次元グリッドデータセットから新しいグリッドをリサンプルする
   grdtrend 2次元グリッドファイルに対する多項式およびフーリエトレンドのあてはめ
   grdtrack 2次元グリッドファイルから1次元の測線に沿ってデータをサンプルする
   grdvector ベクトル場のプロット
   grdview 2次元グリッドデータセットの3次元鳥瞰イメージの作成
   mapproject 座標系の変換
   minmax アスキーファイル内の極端なデータを求める
   nearneighbor 隣接グリッド化手法
   project データの測線/大円への投影
   psbasemap ベースマップの作成
   psclip 多角形ファイルを用いてクリップする
   pscoast 海岸線をプロットし大陸を塗りつぶす
   Pscontour 三角網法による原データのコンタリングおよび画像化
   pshistogram ヒストグラムをプロットする
   psimage ラスターファイルをマップにプロットする
   psmask マップ上の領域をマスクするオーバーレイを作成する
   psmegaplot PostScriptファイルからポースターサイズのプロットを作成する
   psrose セクター/ローズダイアグラムをプロットする
   psscale マップ上にグレイスケールまたはカラースケールをプロットする
   pstext マップ上に文字列をプロットする
   pswiggle マップ上の測線に沿って異常を描く
   psxy マップ上にシンボル,多角形,線をプロットする
   psxyz 3次元イメージ上にシンボル,多角形,線をプロットする
   sample1d 1次元データセットのリサンプル
   spectrum1d 時系列から様々なスペクトル推定を行う
   splitxyz xyzファイルをいくつかのセグメントに分割する
   surface 連続曲率グリッド化アルゴリズム
   trend1d y = f(x)系列に多項式トレンドをあてはめる
   trend2d z = f(x,y)系列に多項式トレンドをあてはめる
   triangulation 最適なDelauneyの三角網化とグリッド化を実行する
   xyz2grd アスキーxyzファイルから2次元グリッドファイルへ変換する

以上54のコマンドがあります(v. 3)が,このほかにオプショナルな.gmtファイルという航海ファイルの操作に関係したコマンドが11個,その他にも数個のコマンドが利用できます.航海ファイルを扱うコマンドがオプショナルになったのは機種依存のバイナリファイルであるためです.将来的には新しいNetCDFに対応した形式の航海ファイルを扱うコマンドが提供される予定のようです.SUNのみで使うことには問題はないので調査所でインストールしてあるGMTv.2にはgmtファイル関係のコマンドも入れてあります.海洋地質地球物理の分野でGMTを使うときにはgmtファイルを利用することが中心になると思います.

これらのコマンドは
   dat2gmt アスキーファイルをgmtに変換する
   gmt2dat gmtファイルをアスキーファイルに変換する
   gmtinfo gmtファイルの航海の情報を与える
   gmtlegs 与えられた領域内の航海名を与える
   gmtlist gmtファイルからあるデータを抽出する
   gmtpath gmtファイルのフルパス名を与える
   gmttrack 航跡をプロットする
   mgd77togmt MGD77形式の航海ファイルからgmtファイルを作成する
   binlegs gmtデータベース作成プログラム
   gmt2bin gmtファイルをデータベース化するためのプログラム
   gmtedit Sunviewで動くgmtファイルのヴィジュアルエディタ

このほか,基本パッケージに入っていないコマンドで
   v1tov2 v1のgrdファイルをv2のgrdファイルに変換する
   v2tov1 v2のgrdファイルをv1のgrdファイルに変換する
   makecpt v2のカラーパレットファイルを作成する
   coastmanager pscoast用のデータベースを作成する
   minmax アスキーテーブルから最大値最小値を抽出する
    psto24 Openwindows V.3の下でPostScriptファイルを24bitSunrasterファイルへ変換する
   grdraster SOESTのデータベース,ETOPO5,HaxbyのFAA,Geoid,Ageからgrdデータを抽出する(*)
   psvelomecha Institut de Physique du Globe de Parisで作成された発震機構解のPostScriptのプロットプログラム
   psnetwork Institut de Physique du Globe de Parisで作成されたDelaunay三角網の作成プログラム
   marvaridj 英Leads Univ.で作成されたOpenwindowsで動くgmtファイルのヴィジュアルエディタ(*)
(*)を付けたコマンドは非公開のプログラムです.

これらのコマンドの使用法はmanコマンドを使ってオンラインマニュアルを表示して調べて下さい.またコマンド名だけを入力したり,不法なパラメータをつけて入力すると同じマニュアルが表示されます.エラーの時にこのマニュアルが表示されると煩雑すぎるのでv. 3では間違ったパラメータを中心として表示されるようになりました.短い形式のオンラインマニュアルがコマンド名の後にスペースとハイフンを付けることで表示できます.

GMTの使用法の例

GMTの入門のためにはThe GMT-SYSTEM Technical Reference and Cookbookを通読し,例を試してみるのが良い方法ですが,この文書にはv.1の主体であったgmtファイルに関する事項が解説されていません.そこで,ここではまず航海ファイルを使う通常の作業の例によってGMTの使用法の概要を示します.
この例でのデータは白嶺丸の1985年のGH85-2航海のデータです.この例ではいくつかのコマンドを使用していますが,それぞれいくつかのパラメータを付けています.そのいくつかは必ず必要なもので,あるものはオプショナルなものです.この例で示したほかにもオプショナルなパラメータがあることが一般的ですが,詳細についてはマニュアルを参照して下さい.これらの例が理解できれば他のコマンドも試しながら修得していけると思います.
v.2ではパラメータ記号は大文字を使い,パラメータの記述中の文字は小文字が普通です.例えば出力の図法を指示するオプションは -Jm1 のように書きます.J は図法指定を意味するパラメータでプロットの際には常に必須のパラメータです.必ず大文字で書きます.m はメルカトル図法を示す記号で小文字を使用します.1は緯度で1度分の距離を出力図上で1インチのスケールで描くという意味です.
パラメータを記述するときはハイフン(-)を書いてからそのパラメータの関係の指定はスペースを入れないで続けて書いて下さい.例えば -G10g5 のようにします.これは,緯度経度などのアノテーションを10゜おきに書いてグリッド線を5度おきに引くという指定です.また,狭い範囲の図で10'おき5'おきにアノテーション,グリッドをそれぞれ描くように指定したいときは -B10mg5m としてください.m は分を示しているわけです.

1)ワークステーションへのデータファイルの取り込み

海洋地質部で使っているCR80形式の航海ファイルGH852A.DATA,GH852B.DATA,GH852C.DATAをRIPSのアーカイブなどからFTP や磁気テープなどを使って取り込みます.

2)GMTで用いるASCIIフォーマットへの変換

これらの白嶺丸の航海ファイルはCRファイルフォーマットで書かれているので,GMT v.1のASCII航海ファイル形式に変換します.この航海ファイル形式についてはAppendixで記述しておきます.
この変換には gh2dat というlocalなコマンドが使えます.
CRファイルGH852A.DATA,GH852B.DATA,GH852C.DATAを次のように,

%gh2dat gh852a.data gh852a.dat
%gh2dat gh852b.data gh852b.dat
%gh2dat gh852c.data gh852c.dat

と入力することによって全部で5分くらいでGMTのASCII航海ファイルgh852a.dat,gh852b.dat,gh852c.datに変換できます.

3)gmt フォーマットへの変換

次にdat2gmtコマンドで,GMTで用いる基本フォーマットであるgmt 形式に変換します.

%dat2gmt gh852a.dat gh852a.gmt
%dat2gmt gh852b.dat gh852b.gmt
%dat2gmt gh852c.dat gh852c.gmt

今度も全部で5分程度かかります.これで*.gmtファイルができました.このgmtファイルと,後で出てくるgrdファイルはバイナリファイルなので直接内容を見ることはできません.それでgmt2dat,grd2xyzのような逆変換コマンドや,gmtinfo,grdinfo等の情報抽出コマンドが用意されています.
現在白嶺丸の地調,石油公団のほとんどの航海と,STARMER関連の「かいよう」,「よこすか」の航海ファイルはgmtファイルにして,/gmt/data/gmtfiles/の中にまとめて入れてあります.

4)航跡図の作成

これらのgmtファイルを使って次のようにgmttrackコマンドを用いると航跡図のプロットができます.pscoastコマンドで陸の部分を黒く描くことができます.

%gmttrack gh852a gh852b gh852c -R129/133/33/36.5 -Jm1.5 -B1g0.5 -P -K > post
%pscoast -R129/133/33/36.5 -Jm1.5 -G0/0/0 -P -O >>post

これは1行目の入力で3つのgmt航海ファイル (gh852a gh852b gh852c)を用いて航跡を描きます.経度129E,133E,緯度33N,36.5Nの範囲で(-R129133/33/36.5),メルカトル図法で1度あたり1.5インチのスケールで(-Jm1.5),枠の数値を1゜毎に書き,0.5゜毎にグリッドを描く(-B1g0.5),それをポートレートモード(用紙を縦に使用する )で出力する(-P),そして,その出力はその後の描画データが追加して重ね書きされる形式である(-K),と指示しています.さらにオプションで指定して一定時間おきに測線にtickを入れたり,時刻をプロットしたりすることもできます.
2行目の入力では同じ範囲に黒く(-G0/0/0)ぬった陸の部分を重ね書き(-O)します.この行は前の行と同じプロットパラメータを使っているので,実際には-R129/133/33/36.5 -Jm1.5の部分は-R -Jmと省略することが可能です.このコマンドのための海岸線のデータは世界中のデータが用意されています.v.2までは狭い範囲で使う場合にはあまり正確とは言えないものでしたが,v3では粗いものから細かいものまで5段階の海岸線データベースが使えるようになりました.
この2行のコマンドを入力すると,航跡図のPostScriptファイルpostが作られます.このファイル名は適当なものでよいのですが求める出力が得られるまで何度もやり直すのが普通なので,簡単で一次的なファイル名を使うのがよいでしょう.
GMTのコマンドのパラメータの与え方に関する一般的な注意に関してAppendixで述べておきます.

5)データの編集

gmttrackと同じような方法で使えるコマンドとして測線に沿って異常値をプロットするpswiggleがあります.このような測線図,異常値のプロットで不良データを発見したら,前述のgmtで用いるASCIIファイル上に戻ってエディタを用いて編集し,またプロットするという作業を繰り返す必要があります.ASCIIファイル上での編集ではviエディタを用いて,検索,置換等が高速に実行できます.その他UNIXのユーティリティAWKが利用できるでしょう.AWKについてはAppendixに簡単に使い方を紹介しています.

6)PostScript ファイルの表示

SUN Solaris 1 (SunOS 4.1.x)のOpenwindows を使っていればポストスクリプトシェル psh を使ってこれを直接画面に表示することができます.
%psh post
また,pageview,imagetool,ghostviewなどのOpenwindows,Xなどのウインドウシステム上のビュアーを使って表示することもできます.これらに関する詳細をAppendixに紹介しておきます.

7)PostScript ファイルのプリント

gsjsant,gsjgmtでは,次のように lpr コマンドでワークステーションに接続されている SPARCprinterに出力することができます.このプリンタはA4用紙しか使えませんが,普通の図であれば20秒から1分程度で400 dpiの高品質のプリントの出力が可能です.

%lpr post

SPARCprinter の場合は,PostScriptファイルやSunrasterファイルは自動的に識別してプリントできます.SUN OS4.1.3ではプリントキューの表示は lpq を使用します.プリントコマンドの取り消しは prm を使用します.
白嶺丸上ではHPのレーザープリンタがsun10からもネットワークを通して利用できます.コマンドは lp を使用しますが,ここで,気を付けなければならないのが,プリンタモデルを指示するオプションです.

%lp -d ps post

とします.この-d psをつけないと,通常のASCIIデータとみなされ膨大なリストを出力します.間違いに気づいたらすぐにプリンタの電源を切り,cancelコマンドを実行します.すでにデータがHPに行ってしまった場合はhp1上でrootでSAMでプリンタマネージャを起動して,プリントキューのリストを出して,出力のキャンセルをします.プリンタを再起動するときには,データを消すためにリセットボタンを押しながら電源を入れ直します.
ファイルをマッキントシュ上に転送して,DropPS,SendPS等のプログラムを利用して,Laerwriter,MicrolinePSのようなPostScriptプリンタにプリントすることもできます.Microlineを使えばB4の出力が可能です.ただし,非常に時間がかかります.

8)航海をまとめて扱う方法

上の例では3つの航海ファイルの航跡をまとめて表示しました.航海数が多いとコマンドの後にgmtファイル名をならべて書くのが大変です.特に決まったディレクトリにgmtファイルをまとめてあるとそのパスも書かなければならないので,扱うファイル数が多くなるとコマンド行が長くなりすぎるてしまうので次の方法が必要です.
航海名を記述したファイルを作り,コマンドの中でUNIXのcat コマンドを利用してその内容を求めて取り込む,というものです.具体的にはファイル852.listをエディタで作ります.Appendixにエディタの簡単な使用法とともにこのファイルの作成方法を示します.
次にgmttrackなどの航海名を入力するタイプのコマンドのパラメータの航海名部分を次のように記述します.ここで使われている引用符は通常のsingle quotation markではなくback quatation markです.
'cat 852.list'
先のgmttrackの場合は次のようになります.

%gmttrack 'cat 852.list' -R129/133/33/36.5 -Jm1.5 -R1g0.5 -P -K >post

9)xyz ファイルの作成

gmtファイルが存在していれば位置と水深,位置と重力異常,位置と地磁気異常と行ったxyzの組み合わせのデータを抜き出してファイルを作ることができます.この次の段階のグリッドファイルを作るためにはこのxyzファイルが必要です.

%gmtlist 'cat 852.list' -Fxyg -R129/133/33/36.5 > 852.xyg

すると,52.listに記述されたgmtファイルから,位置と重力のデータを抽出したファイル852.xygを作成することができます.-Fxyg のかわりに -Fxym で位置と地磁気の,-Fxytで位置と水深のファイルを作成できます.その他船速,進行方向を選択することも可能です.
xyz ファイルはアスキーファイルでエディタで確認,編集が可能です.

10)グリッドファイルの作り方

グリッドファイルの作成には A)グリッドのブロックの平均値または中央値をとる前処理,B)グリッド化処理の2段階の処理が必要です.GMTのグリッド化のコマンドは非常に協力で,ランダムデータ,データ密度の差が大きなデータに有効です.グリッド化の際にはデータの変化の大きな地形などの場合と,変化の小さな重磁力などの場合とでテンションファクターを調整することが可能です.
A)まず blockmean または blockmeddian コマンドでグリッド化の前処理を行います.

%blockmean 852.xyg -R129/133/33/36.5 -I1m > 852.1x1.xyg

これで位置と重力異常のxyzファイルを用いて先に測線図を描いた同じ領域の,1分きざみのグリッドブロックを定義して,そのグリッド毎にデータがあれば位置と値の平均をとりファイル 852.1x1.xyg に出力することができます.blockmedian コマンドも同じような目的に用いるもので平均値の代わりに中央値をとるものです.
B)次にグリッド化のコマンドsurface を用いてグリッドファイル *.grdファイルを作成します.

%surface 852.1x1.xyg -R129/133/33/36.5 -I1m -G852g.grd -T0.2 -C0.1

これで先に作成したブロック内で平均化されたxyzデータを用いてグリッドファイルを作成します.このコマンドで使われる -T -C についてはAppendixで説明します.
このコマンドでは出力ファイルを -G という記号のパラメータを用いて記述します.これはこの出力がバイナリファイルだからです.他の多くのコマンドでは出力がアスキーテキストなので標準出力に出されるため > 記号を使ってファイルにリダイレクトする使い方が一般的です.出力が標準出力であると,何も指示をしない場合には画面に出力されますし,必要なら別のコマンドの入力としてパイプで送り込むことも可能です.なお >> とすると,そのファイルに追加して書き込むことになります.

11)コンター図の作成

グリッドファイル *.grd を使用してコンター図を描きます.

%grdcontour 852g.grd -C50 -Jm1.5 -B1g0.5 -P -K > postg
%psmask 852.1x1.xyg -R129/133/33/36.5 -Jm1.5 -B1g0.5 -I1m -C0.2 -N -P -O -K >> postg
%pscoast -R129/133/33/36.5 -Jm1.5 -G0/0/0 -P -O >> postg
%psh postg

まず grdcontour コマンドで入力を*.grdファイルと指定して,コンター間隔を50(単位はz,ここではmgal)とします.後は投影図法とスケール,枠の記述と緯線経線のグリッド間隔,ポートレートモード,後に PostScript ファイルが重ねられるという指定,です.
次の psmask コマンドはデータのない部分のコンターを隠すためのコマンドで,グリッド化されたxyz ファイルが必要です.範囲,図法とスケール,枠の記述と緯線経線のグリッド間隔,グリッドブロックの大きさを記述した後,影響半径というパラメータ(-C 単位はグリッドの数)をセットします.この半径以内のデータを信頼できるものとします.-Nはこれで指定された領域以外をマスクするというオプションです.これを指定しないと,意味のあるデータ領域のみが隠されます.-Kと-Oでは前のファイルに重ね書きをし,この後で重ね書きをされることを指定しています.
次の pscoast は測線図の場合と同じで,陸地部分を黒く塗りつぶします.
その結果作成されたPostScript ファイルpostgを画面で確認しています.
この例ではコンターの値がかけませんでしたが,描くコンターの値とコンターの値をプロットするかどうかを指定するファイルを作成しておいてそれを取り込むことができます.そのファイル形式はAppendixに示します.

%grdcontour 852t.grd -Ccont.t -Jm1.5 -B1g0.5 -P -K > postt

以下 psmask, pscoastは重力の場合と同様です.ここで -C パラメータはすぐに数値を記述するとコンター間隔を意味しますが,この例のようにファイルにすれば描くコンターを指定できるのです.当然ながら数字だけのファイル名を使わないようにしなければなりません.画像ファイルの作成のところで述べるcptファイルを使ってアノテーションを伴わないコンターを描くこともできます.

12)鳥瞰図の作成

%grdview 852t.grd -Ccol.cpt -Jm1.5 -E230/70 -Qs -B1g0.5 -P > postt

-Qsによりサーフェスモードで鳥瞰図を描きます.色はカラーパレットテーブルファイルcol.cpt で指定しています. 見る角度は230゜方向で,水平から70゜の高度となっています.デフォールトではモードはメッシュプロットモード-Qmです.-Qiとすればイメージをプロットします.さらに,必要なら-Iintensityfileで濃淡を指定できます.
カラーパレットテーブルファイルのフォーマットについてはAppendix (A-11)で述べています.

13)画像ファイルの作成と表示,プリント

%grdimage 852t.grd -Ccol.cpt -Ishade -Jm1.5 -B1g0.5 -P >postt

-Ccol.cpt で指定されたカラーパレットテーブルファイルでカラー画像を作成します.
grdviewと同様に,-Iintensityfileでの濃淡を指定できます.ある方向から光を当てたときの反射強度を与えるコマンドgrdgradientを使って強度ファイルを作成しておきます.
%grdgradient 852t.grd -A270 -Gshade

14)カラー画像の扱い方

grdimage,grdviewなどではカラーポストスクリプト画像を作ることができます.作成した画像はポストスクリプトビュアープログラムで画面に表示したりAppendix (A-9),カラーポストスクリプトプリンタで印刷することができます.
経験上3M-RainbowでPortrait (縦位置)で印刷する場合に-Yのデフォールトの値がおかしいことがあるので,-Yを調整する必要があります.

その他の機能など
塗りつぶしパターンの使い方
-G で指定する.

15)文字のプロット

pstextで座標,座標基準,角度,フォント番号,フォントサイズ,表示文字列を内容とするファイルを読んでプロットします.非常に面倒だと思うでしょうが,できあいの作図,グラフ作成,visualizationなどのソフトでは任意の場所に任意の文字を描くことができないものが多いようですが,GMTでこの機能を使えば任意の場所に必要な文字を入れることができます.また,一旦作図のためのスクリプトを書いておけば再プロット,再利用ができるのでかなり便利です.
(未完の主題)
GMTで使える地図投影法
ETOPO5データの利用法
衛星アルチメトリによるグローバル重力異常データの利用
http://www.aist.go.jp/GSJ/dMG/free/gravityを見て下さい.
グローバル縞状地磁気異常データベースとそのプロット
グローバル海底年代ファイルとそのプロット
プレート境界ファイルとそのプロット
http://www.aist.go.jp/GSJ/dMG/free/platesを見て下さい.
psvelomechaによる発震機構解のプロット
世界の地殻熱流量データベースのプロット
pslibの使い方

Appendix

A-0)GMTが開発された経緯

このシステムはLamont-Doherty Geological Observatory, Columbia University(現在はEarth Observatoryと変更されている)で1960年代にTalwaniらが構築し,その後Walter Brownが作り替えて,多くの人が使ってきたというBrown Bookというデータベースシステムをネットワーク,ワークステーションといったコンピュータ環境の変化に対応させるために主にPaul Wesselが作ったものです.Wesselは重力探査の専門家でLDGOでPhDを取得し,現在University of Hawaii, SOESTのProfessorです.現在の共著者のWalter Smithはやはり海域の重力の研究者でScripps Institute of Oceanographyを経て現在はNOAAのNGDCに在籍しています.LDEO,SIO,SOESTでは日常的に使われていると同時に,利用法,バグ情報の交換がなされ現在もバージョンアップが進められています.現在のバージョンはv.3です.

A-1)GMTはどこから入手するか?

GMTを新しくインストールする場合に必要なファイルは次のところからネットワーク経由で入手できます.
Anonymous FTPサイトkiawe.soest.hawaii.edu (IP 128.171.151.16)の/pub/gmtディレクトリにGMT.tar.Zという名前のファイル名でおいてあります.このディレクトリにはGMT_scripts.tar.Z,GMT_supplemental.tar.Zというオプショナルなスクリプトのサンプルやプログラムもおいてあります.
GMTを動かすためにはnetCDFというグリッドファイルの操作のためのライブラリが必要です.これはAnonymous FTPサイトunidata.ucar.edu(IP 128.117.140.3)にnetcdf.tar.Zという名前のファイルとしておかれています.

A-2)GMTとUNIX

GMTはUNIXワークステーションの上で使います.よくも悪くもUNIXの特徴を持っています.UNIXはMS-DOSと似ています.SUNなどのUNIXワークステーションはマッキントッシュにも似ています.
GMTはUNIXのコマンド群から成り立っていて,それらの使用法はUNIXの通常のコマンド(ツール)を使うのと良く似た使い方になっています.つまりGMTという特定のアプリケーションの中でのみ使用する操作法を学習する必要はありません.MS-DOSは普通「一太郎」や「松」などのアプリケーションプログラムの土台として動いていて,ユーザーがコマンドラインからDOSを使うのは「松」を起動するときにmatuと入力するときだけだったりします.UNIXえは普通ほとんどの仕事をコマンドラインで実行します.アプリケーションの上でも同じようにコマンドラインのように入力するインタフェースを用いるのが一般的です.
GMTのコマンドはUNIXのプロンプト%が出た状態でパラメータとともに入力し,実行が終わるとまたUNIXのプロンプト%が出た状態に戻ります.このような方法のためプログラムが必要なパラメータはすべて最初に与えます.これによって,一つのコマンドの中でプログラムから入力を促したりメッセージを出力する必要がなくなり,コンソールを使用しないですむので,マルチタスクの機能を利用してバックグラウンドで何本もコマンドを実行することができます.書くコマンドはたくさんのオプションを必要として長い入力が必要な場合が多くなりますが,UNIXのcsh (MS-DOSのcommand.comに対応するshellプログラム,キーボードからの入力を解析してプログラムに引き渡すプログラム)の機能や他のツールを使うことによって仕事がしやすくなっています.
UNIXはキャラクタベースのOSですからネットワークにつながったどの端末からでも使うことができます.GMTも所内LANにつながっているか,海洋地質部や白嶺丸のローカルトークのようにLANとゲートウエイを通じてつながっていればどのマッキントッシュからでも利用できます.
いずれにしても,GMTを使うためにはUNIXの初歩的な知識は不可欠です.新しいシステムを学習するのは多かれ少なかれおっくうなものですが,この機会にUNIXを使い始めてみてはいかがでしょう?

A-3)GMTv.1の航海ASCIIファイルフォーマット

この形式は
1行目が航海の始まった年,データレコードの数,航海の名称を書いた行,
2行目以降は次の形式の1行1レコードがデータ数行分続いたものです.
データの形式 yyyy/mm/dd/hh:mm:ss lat lon g m t
緯度経度は度の単位で与えます.gは重力異常値(mgal),mは地磁気異常値(nT),tは測深値(m).
ファイルの先頭部分の例を示せば

1977 7511 c2104
1977/12/17/16:44:00 13.00167 300.57500 NaN NaN NaN
1977/12/17/16:48:00 13.00024 300.58503 NaN NaN -290
1977/12/17/16:50:00 12.99952 300.59005 6.10 NaN NaN

ここで NaN というのはデータが使えない場合を示すものです.測深値は - が必要なことに気を付けて下さい.

A-4)AWKの使い方の例:不良データの検出

AWKはUNIX上でよく使われる簡易言語の一種で,1行のプログラムでもそれなりの仕事ができるというのが特徴です.この例のような表形式のデータに対して,ある条件を設定して1行毎に内容を調べ,条件に合致した行に何らかの操作を行うというフィルタプログラムです.条件を記述したものをパターン,操作を記述したものをアクションと呼びます.普通のAWKのプログラムはパターンとアクションの組み合わせからなります.それを次のように記述します.

awk ' パターン { アクション } ' 適用するファイル名

これで出力はアクションの結果で標準出力である画面に表示されます.ファイルに取り込みたい場合は >temp のように出力をリダイレクトします.
例えば重力異常値として70より大きな値を検出しようとすれば,次のようにしてみます.

%awk '$4 >700 {print $0}' 852.dat > check.g

ここで $4 は4番目の列のデータの内容を示す変数です.また $0 は行全体を意味します.852.datデータは日時,緯度,経度,重力異常,地磁気異常,水深の6列の数値からなていますから,$4 は重力異常値となり,$4>700 で重力異常値が700より大きいという条件を記述したパターンということになります.アクション部分は print $0 なのでパターンに合致した行の内容が全部書き出されることになります.日時と重力異常値だけ書きたいときには print $1, $4 とします.
パターンもアクションも細かい指定が何重にも可能でif文,for文,フォーマット付き印刷などの機能もあります.変数を用いて’ある列の値が前の行の値と比較して差が基準より大きければ印刷する’ということも簡単にできます.プログラムが長くなって1行に書きにくい場合にはエディタを用いてプログラムを書いて次のように使用します.例えば aw.prog を作って実行してみます.

%awk -f aw.prog 852.dat > check.g

aw. progの内容

#test program to extract anomalous data by AWK
$2>90 | $2<-90 | $3<-180 | $3>180
{printf "illegal location data at line %d \n", NR}
$4 > 700 | $4 < -700 | $5 > 1500 | $5 < -1500 | $6< -10000
{printf "possible bad data at line %d contents of GMT are %8.2f %8.2f %8.2f", NR, $4, $5, $6}
このプログラムを実行すると,緯度-90から90,経度-180から180以外のデータの行数,重力異常,地磁気異常の絶対値がそれぞれ700以上または1500以上,または水深が 10,000 m以上のデータの行数と値がプリントされます.

なお,AWKにはマッキントシュ版,DOS版などもあるようです.AWKに関しては開発者が書いた次の教科書があります.
The AWK Programing Language by A.V. Aho, B.W. Kernighan and P.J. Weinberger, 1988, Addison-Wesley Pub. Co.
プログラミング言語AWK 足立高徳訳,1989,アジソン・ウエスレイ・トッパン

A-5)viエディタの使い方の例:航海リストファイルの作成

viはネットワーク上のマッキントッシュなど,どの端末からでも使用できるスクリーンエディタです.入力モードとコマンドモードの区別が画面上でできないので最初は使いにくいかも知れませんが,慣れてくると非常に便利なエディタであることがわかってくると思います.
vi 852.list(リターン)と入力し852.listというファイルを編集するようにvi エディタを起動する.ファイル名無しで起動して書き込みの時に指定してもかまわない. a (キーを1回入力)で入力モードにする.これ以降押したキーはファイルの内容になります.
gh852a(リターン)
gh852b(リターン)
gh852c(リターン)
と入力した後エスケープキーを押して入力モードを終了してコマンドモードにします.コマンドモードではカーソルを動かす,行単位でコピー,削除するといった各種のコマンドが使えます.ここまでの作業で作られた内容はメモリ上と,バッファファイルに存在しているので,最後にファイルに書き込みます.このときはコマンドモードからさらにexモードに入って書き込みの命令を実行します.
:wq (3つのキーを順に入力)とするとファイルに書き込んで,エディタを終了します.ここで : は exモードへはいるための記号,w は書き込み,q はエディタの終了のための ex コマンドです.
viの起動時にファイル名を指定しなかったときには,:w cruise.listのように引数として書き出すファイル名を記述します.その後終了するなら,:q を入力します.
入力モードとコマンドモードの区別が付かなくなったらエスケープキーを押してみて,確実にコマンドモードにして,その後の作業を続けます.コマンドモードでの1文字削除は x,1行削除は yy です.入力時の入力データの変更はバックスペースキーで行います.

A-6)GMTのコマンドのパラメータの与え方に関する一般的な注意

GMTのコマンドのパラメータは与え方の順番は決まっていません.
出力はグリッドファイルの出力をするコマンド以外は標準出力に出されます.図面のプロットの場合には,通常,一連のコマンドを実行して一つのプロットを作成する場合,最初のコマンドで-Kを付けて出力を > を付けてあるファイルにリダイレクトし,引き続くコマンドで -K -O を付けて >> そのファイルに追加していきます.最後のコマンドで -O を付けて最後の出力を >> で追加します.
引数の内で -(ハイフン)がついていないものは入力データファイルとみなされます.
 -(ハイフン)がついたパラメータは大文字で種類が示されます.パラメータの内容はそのパラメータにスペース無しで続けて記述します.
よく使われるパラメータを簡単に示します.
-B プロットする図の枠に付けるアノテーションとグリッドの間隔を指定する.
-C コンター間隔,コンターファイル名,カラーパレットテーブルファイル名を示す.
-G グリッドファイルを出力するタイプのコマンドでは出力ファイル名を示す.
-G 図をプロットするタイプのコマンドではは塗りつぶしのグレイレベル,RGBによる色,またはパターンを示す.
-H 入力アスキーファイルがヘッダを持っていることを示す.行数を付ける.
-I グリッドの間隔を示す.
-J 投影図法を示す.
-K PostScriptの出力がこの後にさらに続くタイプであることを示す.
-O PostScriptの出力が別に存在するファイルにつけ加えられるタイプであることを示す.
-P 用紙を縦で用いることを示す.付けない場合のデフォールトは横である.
-R プロットするデータの範囲を示す.常に-RminX/maxX/minY/maxYの順に指定する.地理的データでは西ー東ー南ー北の順になる.経度ー緯度の順なので注意する.
-U GMTのタイムスタンプをプロットすることを指示する.
-V 処理の途中の情報を表示するように指示する.
-W プロットする線の太さと色,タイプを示す.
-X プロットするデータの原点のX座標のオフセットを示す.デフォールトは1.5である.一度設定するとその後のコマンドでも同じオフセットが有効である.
-Y プロットするデータの原点のY座標のオフセットを示す.デフォールトは1.5である.一度設定するとその後のコマンドでも同じオフセットが有効である.
-: 入力ファイルが(x,y)でなく(y,x)であることを示す.緯度,経度の順になっているデータのプロットなどで使用する.
-# プロットのコピーの枚数を示す.

標準的なパラメータ-B -J -K -O -P -R -X -Y -# -:に関しては一度使用され設定されると,そのディレクトリに.gmtcommandsというファイルが作られ,記録される.引き続くコマンドでも同じパラメータで使用するときには-R -Jmなどとすれば詳細は以前の設定がそのまま用いられる.

A-7)gmtdefaults

GMTで用いるデフォールトパラメータはgmtdefaultsというコマンドで表示することができる.この出力を編集してから,.gmtdefaultという名前でセーブして,ホームディレクトリにおいておけば自分の専用のデフォールトパラメータを設定しておくことができる.

A-8)surfaceコマンドのパラメータ -Tと -C

surface コマンドはランダムに分布する xyz データを入力して,次の式を解くことによって z(x,y)を求めて*.grdファイルを作成します.

(1-T)*L(L(z))+T*L(z)=0

ここでTは0から1のTension Factor,LはLaplacian Operatorです.T=0 の場合はSuperMISP, ISMパッケージの"minimum curvature "に対応する場合となります.その場合には好ましくない振動を生じて異常なデータを作成する可能性があります.そのためT>0を用いるのがよいでしょう.経験的には T=0.25 程度がポテンシャルデータの場合に良さそうであり,急崖のある地形データでは大きく T=0.35 程度がよいようです.T=1では調和面が得られ,データポイント以外では極大,極小値を取らなくなります.デフォールトはT=0です.
surface コマンドをblockmean,blockmedianを適用しないで,実行すると空間的なエリアシングを起こすことになります.
もう一つのパラメータ -Cは convergence limitというもので,繰り返し計算(itaration)を終了させる条件を記述するものです.計算の結果どのグリッドの値もこの値以下の最大値での変化しかしなくなったら収束したとみなします.単位は z と同じです.このパラメータは省略可能です.デフォールトは値の勾配の1%です.

A-9)PostScriptビュアープログラム

Openwindowsv.2に付属している pageview というアプリケーションプログラムを使えばpostファイルを読み込んでプリンタ上のイメージを画面上で確認することができます.
Solaris2.1以降ではOpenwindowsv.3に付属している Imagetool という非常に強力な画像ファイルビュアーを用いて,PostScriptの画像ファイルを表示することができます.Imagetool は主な画像ファイルフォーマットをサポートしており,PostScriptファイルをSunraster,Tiff等のラスターファイルに変換することも可能です.SUN OS4.1.3上では英語版のOpenwindows v.3がgsjsant上で利用でき,psto24のようなPostScript画像からラスターファイルへの変換ができますが,Imagetoolはありません.
X11R5などのXwindowを利用する場合にはPostScript互換のフリーウエアであるGhostScript,Ghostviewが利用できます.ただし,表示速度が遅く,印刷の品質が悪いのが欠点です.GhostScriptはマッキントッシュ用もありますが,表示速度は非常に遅くGMTで作るような大きなPostScript画像の表示には向きません.
NeXTはDisplay PostScriptを使っており,Imagetoolと似たビュアープログラムが利用できるようです.

A-10)cont.tの形式

cont.tの内容は次のようにします.1列目にはコンターを描くz の値.2列目には C か A を記述します.Cはコンターを描く,Aはコンターを描き,アノテーションも書くという意味です.水深 0 から 200 m までを5m 間隔でそれより深い部分は,20m間隔で粗く描き,50m,100m,200mのみには値をプロットするように指示しました.水深値には - を付けること,値は小さい方から書く必要があります.

-----------
-240 C
-220 C
-200 A
-195 C
----------
-55 C
-50 A
-45 C
----------
-15 C
-10 C
-5 C

A-11)カラーパレットテーブルファイルについて

カラーパレットテーブルファイルは z の値を n-1段階に分けるとすると,zの値の小さい方から順に次のように記述する.RGBカラーモデルを使用しているとき,

z0 Rmin Gmin Bmin Z1 Rmax Gmax Bmax
...
zn-2 Rmin Gmin Bmin Zn-1 Rmax Gmax Bmax
F Rfore Gfore Bfore
B Rback Gback Bback

ここでRminなどは例えばz0からz1の範囲での最小値にあたえる赤の成分を示し,Rminは最大値に与える赤の成分の強度を示します.グレイスケールを得たいときにはRGBの値を同じにします.個々のレンジでmaxとminの値を同じにすればそのレンジは一定の色となり,異なる場合にはレンジの途中の値は線形補間された色となります.FとBはテーブルで指定した最大値より大きな,および最小値より小さなデータに対して色を指示するものです.指定しなければデフォールトが使われ普通 F が白,B が黒となています.カラーモデル,デフォールトのForegrouund, Background colorなどは.gmtdefaultsに記述しておくことができます.

A-12)gmthelpメーリングリストへの参加

GMTのユーザーの間で,ノウハウを交換し合うためのgmthelp@soest.hawaii.eduというメーリングリストがあります.たしかsubscribeという内容のメールをこのアドレスあてに送ると登録されたと思います.
以上(未完)