GMTで作る世界地図

Albers図法による日本周辺地図

日本周辺地図

使用データ:世界の陸地と海洋地形1分(解像度約1km)グリッドデータは下記URLからダウンロード可能です。全データは450MBくらいあるので、緯度経度で利用範囲を指定し、部分的にダウンロードすることも可能です。
http://www.bodc.ac.uk/data/online_delivery/gebco/ 

日本〜台湾にかけて陸地から海底までを含めた地図を描くスクリプト(rem文は注釈行です)

set psfile=Output_image.ps
set grddata=GridOne.grd
set cpo=GMT_globe.cpt
set r=0.2
set waku=120/150/20/50
set bp=a10.0f5g0.0
set alb_f=135/35/30/40/1:25000000
psbasemap -Jb%alb_f% -R%waku% -B%bp%/%bp% -P -K > %psfile%
grdimage %grddata% -Jb%alb_f% -R%waku% -B%bp% -C%cpo% -T -P -O -K >> %psfile%
pscoast -Jb%alb_f% -R%waku% -B%bp%/%bp% -W1 -Dh -V -P -O -K >> %psfile%
rem remは注釈文です。
rem set文で出力するポストスクリプトファイル名の変数psfileを設定。このポストスクリプトファイルにイメージが追加されていく。
rem set文でグリッドデータ名の変数grddataにGEBCOからDLした全世界1km解像度のグリッドデータ名を設定
rem set文でカラーパレット名の変数cpoにカラーパレットデータGMT_globe.cptを設定
rem set文で1度あたり何センチで描くかを決める変数rを設定(ここでは関係ない)
rem set文で描く範囲を示す変数wakuを設定。左端/右端/下端/上端の順。
rem set文で緯度経度枠のデザインを指定。 a:何度毎に数値を入れるか。f:何度刻みで枠に目盛を入れるか。g:何度刻みで格子を描くか(ここでは0なので図内に格子が入らない)。
rem psbasemap:画像の下地の設定。-JbのbはAlbers図法で投影されることを示します。
rem grdimage:grddata変数で定義されたグリッドデータを指定された情報で投影します。
rem pscoast:指定された情報で海岸線を投影します。
rem 変数名の前後に%を付けるのはWindowsでシェルスクリプトを動かすときの約束事らしいのですが、詳しいことはよく知りません。ちなみに、カレント ディレクトリの隣のディレクトリ(test)のデータファイルを参照するときはset grddata=..\test\Gridone.grdとなり、単に上位ディレクトリのファイルを参照するときは、set grddata=..\Gridone.grdとなります。

以上のスクリプトを実行すると、Output_image.psというファイルが出来上がります。
このファイルのイメージを見るには、ポストスクリプトファイルを閲覧可能なGsviewやAdobe社製イラストレーターなどのソフトが必要です。

エッケルト図法による日本付近を中心とした世界地図

日本付近を中心とした世界地図

世界地図を描く場合、上記の1km解像度のデータだと処理が長引く恐れがあるので、5km解像度くらいのものが適当だと思います。それでも、下記スクリプトで得られる結果は350MBくらいのサイズになります。5km解像度のデータはNOAAの以下のサイトからDLできます。世界の地形と海洋のデータはこのサイト中、カバレージがGLOBAL_LAND_and_OCEANSのもので、Both(両方)となっているリンクところをクリックするとDLできます。(ただし、ftp経由なのでネット環境によってはDL出来ないこともあります)
http://www.ngdc.noaa.gov/seg/fliers/se-1104.shtml

rem 日本付近を中心とした地球展開図
set waku=-30/330/-90/90
set ekrp=kf150/0.06c
set ofn=glb_ekr.eps
set grd=tbase.grd
grdimage %grd% -R%waku% -J%ekrp% -CGMT_globe.cpt -P -T -K > %ofn%
pscoast -R%waku% -J%ekrp% -Bg30 -Ia/0/0/0/255 -W0 -P -O -K >> %ofn%

以上のスクリプトを実行すると、glb_ekr.epsというポストスクリプトファイルができます。PCのスペックにもよりますが、gsviewでイメージを開くのに数分〜10分程度かかります。イメージが完全に展開されたら、file -> convert でjpgなどの閲覧しやすい画像に変換しておくことで、その後の処理が楽です。

都市名をプロットするスクリプトは上記のスクリプトに以下の行を付加し、下記の都市データを同一フォルダ内に置けばOKです。コマンドの詳しい説明は省略します。gawk文のあたりはちょっと複雑そうに見えますが、パイプコマンドを示す"|"のところで左と右に分かれています。つまり、パイプの左で作った pstext用のデータをパイプで右に受け流しているのです。左からなにかがやってきて、右にそのまま受け流すワケですね。!?どこかで聞いたような、、、。

set wcn=w_city_m.txt
gawk "{print $3 , $2}" %wcn% | psxy -R%waku% -J%ekrp% -Sc0.1 -W1/0/0/0 -G255/0/0 -P -O -K >> %ofn%
gawk "{print $3 , $2, 9, 45, 1, 1, $1}" %wcn% | pstext -R%waku% -J%ekrp% -P -O >> %ofn%
rem pstextが必要なデータ並び: 経度,緯度,ポイント,角度,フォント,書き位置

※主要都市名をプロットする場合のデータ形式サンプル。  
以下の様な並び(都市名 緯度 経度)で都市のデータがある場合、地図上に都市名をプロットできます。都市名ファイルをテキスト形式でw_city_m.txt とし、以下の並びでデータが入っているものとします。
Abadan 30.350000 48.266667
Abidjan 5.316667 355.983333
Abu_Dhabi 24.280000 54.220000
Abuja 9.050000 7.320000
Accra 5.350000 359.940000