WANtaro雑記帳（Win_GMT：東アジア震央マップ）

GMT の pscoast と psxy コマンドを使って，東アジアの震央マップを作成した事例です．

準備

地震データは米国地質調査所のデータをダウンロードさせていただきました．

ダウンロードしたデータは以下のようになっています．


      Link to USGS Home Page 


  NEIC: Earthquake Search Results


                    U.  S.  G E O L O G I C A L  S U R V E Y

                     E A R T H Q U A K E  D A T A  B A S E



 FILE CREATED:  Sun Oct 30 05:53:38 2011
 Global Search   Earthquakes=    383175
 Catalog Used: PDE
 Date Range:   1900/01/01   to    2011/10/31
 Magnitude Range:   3.0  -  10.0
 Depth Range:    0  -  100
 Data Selection: Historical & Preliminary Data


 CAT    YEAR  MO DA  ORIG TIME   LAT    LONG  DEP  MAGNITUDE  IEM DTSVNWG DIST
                                                              NFO          km
                                                              TF 

 PDE    1973  01 01 034609.80  -9.21  150.63  41  5.3 mbGS    ... .......      
 PDE    1973  01 01 052229.80 -15.01 -173.96  33  5.0 mbGS    ... .......      
 PDE    1973  01 01 114237.50 -35.51  -16.21  33  6.0 MsGS    ... .......      
 PDE    1973  01 01 213558.60  13.90  -44.88  33  4.4 mbGS    ... .......      
 PDE    1973  01 01 214117.10  13.74  -45.27  33  4.4 mbGS    ... .......      
 PDE    1973  01 02 005320.30  -9.85  117.43  66  5.5 mbGS    ... .......      
 PDE    1973  01 02 015028.50  38.22   20.17  33  4.0 MLATH   ... .......      
 PDE    1973  01 02 022709.20   1.03  126.21  61  5.4 mbGS    ... .......      
・・・・・・・・・・
・・・・・・・・・・
・・・・・・・・・・
 PDE-Q  2011  10 28 081429.60   5.09   94.22  57  5.0 mbGS    ... .......      
 PDE-Q  2011  10 28 084809.40  38.49   43.77  10  4.6 mbGS    ... .......      
 PDE-Q  2011  10 28 090837.95 -40.63  126.42  10  4.7 mbGS    ... .......      
 PDE-Q  2011  10 28 091846.27  35.53  -97.37   5  3.5 LgGS    ... .......      

------------------------------------------------------------------------

*USGS National Earthquake Information Center*
USGS Privacy Statement  | Disclaimer
 	FirstGov

GMTへの入力用データとしては，以下のようにデータ部分だけに加工し，これを「inp_eqdat.txt」とします．

また，同様にM7以上の地震のみを取り出した「inp_eqdat_M7.txt」を作成します．

 PDE    1973  01 01 034609.80  -9.21  150.63  41  5.3 mbGS    ... .......      
 PDE    1973  01 01 052229.80 -15.01 -173.96  33  5.0 mbGS    ... .......      
 PDE    1973  01 01 114237.50 -35.51  -16.21  33  6.0 MsGS    ... .......      
 PDE    1973  01 01 213558.60  13.90  -44.88  33  4.4 mbGS    ... .......      
 PDE    1973  01 01 214117.10  13.74  -45.27  33  4.4 mbGS    ... .......      
 PDE    1973  01 02 005320.30  -9.85  117.43  66  5.5 mbGS    ... .......      
 PDE    1973  01 02 015028.50  38.22   20.17  33  4.0 MLATH   ... .......      
 PDE    1973  01 02 022709.20   1.03  126.21  61  5.4 mbGS    ... .......      
・・・・・・・・・・
・・・・・・・・・・
・・・・・・・・・・
 PDE-Q  2011  10 28 081429.60   5.09   94.22  57  5.0 mbGS    ... .......      
 PDE-Q  2011  10 28 084809.40  38.49   43.77  10  4.6 mbGS    ... .......      
 PDE-Q  2011  10 28 090837.95 -40.63  126.42  10  4.7 mbGS    ... .......      
 PDE-Q  2011  10 28 091846.27  35.53  -97.37   5  3.5 LgGS    ... .......

実行用バッチファイルと解説

gmtset ANOT_FONT_SIZE 10
set fig_out=fig_eq_east.eps
gawk -f awk_klegend.awk inp_legend2.txt > legend.bat
pscoast -R90/150/0/45 -JM12c -Ba15g15 -G200 -S255/255/255 -N1 -P -K > %fig_out%
gawk "{print $6,$7}" inp_eqdat.txt | psxy -: -R -J -B -Sc0.1 -G255/255/255 -W1/255/0/0 -P -K -O >> %fig_out%
gawk "{print $6,$7}" inp_eqdat_M7.txt | psxy -: -R -J -B -Sc0.2 -G255/0/0 -P -K -O >> %fig_out%
pscoast -R -J -W -P -K -O >> %fig_out%
call legend.bat
del _*
del .*
rem ***************************************************************************
rem * gawk -f awk_klegend.awk inp_legend2.txt > legend.bat
rem *     凡例描画用バッチファイルを作成するqwkスクリプトの実行
rem *     (入力ファイルはinp_legend2.txt，出力ファイルはlegend.bat)
rem * pscoast -R90/150/0/45 -JM12c -Ba15g15 -G200 -S255/255/255 -N1 -P -K > %fig_out%
rem *     -R90/150/0/45 作図範囲を東経90〜150度，北緯0〜45度に指定
rem *     -JM12c        メルカトル図法で幅12cmの図を作成
rem *     -Ba15g15      座標数値15度ピッチ，グリッド線15度ピッチで描画
rem *     -G200         陸域を灰色で着色
rem *     -S255/255/255 水域を白色で着色
rem *     -N1           陸域の国境を描く
rem *     -P            縦長用紙指定
rem *     -K            追加するコマンドがあることを示すオプション
rem * gawk "{print $6,$7}" inp_eqdat.txt | psxy -: -R -J -B -Sc0.1 -G255/255/255 -W1/255/0/0 -P -K -O >> %fig_out%
rem *     gawk "{print $6,$7}" inp_eqdat.txt
rem *             gawkにより入力ファイルinp_eqdat.txtの6列目と7列目を
rem *             psxyの入力データとする
rem *     -:      入力ファイルinp_eqdat.txtの6・7列は緯度・経度の順であるが，
rem *             GMTでの指定は標準が経度・緯度の順なので，これを弱点させる．
rem *             (x,y)の入力順序を(y,x)に変更する．
rem *     -Sc0.1 -G255/255/255 -W1/255/0/0
rem *             赤縁白抜き直径0.1cmの円を描画  
rem *     -P      縦長用紙指定
rem *     -K -O   中間コマンドでの約束ごと
rem * gawk "{print $6,$7}" inp_eqdat_M7.txt | psxy -: -R -J -B -Sc0.2 -G255/0/0 -P -K -O >> %fig_out%
rem *     gawk "{print $6,$7}" inp_eqdat_M7.txt
rem *             gawkにより入力ファイルinp_eqdat_M7.txtの6列目と7列目を
rem *             psxyの入力データとする
rem *             (この入力ファイルにはM7以上の自信データが格納されている)
rem *     -Sc0.2 -G255/0/0
rem *             直径0.2cmの赤丸を描画
rem * pscoast -R -J -W -P -K -O >> %fig_out%
rem *     プロットが多く，陸域の境界が見えなくなるので，陸域境界線を再描画．
rem * call legend.bat  awkで作成した凡例描画用バッチファイル実行
rem ***************************************************************************
rem

ダウンロードと画像確認

ファイル名	概要
awk_klegend.txt	凡例作成用awkスクリプト(拡張子はtxtにしています)
inp_legend2.txt	凡例作成用awkスクリプト実行用入力ファイル
legend.txt	凡例作成用awkスクリプト実行結果バッチファイル(拡張子はtxtにしています)
bat_wg_eqcoast.txt	制御用バッチファイル(拡張子はtxtにしています)
fig_eq_east.png	出力画像(ImageMagickでpngに変換)