個人メモ

松葉が寄せ集めた個人メモです.不定期に更新します. 現在は分類が不十分なので,徐々に整理していきます. 基本的には動作確認をしたものだけを掲載しています.

GPhys/GGraph関連

基本設定

作図するときの毎回の設定(お作法).そろそろ描画準備のこの一連の設定はメソッドにしたほうがよいかもしれないと思うようになった.

require 'numru/ggraph'
require '/home/matsuba/define' # 自作メソッド集
include NumRu

if __FILE__ == $0 then
  # ファイル名の設定
  psname = "Fig/" + $0.split('.')[0]
  DCL.swcset('fname', psname)

  # 描画準備
  DCLExt.sw_set_params 'rlwfact'=>999 # PDF 出力のときに画面に似せる
  iws = (ARGV[0]||(puts 'Workstation ID(I)?'; DCL.sgpwsn; gets)).to_i
  DCL.gropn(iws)
  GGraph.margin_info(nil) # 実行スクリプトを描画の隅に描く
  DCL.sldiv('y', width=1, height=1)
  DCLExt.sg_set_params 'lcntl'=>true, 'lfull'=>true, 'isub'=>96,\
  'lfprop'=>true, 'ifont'=>2
  DCLExt.uz_set_params 'indext1'=>3,'indext2'=>5,'indexl1'=>5,\
  'indexl2'=>5,'inner'=>-1
  DCLExt.ud_set_params 'lmsg'=>false, 'indxmj'=>5, 'indxmn'=>3
  DCL.uzfact(0.7)

  # ----- Main Routine ----- #
  DCL.grcls
end 

GGraph.margin_info は DCL.sldiv の前に呼ぶこと. 順番が逆だと文字列の下半分が切れてしまう(バグ?).

GPhysオブジェクトの軸の型を変更する

VArray#to_f メソッドと Axis#set_pos メソッドを組み合わせる.

vgp = gp.axis('pressure').pos.to_f # 座標値を実数化
gp.axis('pressure').set_pos(vgp)

bash関連

デフォルトのエディタを変える

less や crontab などで起動するデフォルトのエディタは vim であるが, emacs ユーザーにとっては不便極まりないのでなんとか変更したい. ~/.bashrc で以下のように記入する.

export EDITOR=/usr/bin/"emacs -nw"

普段 no window のオプションを使用しているので -nw を付けている.

pdfファイルの総ページ数を取得する (2016-07-19 追記)

pdfinfo と gawk を利用する.対象とする PDF を hoge.pdf として以下のように使う(動作確認済).

pdfinfo hoge.pdf 2>/dev/null | gawk '/Pages/ {print $2}'

なお,pdfinfo, gawk が入っていなければ,debian であれば sudo apt-get install で取得できる.

sudo apt-get install gawk
sudo apt-get install poppler-utils

poppler-utils のなかに pdfinfo が入っている(apt-cache search pdfinfo の検索結果でこのパッケージに行き着く).

psファイルの総ページ数を取得する

pdf の総ページ数の取得は割とよく知られているが, ps ファイルの総ページ数の取得はあまり書かれていないのでメモ.

$ grep -c 'showpage' hoge.ps

-c はカウントするオプション.これと ps ファイルを分割する psselect とを利用すると, 複数ページある ps ファイルを 1 枚 1 枚に分割することができる.

psfile=hoge
pages=`grep -c 'showpage' ${psfile}.ps`
for i in `seq 1 ${pages}`; do
   num=`printf "%03d" ${i}` 
   outfile=${psfile}_${num} # hoge_001, hoge_002 のように命名
   psselect ${i} ${psfile}.ps ${outfile}.ps   
done

2016-05-11 追記

どうやらこれが使えるのはDCL が吐き出したpsファイルのみの可能性が出てきた.DCL 6.x からはPDF出力がデフォルトになったのでこの方法が使えなくなった.pdf2ps で変換かけて作成したpsファイルにはshowpageの情報が書かれないようだ.詳しく調査します.

ディレクトリ名一覧を取得する (2016-07-19 追記)

awk を使うことでできる.

ls -l | awk '$1 ~ /d/ {print $9 }'

データ置き場

気象庁ベストトラックデータ(txt形式)

以下のURLから取得することができる(気象庁HP英語版).

http://www.jma.go.jp/jma/jma-eng/jma-center/rsmc-hp-pub-eg/besttrack.html

Ruby関連

月末の日にちを取得する

Dateクラスを用いる.日にちの引数に-1を指定する.

require 'date'
Date.new(2007, 2, -1).day
 #=> 28

TeX 関連

プリアンブルの設定(2016-07-27 現在)

いろいろと試行錯誤した結果,以下に落ち着いている.\partd, \dpartd はよく使うので常時自作メソッドとして呼んでおく.一番最後のおまじないは,作成したPDFにしおりをつける方法.

\documentclass[11pt, a4j]{jarticle}
\newcommand{\partd}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\dpartd}[2]{\dfrac{\partial #1}{\partial #2}}
\usepackage[dvipdfmx,hiresbb]{graphicx}
\usepackage[dvipdfmx]{color}
\usepackage[dvipdfmx,colorlinks,linkcolor=blue]{hyperref}
\usepackage{amsmath, amssymb}
\usepackage{bm}

\setlength{\abovecaptionskip}{5pt}
\setlength{\belowcaptionskip}{5pt}
\makeatletter
\renewcommand{\theequation}{%
  \thesection.\arabic{equation}}
  \@addtoreset{equation}{section}
  \makeatother

\AtBeginDvi{\special{pdf:tounicode EUC-UCS2}}

未分類

2次元座標変換

本家のサイトには座標変換に関する説明(補間,座標変換)があるが,2次元の座標変換に関しては説明が少ないというかほとんどなされていないので,補足できる範囲でしてみる.なお,理解が足らず間違いが含まれているかもしれないので,ご容赦ください.

たとえば,極座標のデータ gp(r, theta) を直交座標系にする(よくあるのは直交座標を極座標系にするだろうが,その説明は後日)ことを考える.手順は主には以下である.

gx(r, theta) を作成するときにちょっと工夫がいる.gr(r)*gtheta(theta) のような計算では GPhys は次元を拡張してはくれないためである.

grad = gp.axis('radius').to_gphys
gtheta = gp.axis('theta').to_gphys.convert_units('rad')

zero = gp*0  # ここがミソ.次元拡張のために (r, theta) の次元を持ったゼロ配列を用意する
gx = grad*(zero + gtheta.cos) # gx の属性は grad のものが引き継がれる
gx.long_name = 'x-coordinate'

gy = grad*(zero + gtheta.sin)
gy.long_name = 'y-coordinate'

カラーレベルの作成

GGraph.tone は内部で UEPACK を呼んでおり,DCL.uegtlb(min, max, dx) でカラーレベルを作成している(さらに内部で DCL.uegtla を呼んでいる).切りのよい目盛り打ちになるのは嬉しいが,端の最大・最小値が中途半端になるのがどうもいやだった.そこで,中途半端な値をそぎ落として,無限大の三角にするカラーレベル作成メソッドを作った.それが以下である.最小値,最大値は GGraph.tone では内部で計算するが,ここでは陽に与える必要がある.なお,このメソッド作成にあたって DCL.uegtla のソースを参考にしている.

def clrbar_rmiss_uegtla(min, max, nlev)
  # 実数化しておく
  min = min.to_f
  max = max.to_f
  # UEPACK に関する設定
  DCL.ueiset('nlev', nlev)
  lbound = DCL.uelget('lbound')
  DCL.uegtla(min, max, 0)
  # ここから uegtla で設定されたパターンを解析する
  num = DCL.ueqntl
  dz = DCL.rgnge((max-min)/nlev)
  zmin = DCL.irle(min/dz)*dz
  zmax = DCL.irge(max/dz)*dz

  tlev1_first_flag = false
  tlev2_last_flag = false
  for i in [1, num]
    tlev1 = dz*((zmin + (i-1)*dz)/dz).round
    tlev2 = tlev1 + dz
    if lbound && i == 1 then
      tlev1_first_flag = true if tlev1 == min
    elsif lbound && i == num then
      tlev2_last_flag = true if tlev2 == max
    end
  end

  tlev1 = Array.new(num) # for tlev1
  tlev2 = Array.new(num) # for tlev2
  for i in 0..num-1
    tlev1[i], tlev2[i], ipat = DCL.ueqtlv(i+1)
  end
  # レベルを設定していく
  levels = Array.new
  tlev1[1..-1].each{|val| levels.push(val)}
  levels.unshift(tlev1[0]) if tlev1_first_flag
  levels.push(tlev2[-1]) if tlev2_last_flag
  # 欠損値
  rmiss = DCL.glrget('rmiss')
  levels.unshift(rmiss)
  levels.push(rmiss)
  p levels
  return levels
end

(追記 2017-03-28)

min と max が一致していると,dz=0 となるのでゼロ割りが発生して落ちるみたいです(RuntimeError: [UEGTLA] XMIN SHOULD BE LEAST THAN XMAX.).dcdvlopに投稿された内容とたぶん同質のものだと思うのでご注意.なお,gphys-1.5.1 から?同様の機能を持つメソッドが導入されていると推測されるので,そっちを使えばよいのでは?

sldiv でフレームを分割したときに中途半端なフレームから図を描く.

DCL.sldiv('y', width=3, height=4) のようにフレームを分割したとき, 最初の枠からではなく途中から描きたいときのやり方.改ページを何度か呼べば実現できそう.改ページに相当するのは DCL.grfrm であるので,これを必要な回数だけ呼べばよい.以下は,自分が必要とするやり方(1ページに00LST~11LST, 12LST~23LST の絵が並ぶようにしたいとき).当然 DCL.gropn は呼んでいるものとする.

hour = (Define::basetime2date(tarray[0] + 3600*7)).slice(8..9).to_i
i = hour%(width*height)
i.times{|num| DCL.grfrm}

以下,通常どおりの描画.

DCL でパラパラアニメーションにする

次の設定でOK.ただし,swpack の設定なので gropn の前に呼ぶこと.

DCLExt.sw_set_params 'lalt'=>true, 'lwait'=>false 

NArray と Numeric との2項演算の仕様

Numeric に合わせて NArray を upcast する演算は行なわれない (upcast とは,整数と浮動小数点との演算なら,整数を浮動小数点に変換した上で演算すること).なので,GPhys オブジェクトの座標値や物理量データが整数値だと, Numeric との演算のときに不具合が生じる可能性があるので注意する.

書式付き文字列

Rubyの場合

day = sprintf("%02d", i)

bashの場合

day=`printf "%02d" ${i}` # printf で書式付きで出力した結果を代入する

高品位のフォントを使用する

DCLでフォントを変更するときのおまじない

DCLExt.sg_set_params 'lfprop'=>true, 'ifont'=>2

ベクトル描画自由自在

自分なりに調べてたどりついた結果ではあるが,間違っているかも.

# unit vector の大きさを陽にスケーリングするときの雛型
# viewport を陽に与える (ユニットベクトルを描く位置を調整するため)

vxunit = 0.05; vyunit = 0.05
uxunit = 30.0; uyunit = 30.0

xfact1 = vxunit/uxunit
yfact1 = vyunit/uyunit

hash_vector = {
  'lnrmal'=>false, 'xfact1'=>xfact1, 'yfact1'=>yfact1, 
  'lmsg'=>false, 'lunit'=>false, 
  'uxunit'=>uxunit, 'uyunit'=>uyunit, 'lumsg'=>false, 
  'vxuloc'=>viewport[1]+0.05, 'vyuloc'=>viewport[2]+0.02,
  'index'=>15, 'iuindx'=>15
}

DCLExt.ug_set_params hash_vector
DCL.ugsut('x', sprintf("%#.1f",uxunit))
DCL.ugsut('y', sprintf("%#.1f",uyunit))

以上の設定をしたうえで,GGraph.vector を呼ぶ. 各種変数の意味については ugpack の内部変数のマニュアルを参照のこと (UGPACKの内部変数). そのときに,ユニットベクトルを描かないといけないので 'unit_vect'=>true とし,上記の設定を反映させるために 'flow_vect'=>false とする.具体的にはたとえば以下のように記述する.

opt_vector = {
   'xintv'=>10, 'yintv'=>10, 'unit_vect'=>true, 'flow_vect'=>false,
   'transpose'=>false, 'title'=>'', 'annotate'=>false
 }

 GGraph.set_fig 'viewport'=>viewport
 GGraph.vector uc, vc, true, opt_vector

(注): DCL のマニュアルのなかにスケーリングの説明がなされているページを見たことがあるのだが,どこで見たのか思い出せないでいる.思い出せばここにリンクを記そうと思う.ぜひとも情報を寄せていただきたい.

東西・南北風速成分から風向を計算する

東西風速を u,南北風速を v とする.u>0 で西風,v>0 で南風を表わすものとする. このとき風が吹き込んでくる角度は,北からを0度として時計回りに数えていくものとするとき (すなわち東風は90度,南風は180度,西風は270度),

(u.atan2(v) + Math::PI)*180.0/Math::PI

で求められる.GPhys#atan2(*args) の使い方に注意.

カラーパターンを逆順に指定する

GGraph.tone のオプションで 'clr_min' => 87,'clr_max' => 13 のように,通常とは逆に指定する.内部では上限と下限を設定したのちにトーンレベルに応じて線形にパターンを指定しているだけなのでこのような応用が可能.

鉛直座標変換

モデル面データを気圧面データに変換するなど.GPhys#interpolate メソッドを用いて線形内挿する.JMA-NHM の場合,対数圧力座標を対象として線形内挿している.自作した変換メソッドを以下に示す(あとで掲載する予定).

欠損値も含めて NetCDF に書き出すとき

データの型を NArrayMiss にするだけではダメ.書き出したい GPhys オブジェクトの属性に missing_value を付け足す.

rmiss = DCL.glrget('rmiss')
gp.put_att('missing_value', NArray[rmiss])

GPhys#interpolate の返り値である GPhys オブジェクトにはこの属性は付かない仕様になっているので,NetCDF に書き出す前に手動で設定する必要がある.

GGraph のメソッドで無理矢理対数軸のラベルを変えるとき

1000, 500, 200, 100,・・・のような表記にしたいときの対処法. 本家の方によると現状ではおまかせにせずに,

GGraph.set_axes 'yside'=>'r'

で左側の y 軸を書かないようにしてから,あとから自分で

DCL.ulylog 

を呼んで対数軸を書くとよいとのコメントを頂いている. しかしここでは敢えて GGraph の描画メソッドに頼りきって描画を試みようと思う.

内部では DCL.usdaxs で軸を描いている.'itr'=>2 とすると y 軸が対数軸になるが,このとき y 軸は DCL.ulylog で描いている.ラベルの設定は uspack の 'itypey' を参照していることがソース解読によりわかったので,

DCL.usiset('itypey', 3) 

とすれば望みの絵がかけるはずだった.

ところが,GGraph の描画メソッドを呼ぶと,内部で DCL.grfrm/DCL.grfig を呼んでいるので,uspack の設定は一度リセットされてしまうらしい.そこで

GGraph.next_fig 'no_new_fig'=>true

として内部で呼ばないようにし,代わりに自分で DCL.grfrm を呼ぶ (新たなフレームで描く場合).そのあとで DCL.usiset('itypey', 3) を呼ぶ.こうするとリセットされることなく設定が反映されるので,ラベルが 1000, 500, 200, 100 に置き換わる.まとめると,

GGraph.next_fig 'no_new_fig'=>true
DCL.grfrm
DCL.usiset('itypey', 3)
GGraph....

GGraph メソッドで描画するとき uzinit でリセットされてしまうパラメータの設定の仕方

上の対数軸を描くやり方の一般形.GGraph.tone などは内部で DCL.grfrm/grfig を呼んでいる.さらに内部で DCL.uzinit を呼んでいるため,該当するパラメータはリセットされてしまう.リセットされないように反映させるには,GGraph.set_fig 'no_new_fig'=>true とする.そのあとで自分で DCL.grfrm/grfig を呼び,パラメータを設定する.そうして GGraph の描画メソッドを呼べばよい.以下,GGraph の描画メソッドを呼んだときに内部で呼ばれる GGraph.fig の該当するソースを以下に掲載する.

if opts['no_new_fig']
  # do nothing
elsif opts['new_frame']
  DCL.grfrm
else
  if opts['keep_axis_offset'] # SAVE OFFSET VALUES
    offset_org = (%w(XT XB YL YR)).map{|zs| DCL::uzpget("ROFF#{zs}")}
  end
  DCL.grfig
end

経度線・緯度線を描かないようにする

地図投影を設定しているとき.umpack のパラメータで設定する

DCLExt.um_set_params 'lgridmj'=>false, 'lgridmn'=>false

ffmpegのインストール

debian パッケージなら以下のコマンドで一発.ただし事前に/etc/sources.list を編集する必要があるかもしれない.一度 apt-cache search でパッケージの有無を確認すること.

$ sudo apt-get install ffmpeg

DCL の出力から動画を作成する

未整理.DCL ユーティリティが使えなくなったのでとりあえずの代替案.

$ ruby anim_850vor.rb #=> dcl.pdf を出力
$ convert -quality 100 -rotate 90 dcl.pdf dcl_%04d.png # 連番でファイルを生成
$ ffmpeg -r 20 -i dcl_%04d.png -qscale 0 out.wmv

この方法だと画像の余分な縁は切り取れないし,画質も悪くなる(将来的に改善する).また,オプションの順番も重要で,-r の場所が違っていたりすると設定がうまく反映されない場合がある.

PDFで出力するときの線の太さを変える

DCL 6.x から導入された swpack の制御パラメータ.PDF で出力するときの線幅のファクター.デフォルト値は 0.2 である.999. は特殊な設定で画面に似せるようだ.わたしは画面表示に似せたい方なので,gropn の前に以下のように宣言している.

DCL.swpset('rlwfact', 999)

この情報は以下を参照した(2016-05-13 現在).

ビッグエンディアンで処理するときのコンパイラのオプション

big endian で書かれたデータを little endian で読もうとしてドツボに嵌った経験からメモしておく.

gfortran の場合

$ gfortran -fconvert=big-endian mfhm2gr.f90

ifort の場合

$ ifort -convert big_endian mfhm2gr.f90