matplotlib (python) による解析・可視化¶

著者松本洋介（千葉大学）

はじめに¶

matplotlib（とnumpy）の基本的な使い方及びアニメーションの作り方を説明します。matplotlibはpythonプログラム言語をベースとした解析・可視化ライブラリで、IDLと比べてライセンスフリーで使えるため、手頃な解析環境として人気があります。

jupyter notebookは、pythonプログラミング及びmatplotlibによる可視化などのチュートリアル／ハンズオンのための便利なツールとして知られています。本notebookはjupyter notebookのセットアップができている環境下で動作します。

以下ではCANS+ のデータ読み込み、それらの可視化機能、解析手法の紹介をします。

本解説は$CANSPLUS_DIR/notebook内にmatplotlib_hands-on.ipynbとして用意されています。jupyter notebookを使って読んで下さい。

以下で使われるプロシージャの一部は$CANSPLUS_DIR/python内に用意されています。これらは各課題から共通して使えるものです。

python環境の準備¶

インストール¶

anaconda3をインストールするのがお手軽です。
https://www.anaconda.com/download/ よりバイナリファイルをダウンロードして、圧縮ファイルをインストール先に展開してください。
本サマーセミナーでは、すでにanaconda3がインストールされています。

環境設定¶

anaconda3をインストールした先を以下では ~/anaconda3 とします。

bashの場合、

export PATH=~/anaconda3:$PATH

cshの場合、

setenv PATH ~/anaconda3:${PATH}

とすることで、anaconda3内のpythonを優先的に使うように設定します。これにより、OSに別途インストールされているpython環境と分けて使用します。以上の設定は.bashrcや.tcshrcに書いておくとよいです。

[1]:

import sys
print(sys.version)

3.6.6 |Anaconda custom (64-bit)| (default, Jun 28 2018, 17:14:51)
[GCC 7.2.0]

のように出力されたら、インストールは正常に完了しています。試しにnumpyをインポートして利用してみましょう。

[2]:

import numpy as np
data = np.array([1000,100.0,1])

[3]:

print(data)

[1000.  100.    1.]

[4]:

print(data.dtype)

float64

jupyter notebookの起動¶

matplotlibを使った可視化・解析方法の例は$CANS_PLUSDIR/notebook/内のmatplotlib_hands-on.ipynbにあります。本notebookの利用方法は、

$ jupyter-notebook $CANSPLUS_DIR/notebook/matplotlib_hanbs-on.ipynb

で起動します。ブラウザが立ち上がり、ブラウザのページの中でコマンドをインタラクティブに実行していきます。

セル内の改行は Enter キー、セル内の命令の実行は Shift + Enter キーで行います。

データの読み込み¶

テキストファイルデータの読み込み¶

[5]:

import numpy as np
import glob
list = glob.glob('*_vyi.dat')
data = np.loadtxt(list[0])

IDLバイナリデータ（IDL save file）の読み込み¶

[6]:

import scipy.io as sio
data_set= sio.readsav('sample.dat')
data = data_set['data']
xax = data_set['x']
yax = data_set['y']

DACバイナリデータの読み込み¶

CANS+はCANSから使われる独自のDACバイナリ形式を採用しています。CANS+で出力したデータの読み込みルーチンの使い方を説明します。

bashの場合:

$ export PYTHONPATH=$CANSPLUS_DIR/python

cshの場合:

$ setenv PYTHONPATH $CANSPLUS_DIR/python

のように、pythonコマンドラインから読み込めるpathの設定を行います。その後

$ python

として、

from dac_read import dac_read
#./data/にデータがある場合
#2D計算の場合。10スナップショット目の全ての物理変数を読み込む例
data, x, y = dac_read('./data/0010_*_rank=*.dac',dimension=2)

#3D計算の場合。質量密度のみ全ての時間ステップを読み込む例
data, x, y, z = dac_read('./data/0???_ro_rank=*.dac',dimension=3)

の様にして計算結果を読み込むことができます。

注意¶

多次元配列の場合、pythonでは、行・列の順に格納されます。これは、それぞれの次元がY軸、X軸に相当していて、IDLやFotranと異なる（C言語と同じ）ことに注意が必要です。

matplotlibを使った可視化¶

matplotlibの使い方参考サイト¶

以下では、解析・可視化について基本的なところしか紹介しません。細かい調整等が必要な場合は下記のサイトが参考になると思いますので、適宜参照してください。

モジュールのインポート¶

まずはmatplotlibのモジュールをインポートします。おまじないだと思ってください。

[7]:

import matplotlib.pyplot as plt

描画ウィンドウの表示方法¶

notebook内では、描画内容はセル内に表示されていますが、pythonのコマンドライン上では、最後に

plt.show()

もしくは

plt.pause(0.1)

のようにして、新しく立ち上がるウィンドウ内に描画します。

前者と後者の方法の違いは、前者はウィンドウにフォーカスされて、ウィンドウを閉じない限りコマンドラインに戻れないのに対して、後者ではウィンドウを立ち上げたままコマンドラインに戻ることができます。

ipythonを使うとjupyter notebookと同じように使えて、

$ ipython --pylab

のように–pylabオプションをつけてipythonを起動すれば、描画時に自動的にウィンドウが立ち上がります。

1次元プロット¶

線プロットの作り方です。１次元plotの詳しい使い方は?plt.plotより確認できます。

[8]:

xpos1 = np.argmin(np.abs(xax-np.max(xax)/2))
plt.plot(yax,data[7,:,xpos1])

[8]:

[<matplotlib.lines.Line2D at 0x7fcc364debe0>]

軸にタイトルを入れたりして見栄えを変えます。

[9]:

#デフォルトフォントサイズの設定
#個別に指定しない限りは以降この設定が反映される
plt.rcParams['font.size'] = 18
#軸のタイトル、タイトルの追加
plt.plot(yax,data[7,:,xpos1])
plt.xlabel('$y$')
plt.ylabel('$V_x$')
plt.title('K-H instability',fontsize=24)

[9]:

Text(0.5,1,'K-H instability')

別のデータのプロットを重ね描きして、レジェンドを追加します。

[10]:

plt.plot(yax,data[7,:,xpos1],'ro-',markersize=5,label='$V_x$ at $x=L_x/2$')
xpos2 = np.argmin(np.abs(xax-np.max(xax)/8.))
plt.plot(yax,data[7,:,xpos2],'b+-',markersize=5,label='$V_x$ at $x=L_x/8$')
plt.title('K-H instability',fontsize=24)
plt.xlabel('$y$')
plt.ylabel('$V_x$')
plt.legend(loc='upper right',fontsize=14)

[10]:

<matplotlib.legend.Legend at 0x7fcc363565f8>

2次元マップ¶

pcolormeshを使った方法¶

[11]:

#x軸、y軸データから2次元データを作る
xax2d, yax2d = np.meshgrid(xax,yax)
plt.pcolormesh(xax2d,yax2d,data[5,:,:])

[11]:

<matplotlib.collections.QuadMesh at 0x7fcc362c8390>

軸タイトルやデータに基づくアスペクト比の変更、カラーバーを追加するなどして見栄えを整えます。

[12]:

plt.pcolormesh(xax2d,yax2d,data[5,:,:],cmap='bwr')
plt.xlabel('$X (\lambda)$')
plt.ylabel('$Y (\lambda)$')
##図のアスペクト比を軸データを元に決める
plt.gca().set_aspect('equal', adjustable='box')
plt.colorbar()

[12]:

<matplotlib.colorbar.Colorbar at 0x7fcc36260668>

カラーマップは https://matplotlib.org/examples/color/colormaps_reference.html を参考に指定してください。

ベクトル場の表示¶

ベクトル矢印による表現 (quiver）¶

上記で作成した2次元カラー図に速度場ベクトルを付け加える場合を例に示します。

[13]:

plt.pcolormesh(xax2d,yax2d,data[5,:,:],cmap='jet')
plt.xlabel('$X (\lambda)$')
plt.ylabel('$Y (\lambda)$')
plt.gca().set_aspect('equal', adjustable='box')
plt.colorbar()
plt.quiver(xax2d,yax2d,data[7,:,:],data[8,:,:])

[13]:

<matplotlib.quiver.Quiver at 0x7fcc361fb240>

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

[14]:

plt.pcolormesh(xax2d,yax2d,data[5,:,:],cmap='jet')
plt.xlabel('$X (\lambda)$')
plt.ylabel('$Y (\lambda)$')
plt.gca().set_aspect('equal', adjustable='box')
plt.colorbar()
#配列要素を5つ飛ばしで描画（::5）
plt.quiver(xax2d[::5,::5],yax2d[::5,::5],data[7,::5,::5],data[8,::5,::5],scale=np.max(xax),width=0.0075)

[14]:

<matplotlib.quiver.Quiver at 0x7fcc34da98d0>

流線による表現（streamplot）¶

矢印だとうまく表現できない場合もあるので、別の方法を紹介します。

[15]:

plt.pcolormesh(xax2d,yax2d,data[5,:,:],cmap='jet')
plt.xlabel('$X (\lambda)$')
plt.ylabel('$Y (\lambda)$')
plt.gca().set_aspect('equal', adjustable='box')
plt.colorbar()
#配列要素を3つ飛ばしで描画（::3）
plt.streamplot(xax2d[::3,::3],yax2d[::3,::3],data[7,::3,::3],data[8,::3,::3],color=np.sqrt(data[7,::3,::3]**2+data[8,::3,::3]**2),cmap='binary_r')
#ベクトル場の大きさで流線の色をつける
plt.colorbar()

[15]:

<matplotlib.colorbar.Colorbar at 0x7fcc34517278>

等高線¶

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

[16]:

plt.pcolormesh(xax2d,yax2d,data[5,:,:],cmap='jet')
plt.xlabel('$X (\lambda)$',fontsize=18)
plt.ylabel('$Y (\lambda)$',fontsize=18)
plt.tick_params(labelsize=16)
plt.gca().set_aspect('equal', adjustable='box')
plt.colorbar()
plt.contour(xax2d,yax2d,data[5,:,:])

[16]:

<matplotlib.contour.QuadContourSet at 0x7fcc344add68>

等高線の数を増やす場合は、

[17]:

plt.pcolormesh(xax2d,yax2d,data[5,:,:],cmap='jet')
plt.xlabel('$X (\lambda)$')
plt.ylabel('$Y (\lambda)$')
plt.gca().set_aspect('equal', adjustable='box')
plt.colorbar()
plt.contour(xax2d,yax2d,data[5,:,:],15,colors='black')

[17]:

<matplotlib.contour.QuadContourSet at 0x7fcc344631d0>

結果のファイルへの出力¶

上記2次元マップで作成した図をファイルに保存するには、それぞれのクラスに共通したsaveメソッドを使ってファイルに保存します。ファイル名の拡張子からファイルフォーマットが自動的に決まります。

[18]:

plt.pcolormesh(xax2d,yax2d,data[5,:,:],cmap='jet')
plt.xlabel('$X (\lambda)$')
plt.ylabel('$Y (\lambda)$')
plt.gca().set_aspect('equal', adjustable='box')
plt.colorbar()
plt.contour(xax2d,yax2d,data[5,:,:],15,colors='black')
#描画内容を保存 (解像度をDPI=150で指定)
plt.savefig('result_matplotlib.png',bbox_inches='tight',dpi=150)
plt.savefig('result_matplotlib.eps',bbox_inches='tight')
plt.close()

アニメーションの作り方¶

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

[19]:

data_set = sio.readsav('sample_ro-time.dat')
ro = data_set['ro']
nt = ro.shape[0]
for l in range(0,nt):
    plt.pcolormesh(xax2d,yax2d,ro[l,:,:],cmap='jet')
    plt.xlabel('$X (\lambda)$')
    plt.ylabel('$Y (\lambda)$')
    plt.gca().set_aspect('equal', adjustable='box')
    plt.colorbar()
    plt.savefig('result_matplotlib'+str(l).zfill(3)+'.png',bbox_inches='tight',dpi=150)
    plt.close()

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

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

$ convert result???.png result.gif

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

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

$ ffmpeg -r 10 -i result_matplotlib%03d.png -qscale 0 -vcodec mpeg4 result.mp4

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

[ ]: