GMT でヒストグラムを描く事例です.
プロットするデータは,気象庁のホームページに掲載されている,前橋市の過去 114 年間の年最大日雨量です.
並び替えを自動で行ってくれ,プロットデータも出力できるので便利です.
実行用バッチファイルと解説
rem **************************
rem ヒストグラム描画
rem **************************
set inp_org=inp_PRB_DM.txt
set out_text=out_DM_HISTO.txt
set wd=10
set xlabel="Maximum daily rainfall (mm/day)"
set fig_out=fig_DM_HISTO.eps
rem
minmax %inp_org% -H1 > %out_text%
pshistogram %inp_org% -W%wd% -IO -H1 >> %out_text%
rem ***************************************************************************
rem * minmax %inp_org% -H1 > %out_text% *
rem * ファイル%inp_org%の内容より,ファイル名,データ数,最小/最大を *
rem * ファイル%out_text%に新規出力 *
rem * -H1 により先頭1行読み飛ばし *
rem * minmax の出力形式は以下のとおり. *
rem * (ファイル名): N = (データ数) <最小値/最大値> *
rem * 事例== inp_PRB_DM.txt: N = 114 <43.2/357.4> *
rem * *
rem * pshistogram %inp_org% -W%wd% -IO -H1 >> %out_text% *
rem * ファイル%inp_org%の内容を値域幅%wd%でx値,y値(頻度)を *
rem * ファイル%out_text%に追加出力 *
rem * -W%wd% === 値域幅wd(頻度を集計するx軸幅) *
rem * -IO ===x値,y値(頻度)を出力 *
rem * -H1 ===先頭1行読み飛ばし *
rem * WARNING: 269046071 multiple modes found という警告がでるが *
rem * ファイル出力は行われる. *
rem * -IOオプションでxy値の出力は行われるが-Iオプションのみでは同様の *
rem * 警告が出て画面出力はするがファイル出力はしない *
rem * このため,最大・最小をファイルに落とすにはminmaxコマンドを *
rem * 使用する必要有. *
rem ***************************************************************************
rem
gawk "{if(NR==1){nd=$4;str0=$5}}END{sub(/,\"\",str0);sub(/>/,\"\",str0);split(str0,a,\"/\");print nd,a[1],a[2]}" %out_text% > _tempdat.txt
rem ***************************************************************************
rem * ファイル%out_text%の1行目(minmaxの出力)からデータ数,最小値,最大値を *
rem * 抜き出しファイル_tempdat.txtに書き出す *
rem ***************************************************************************
rem
gawk "{printf \"0.95 0.95 12 0 0 MR n=%%g\n\",$1}" _tempdat.txt > _txtdat.txt
gawk "{printf \"0.95 0.90 12 0 0 MR min=%%g\n\",$2}" _tempdat.txt >> _txtdat.txt
gawk "{printf \"0.95 0.85 12 0 0 MR max=%%g\n\",$3}" _tempdat.txt >> _txtdat.txt
rem ***************************************************************************
rem * ファイル_tempdat.txt(1行のファイル)からデータ数,最小値,最大値の描画用 *
rem * ファイル_txtdat.txtに書き出す *
rem ***************************************************************************
rem
gmtset ANOT_FONT_SIZE 14
gmtset LABEL_FONT_SIZE 14
gmtset TICK_LENGTH -0.2c
pshistogram %inp_org% -JX15/10 -B:%xlabel%:/:"Frequency":WSen -H1 -W%wd% -L0.5p -G230 -Z0 -P -K > %fig_out%
pstext _txtdat.txt -R0/1/0/1 -J -N -P -O -K >> %fig_out%
gawk "BEGIN{FS=\",\"}{if(NR==1)printf \"0.2 10.2 12 0 0 BL %%s\n\",$0}" %inp_org% | pstext -R -J -N -P -O >> %fig_out%
rem ***************************************************************************
rem * pshistogram *
rem * -JX15/10 === 描画領域幅15cm,高さ10cm *
rem * -B:"Maximum daily rainfall (mm/day)":/:"Frequency":WSen *
rem * === 軸ラベルと数値を描画する軸を指定 *
rem * -H1 === 先頭1行読み飛ばし *
rem * -W%wd% === 値域幅wd(頻度を集計するx軸幅) *
rem * -L0.5p === 指定した太さのペンで棒の輪郭を描く *
rem * [デフォルトでは輪郭なし] ここでは0.5ptを指定 *
rem * -G200 === 領域をグレースケール230で塗りつぶし *
rem * -Z0 === 縦軸の選択.デフォルト(指定無しor -Z0)は個数. *
rem * -Z1を指定すれば頻度% *
rem * -P === 縦長紙を指定 *
rem * -Rを指定していないので軸範囲は自動選定 *
rem * *
rem * pstext *
rem * -R0/1/0/1 === pshistogramでは軸範囲はGMTにまかせたので, *
rem * データ数・最小値・最大値描画時は軸範囲を再設定する *
rem * 最後に左上にコメントを追加 *
rem ***************************************************************************
set inp_org=
set out_text=
set wd=
set xlabel=
set fig_out=
del .gmt*
del _*
|
バッチファイルのダウンロードと画像確認