TL;DR np.meshgrid を呼び出すときにオプション indexing='ij' を指定するのが良いです。
状況
Pyvista で、RectilinearGrid や StructuredGrid を初期化するときに、奇妙な振る舞いを見つけました。
ちなみに Windows11, Python 3.9, Pyvista 0.37, Numpy 1.14 で確認しました。
以下のコードを実行したときの RectilinearGrid と StructuredGrid の Dimension を見てみます。
import numpy as np
import pyvista as pv
xrng = np.arange(-10, 10, 2, dtype=np.float32)
yrng = np.arange(-10, 10, 5, dtype=np.float32)
zrng = np.arange(-10, 10, 1, dtype=np.float32)
grid_rl = pv.RectilinearGrid(xrng, yrng, zrng)
x, y, z = np.meshgrid(xrng, yrng, zrng)
grid_st = pv.StructuredGrid(x, y, z)
このコードは、https://docs.pyvista.org/version/stable/api/core/_autosummary/pyvista.StructuredGrid.html を参照して書きました。
RectilinearGrid grid_rl は 10 x 4 x 20 になります。
In [2]: grid_rl
Out[2]: RectilinearGrid (0x1e7a21991c0)
N Cells: 513
N Points: 800
X Bounds: -1.000e+01, 8.000e+00
Y Bounds: -1.000e+01, 5.000e+00
Z Bounds: -1.000e+01, 9.000e+00
Dimensions: 10, 4, 20
N Arrays: 0
StructuredGrid grid_st は 4 x 10 x 20 です。
In [3]: grid_st
Out[3]: StructuredGrid (0x1e7a21b1760)
N Cells: 513
N Points: 800
X Bounds: -1.000e+01, 8.000e+00
Y Bounds: -1.000e+01, 5.000e+00
Z Bounds: -1.000e+01, 9.000e+00
Dimensions: 4, 10, 20
N Arrays: 0
cell_data の先頭に値を入れて plot で様子を見ます。
grid_st.cell_data["idx"] = np.zeros(grid_st.n_cells)
grid_st.cell_data["idx"][:3] = np.arange(3)
grid_st.plot(show_edges=True)
最初に y 軸方向に変化するようになっていることがわかります。この状況は「わかりやすくない」ので、初に x 軸方向に変化するようにしたくなります。
RectilinearGrid の場合は、x軸、y軸、z軸の順に変化します。こちらのほうが「わかりやすく」て使い勝手の良いふるまいです。
どうしてこうなるのか
少し調べてみると、これは pyvista の問題ではなくて、numpy の mesh_grid の挙動に沿ったものです。https://numpy.org/doc/stable/reference/generated/numpy.meshgrid.html Returns: のところを見ると
For vectors x1, x2,…, xn with lengths
Ni=len(xi), returns(N1, N2, N3,..., Nn)shaped arrays if indexing=’ij’ or(N2, N1, N3,..., Nn)shaped arrays if indexing=’xy’ with the elements of xi repeated to fill the matrix along the first dimension for x1, the second for x2 and so on.
と書いてあります。つまり、デフォルトの indexing = 'xy' だと、(N2, N1, N3) の配列を返すと書いてあります。それで、y-x-z の並び順になるようです。
対応法
x, y, z = np.meshgrid(xrng, yrng, zrng, indexing='ij')
grid_st = pv.StructuredGrid(x, y, z)
とします。
In [18]: grid_st
Out[18]:
StructuredGrid (0x1edb3031340)
N Cells: 513
N Points: 800
X Bounds: -1.000e+01, 8.000e+00
Y Bounds: -1.000e+01, 5.000e+00
Z Bounds: -1.000e+01, 9.000e+00
Dimensions: 10, 4, 20
N Arrays: 0
StructuredGrid grid_st も 10 x 4 x 20 となり、x 軸方向が最初の変化方向になっていることが確かめられました。
今後
Pyvista の StructuredGrid の説明ページを修正するようプルリクを書きました。