【pyevtk】PythonでParaView向けのファイル入出力

この記事ではPythonからCAE(Computer Aided Engineering)分野でよく使われている可視化ツールParaView向けのファイルを入出力できるpyevtkについて解説します。

この記事を読むことであなたは

・ParaViewとは何か？
・vtkフォーマットとは何か?
・pyevtkによるvtkフォーマットの出力方法
について理解することができます。

ParaView

ParaViewとは?

　ParaView(https://www.paraview.org/)とは、科学技術計算の計算結果を可視化するためのオープンソースの可視化アプリケーションです。フリーで使えるにもかかわらず、非常に多くの可視化機能があるため、大学、企業問わず非常に多くのユーザーから人気を集めています。オープンソースのCAEソフトの可視化ツールとしてはほぼデファクトスタンダード化しているといえます。例えば、ParaViewを可視化ツールとして使ったオープンソースの流体解析ソフトウェアとしてはOpenFOAMが有名です。ParaViewは2000年ごろから、ロスアラモス国立研究所とKitware社によって開発プロジェクトが始まり、現在に至るとのことです。

　Pythonにも可視化のライブラリとしては、matplotlibやseaborn, plotlyなどがありますが、これらのライブラリは統計データや2次元までデータの可視化に向いているため、3次元データや科学技術計算の結果はこちらのParaViewを使う方が良いと思います。

ParaViewのインストール

　インストールは公式ホームページからダウンロードします。

ParaView向けのフォーマットvtk

　ParaViewでは基本的にvtk(Visualization Toolkit)というフォーマットのファイルを読み込みます。vtkファイルは後でpyevtkによる出力方法を解説しますので、自分で書く必要はありませんが次のようなフォーマットです。(Legacyフォーマットといわれます。)

# vtk DataFile Version 2.0
vector
ASCII
DATASET UNSTRUCTURED_GRID
POINTS 4 float
0 0 0
1 0 0
1 1 0
0 1 0
CELLS 1 5
4 0 1 2 3
CELL_TYPES 1
9
POINT_DATA 4
VECTORS point_vectors float
1 1 0
-1 1 0
-1 -1 0
1 -1 0
CELL_DATA 1
VECTORS cell_vectors float 
1 1 0

vtkは更にいくつかの種類に分けられます。よく使うフォーマットを以下に纏めます。これ以外にもありますが、これらを覚えておけば実用上は問題ないと思います。

ImageData(.vti): 構造イメージデータ。等間隔格子などはvtiを使うと便利です。
StructuredGrid(.vtr): 構造格子データ。こちらは必ずしも等間隔でなくても構いません。
UnstructuredGrid(.vtu): 粒子などの非構造のグリッドデータに使います。

これ以外のフォーマットやvtkに関する詳しい情報は公式のドキュメントをご覧ください。

【pyevtk】Pythonのvtkフォーマット出力ライブラリ

pyevtkとは / インストール方法

　pyevtk(Python Export VTK)はPythonのサードパーティーライブラリで、Numpyアレーを所定のvtkフォーマットで簡単に出力することができます。pyevtkは標準ではインストールされていないので、pipでインストールします。

pip install pyevtk

pyevtkの基本の使い方

　ここでは上で挙げたvti, vtr, vtuを例にサンプルプログラムを示します。次のプログラムはvti用の出力プログラムです。

from pyevtk.hl import imageToVTK
import numpy as np


x = np.linspace(0, 10, 50)
y = np.linspace(0, 10, 50)
z = np.linspace(0, 10, 50)
X, Y, Z = np.meshgrid(x, y, z, indexing="ij")

func = np.sin(X) * np.cos(Y) * Z

cell_data = {"data" : func}
imageToVTK("./image_res", pointData=cell_data)

　まず1行目でpyevtk.hlからimageToVTKをimportしています。ここでは、2行目以下で適当な3次元等間隔格子のNumpyアレーを発生させています。実際に、vtiを出力しているのはimageToVTKになります。この関数には“ファイルのパス(拡張子を除いたもの)”、”データ(等間隔格子)”の２つを与えます。データはどの関数でも共通ですが、辞書で与えます。ここではcell_data = {“data”: func}という風に指定していますね。

次にvtuのサンプルプログラムを示します。

from pyevtk.hl import pointsToVTK
import numpy as np

points = 100
x = np.random.rand(points)
y = np.random.rand(points)
z = np.random.rand(points)

pressure = np.random.rand(points)
temp = np.random.rand(points)
point_data = {"temp": temp, "pressure": pressure}
pointsToVTK("./point_res", x, y, z, data=point_data)

　流れは先ほどのvtiと同じです。vtiとの違いはデータが等間隔格子ではなくx, y, zの座標を持った点であるということです。ここでは乱数から100点ランダムに発生させています。vtuの場合は、pointsToVTKという関数を使います。pointsToVTKには“ファイルのパス”、”ポイントの座標”、”ポイントの持つ値”を渡します。

最後にvtrについてです。次のプログラムを見てください。

import numpy as np
from pyevtk.hl import gridToVTK

nx, ny, nz = 6, 6, 2 
lx, ly, lz = 1.0, 1.0, 1.0 
dx, dy, dz = lx/nx, ly/ny, lz/nz 
ncells = nx * ny * nz 
npoints = (nx + 1) * (ny + 1) * (nz + 1) 

x = np.linspace(0, 10, nx) + 0.0001 * np.random.randn(nx)
y = np.linspace(0, 10, ny) + 0.0001 * np.random.randn(ny)
z = np.linspace(0, 10, nz) + 0.0001 * np.random.randn(nz)

pressure = np.random.rand(ncells).reshape( (nx, ny, nz)) 
temp = np.random.rand(npoints).reshape( (nx + 1, ny + 1, nz + 1)) 
gridToVTK("./structured", x, y, z, cellData = {"pressure" : pressure}, pointData={"temperature": temp})

　vtrの出力にはgridToVTKを使います。ここでは等間隔格子に少しだけノイズを加えて、その格子に対して出力するようにしています。gridToVTKには“ファイル名”、”格子のx”、 “格子のy”、 “格子のz”、”データ”を与える必要があります。ここで一点だけ注意点なのはセルデータの場合、軸毎に２つの格子点で挟まれるセルに対して、物理量を定義するので、格子点-1のサイズのNumpyアレーになるということです。格子点の上に定義されるポイントデータは格子点のディメンションと同じになります。