IDLによる解析・可視化

著者

松本洋介（千葉大学）、簑島敬（海洋開発研究機構）

本章では、 CANS+ での解析・可視化ツールであるIDL（Interactive Data Language）の基本的な使い方と、用意されているIDLプロシージャ及びアニメーションの作り方を説明します。IDLのバージョンが8.0以降では、グラフィックス環境が大きく変わり、可視化が容易かつ洗練されました。 CANS+ のデータ読み込み、それらの可視化機能、解析手法の紹介をします。さらなる詳細は IDLのドキュメントサイト を参照してください。

$CANSPLUS_DIR/idl内には、以下のようなIDLのプロシージャが用意されています。これらは各課題から共通して使えるものです。

課題に固有のプロシージャは各課題ディレクトリ内に用意されており、共通プロシージャを内部で利用しています。

環境変数の設定

環境変数$IDL_STARTUPに$CANSPLUS_DIR/idl/init.proを設定します。

bashの場合、

export IDL_STARTUP=$CANSPLUS_DIR/idl/init.pro

tcshの場合、

setenv IDL_STARTUP $CANSPLUS_DIR/idl/init.pro

$CANSPLUS_DIR/idl/init.pro内には、pathの設定、IDL内の環境の設定等が含まれており、IDL起動時に自動的に設定されます。各自の好みに合わせて修正してください。

データの読み込み

オリジナル CANS で使われていたバイナリデータフォーマットであるDAC形式データを読み込む共通プロシージャが用意されています。各MPIプロセスごとに出力されたデータをまとめて一つの配列に格納します。

dac_read.pro

Name
  dac_read
  DACデータを読み込む

Syntax
  dac_read,data,x,y,filename 　　; 2D データ
  dac_read,data,x,y,z,filename ; 3D データ

Arguments
  data
    記録するデータ配列

  x, y(, z)
    シミュレーションのセル幅から計算した、ｘ，ｙ(，ｚ)軸

  filename
    読み込むデータのファイル名。ワイルドカードが使え、パターンに
    マッチしたファイルを一度に読み込む。ただし、読みこむ全ての数
    値データ配列数は等しい必要がある。

Example1
  IDL> dac_read,bx,x,y,z,'010000_bx_rank=*.dac'

Example2
  IDL> dac_read,allvars,x,y,'010000_*_rank=*.dac'

Example3
  IDL> dac_read,bxall,x,y,z,'0[1-5]0?00_bx_rank=*.dac'

1次元プロット

1次元のプロットを作成する方法を以下に簡単に示します。

IDL> dac_read,data,x,y,z,'0011_*_rank=*.dac' ;; データ読み込み
IDL> pl=plot(x,data[*,6],xtitle='x',ytitle='$\rho$') ;; 軸のタイトルにLatexコマンドが使える！
_images/idl1dplot1.png

1次元プロットの例。左は修正前、右は後から修正を加えたもの。

これに修正を加える場合は、plのpropertyに値を代入します。

IDL> pl.symbol='circle' ;; ○の追加
IDL> pl.color='red'     ;; 線を赤に
IDL> pl1.xstyle=3       ;; x方向の描画範囲を少し広げる
IDL> pl1.ystyle=3       ;; y方向の描画範囲を少し広げる
IDL> pl.font_size=16    ;; フォントサイズを16に
IDL> pl.font_name='Times' ;; フォントをTimesに
IDL> leg=legend(target=pl,position=[0.85,0.85],/data,/auto_text_color)
IDL> leg.font_name = 'Times'
IDL> leg.font_size = '14'
IDL> leg.label = '$\rho$'

のようににして後からプロットの属性に変更を加えることができます。GUIでの操作も可能です。凡例を加えるには、legend関数を使います。

複数のプロットを並べて表示する場合はlayoutオプションを使います。この例では2x2のパネル状に並べて表示します。

IDL> pl1=plot(x,data[*,6],xstyle=1,ystyle=2,color='red',$
IDL>          layout=[2,2,1])
IDL> pl2=plot(x,data[*,5],xstyle=1,ystyle=2,color='blue',$
IDL>          layout=[2,2,2],/current)
IDL> pl3=plot(x,data[*,7],xstyle=1,ystyle=2,color='purple',$
IDL>          layout=[2,2,3],/current)
IDL> pl3=plot(x,data[*,8],xstyle=1,ystyle=2,color='green',$
IDL>          layout=[2,2,4],/current)
IDL> pl1.title='$\rho$'
IDL> pl2.title='$P$'
IDL> pl3.title='$V_x$'
IDL> pl4.title='$V_y$'
_images/idl1dplot2.png

上記の設定の結果

するとこのように表示されます。

その他の応用は、サンプルプロシージャ md_shktb/batch.proをご覧ください。また、plot, legend関数に関する情報は

IDL> ?plot
IDL> ?legend

とすることで、ブラウザで見ることができます。

2次元可視化

カラーマップ図

物理量の２次元プロファイルを表示する方法を以下に示します。

IDL> dac_read,data,x,y,z,'0013*.dac'
IDL> dx = abs(x[1]-x[0]) ;; 軸の位置の修正のため
IDL> dy = abs(y[1]-y[0]) ;; 軸の位置の修正のため
IDL> img = image(data[*,*,5],x-dx/2.,y-dy/2.,axis_style=2,$
IDL              xtickdir=1,ytickdir=1,rgb_table=33)
IDL> cb  = colorbar(target=img,orientation=1,tickdir=1,$
IDL                 position=[0.85,0.1,0.9,0.9])

imageで２次元配列データを画像化し、colorbarでカラーバーを指定の場所に配置します。すると

_images/idl2dplot1.png

2次元カラー図の例。左が修正前、右が修正後。

これでは画像が小さいので拡大、フォント名・サイズを変更、軸にタイトルを付け加え見栄え良くします。

IDL> img.font_name='Times'
IDL> img.font_size='16'
IDL> img.xtitle='$X (\lambda)$' ;;CAN USE LATEX COMMANDS FOR TITLE
IDL> img.ytitle='$Y (\lambda)$'
IDL> img.xstyle=1
IDL> img.ystyle=1
IDL> img.scale,2.0,2.0 ;; ENLARGE IF IMAGE IS SMALL
IDL> cb.font_name='Times'
IDL> cb.font_size='14'

これらの属性の修正は、関数の引数として最初から渡すことも可能です。その他の応用は、サンプルプロシージャ md_OTvortex/batch.pro, md_kh/batch.pro をご覧ください。また、image、colorbar関数に関する情報は

IDL> ?image
IDL> ?colorbar

とすることで、ブラウザで見ることができます。

ベクトル図

ベクトル場を矢印で表現することも可能です。例えば、上記で作成した2次元 カラー図 に速度場ベクトルを付け加える場合は、

IDL> vct = vector(data[*,*,7],data[*,*,8],x-dx/2.,y-dy/2.,/overplot)

のようにします（単体で表示する場合はoverplotオプションを外す）。

_images/idl2dplot2.png

2次元カラー図+速度場ベクトルの例。左が修正前、右が修正後。

しかし、そのままだと全ての情報を矢印で描くことになってしまうので、サンプリングして描画するように修正する。

IDL> vct.x_subsample = 10
IDL> vct.y_subsample = 7
IDL> vct.transparency = 50 ;; ついでに矢印を透過（50%）

その他の応用は、サンプルプロシージャ md_kh/batch.proをご覧ください。また、vector関数に関する情報は

IDL> ?vector

とすることで、ブラウザで見ることができます。

等高線図

磁力線を描く（ベクトルポテンシャルの等高線）ためには、等高線を図示するのが便利です。例えば、上記で作成した2次元 カラー図 に質量密度の等高線を描き加える（この例では圧力分布とほぼ同じ）場合は、

IDL> cnt = contour(data[*,*,6],x-dx/2.,y-dy/2.,/overplot,c_color=0)

のようにして行います。

_images/idl2dplot3.png

2次元カラー図+等高線の例。左が修正前、右が修正後。

等高線ラベルを削除、等高線の数を増やす場合は、

IDL> cnt.c_label_show=0
IDL> cnt.n_levels=30

のようにして後から修正ができます（関数に引数として渡すことも可能）。contour関数に関する情報は

IDL> ?contour

とすることで、ブラウザで見ることができます。

ベクトル場の流線

ベクトルデータから流線を描くには、以下の共通プロシージャが利用できます。

field_lines_2d.pro

Name
  field_lines_2d
    ２次元のベクトル量から流線（磁力線）を描く

Syntax
  fl = field_lines_2d(ux, uy, [npos=npos,] [nsteps=nsteps,] [len=len,]
                      [r0=r0,] [dir=dir])

Return Value
  fl
    描く流（磁力）線の２次元座標情報。位置を各方向のグリッド数で規格化。

Arguments
  ux, uy
    描画する２次元ベクトル量

Keywords
  npos
    線の数。デフォルトは100。

  nsteps
    トレースするステップ数。デフォルトは10000。

  len
    各方向１ステップあたりの長さ（グリッド単位）を指定する2次元配列

  r0
    各線のスタート位置（x、y）。与えられなければ、ランダムに決められる。
    r0=r0[2,npos]である必要がある。

  dir
    トレースする方向。+1ならばベクトルの向き。-1ならば反対向き。デフォルトは+1。

Example
  IDL> dac_read,data,x,y,z,'0011*_rank=*.dac'
  IDL> dx = abs(x[1]-x[0])
  IDL> dy = abs(y[1]-y[0])
  IDL> lx = max(x)-min(x)+dx
  IDL> ly = max(y)-min(y)+dy
  IDL> img = image(data[*,*,5],x-dx/2.,y-dy/2.,axis_style=2,$
  IDL>             xtickdir=1,ytickdir=1,rgb_table=33)
  IDL> npos=10 ;; Number of starting points
  IDL> fl=field_lines_2d(data[*,*,0],data[*,*,1],npos=npos)
  IDL> .run
  - for k=0,npos-1 do begin
  - pl = plot(fl[0,*,*,k]*lx,fl[1,*,*,k]*ly,/overplot,color='black')
  - endfor
  - end
  IDL>

フーリエ変換

現象の発展に伴ってどのようなモードが成長しているのかを把握するためには、フーリエ変換を行います。 1次元データと2次元データをフーリエ変換する共通プロシージャを用意しています。

fftf.pro

Name
  fftf
    1次元データをフーリエ変換し、振動数や波数に対して昇順にソートする。

Syntax
  f1d = fftf(data, time, direction, freq=freq)

Return Value
  f1d
    入力された1次元データ data のフーリエ変換。振動数や波数に対して昇順にソートされている。
    関数の平方の総和（積分）が、そのフーリエ変換の平方の総和（積分）と等しくなるように規格
    化されている（パーセバルの定理）。

  freq
    独立変数 time に対応する振動数や波数が出力される。2$\pi$が掛かっていないことに注意。

Arguments
  data
    1次元データ。

  time
    時刻や座標といった独立変数。

  direction
    フーリエ変換の方向。-1ならフーリエ変換、+1なら逆変換。指定しない場合は-1。

Example
  IDL> data = randomu(systime(/sec), 512)
  IDL> time = findgen(512)
  IDL> f1d = fftf(data, time, freq=freq)
  IDL> plot, freq, abs(f1d), /yl

fft2d.pro

Name
  fft2d
    2次元データをフーリエ変換し、振動数や波数に対して昇順にソートする。

Syntax
  f2d = fft2d(data, x, t, direction, wnum=wnum, freq=freq)

Return Value
  f2d
    入力された2次元データ data のフーリエ変換。振動数や波数に対して昇順にソートされている。
    関数の平方の総和（積分）が、そのフーリエ変換の平方の総和（積分）と等しくなるように規格
    化されている（パーセバルの定理）。

  wnum
    独立変数 x に対応する振動数や波数。2$\pi$が掛かっていないことに注意。

  freq
    独立変数 t に対応する振動数や波数。2$\pi$が掛かっていないことに注意。

Arguments
  data
    2次元データ。

  x, t
    時刻や座標といった独立変数。

  direction
    フーリエ変換の方向。-1ならフーリエ変換、+1なら逆変換。指定しない場合は-1。

Example
  IDL> data = randomu(systime(/sec), 128, 128)
  IDL> x = findgen(128)
  IDL> t = findgen(128)
  IDL> f2d = fft2d(data, x, t, wnum=wnum, freq=freq)
  IDL> plot_clcnt, alog10(abs(f2d)), xax=wnum, yax=freq

図の出力

上記 2次元カラー図 等で作成した図をファイルに保存するには、

IDL> img = image(data[*,*,5],x-dx/2.,y-dy/2.,axis_style=2,$
IDL>             xtickdir=1,ytickdir=1,rgb_table=33)
IDL> img.save,'result.png' ;; 画像で保存
IDL> img.save,'result.eps' ;; ポストスクリプトで保存

のようにして、saveメソッドを利用して保存します（plotでも同様）。ファイル形式は拡張子で自動判定します。

サンプルプロシージャ

IDLによる可視化のサンプルプロシージャ batch.pro が課題ディレクトリ "md_shktb/" , "md_OTvortex/" , "md_kh/" にありますので、適宜修正してお使いください。 batch.proは

IDL> .run batch

で実行し、書かれている内容をバッチ処理します。使用例は各課題（ 衝撃波管問題 、 Orszag-Tang渦問題 、 ケルビン・ヘルムホルツ不安定 ）のページをご覧ください。

アニメーションの作り方

以下の様に、IDL上で結果をpngファイルとして保存します。

IDL> dac_read,ro,x,y,z,'*_ro_rank=*.dac'
IDL> info=size(ro,/dimension)
IDL> nt = info[2]
IDL>.run
- for i=0,nt-1 do begin
-   if(i eq 0)then begin
-     img = image(ro[*,*,i],x-dx/2.,y-dy/2.,axis_style=2,$
-                 xtickdir=1,ytickdir=1,rgb_table=33)
-   endif else begin
-     img.setdata,ro[*,*,i]
-   endelse
-   img.save,'result'+strcompress(string(i,format='(i3.3)'),/remove)+'.png'
- endfor
- end

上記のようにすると、作業ディレクトリには"result005.png"のように連番でファイル名が付けられた画像ファイルが作成されます。

結果をpngファイルに出力後、LinuxにインストールされているImagemagick (convert)やffmpegを使って、結果の図をアニメーション化します。まず、convertを使った方法を説明します。

$ convert result???.png result.gif

のようにすると、gifアニメーションが作成されます。convertの際の詳しいオプションは"man convert"に委ねます。

次に、ffmpegを使ったアニメーション作成例を示します。

$ ffmpeg -i result%03d.png -sameq -vcodec mjpeg result.avi

のように、gifアニメーション以外のフォーマット（この例では、元の画像と同じ画質で、エンコードにmjpeg、格納形式としてaviを指定している）の動画ファイルが作成可能です。更なるオプションは"man ffmpeg"に委ねます。

上記の様に、IDL上での結果を連番で名前を付けた画像ファイルに保存すれば、Windows, Mac-OSでのその他のツールも利用可能だと思います。