水環境保全チーム

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した全世界１km解像度のグリッドデータ名を設定 rem set文でカラーパレット名の変数cpoにカラーパレットデータGMT_globe.cptを設定 rem set文で１度あたり何センチで描くかを決める変数ｒを設定（ここでは関係ない） rem set文で描く範囲を示す変数wakuを設定。左端/右端/下端/上端の順。 rem set文で緯度経度枠のデザインを指定。 a：何度毎に数値を入れるか。ｆ：何度刻みで枠に目盛を入れるか。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などの閲覧しやすい画像に変換しておくことで、その後の処理が楽です。　都市名をプロットするスクリプトは上記のスクリプトに以下の行を付加し、下記の都市データを同一フォルダ内に置けばＯＫです。コマンドの詳しい説明は省略します。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 　　　　　　　　・　　　　　　　　・

寒地土木研究所ホームページへ

サイトマップ