3次元配列を座標変換してmatplotlibで3D化したい

Question

実現したいこと

物理シミュレーションして出てきた結果の配列をmatplotlibで３次元のベクトル図に可視化したいです。
###前提
シミュレーションソフト(CANS+)が出した物理量の結果を時間と位置座標[t,x,y,z]に対応した4次元配列で出力で打ち出すことはできています。このデータは位置座標が円柱座標になっているので3次元化するには、x'=xcos(y) y'=ysin(y)と座標変換した上で出力する必要があります。
ここでおきた問題は

• 座標変換がうまくいかない
• 元々のデータを直交座標として3Dにしようとしてもうまく行かない。

解決方法を探しています

発生している問題・エラーメッセージ

エラーメッセージ(座標変換しない場合)
File "main.py", line 28, in <module>
    three_pict.cm()
  File "/Users/Username/Desktop/picture/three_pict.py", line 28, in cm
    ax.quiver(x1[::4], y1[::4], z1[::4], data1[ts,::4,:,::4], data2[ts,::4,:,::4], data3[ts,::4,:,::4], length=0.1, normalize=True)
  File "/usr/local/lib/python3.7/site-packages/mpl_toolkits/mplot3d/axes3d.py", line 2677, in quiver
    bcast = np.broadcast_arrays(*input_args, *masks)
  File "<__array_function__ internals>", line 6, in broadcast_arrays
  File "/usr/local/lib/python3.7/site-packages/numpy/lib/stride_tricks.py", line 264, in broadcast_arrays
    shape = _broadcast_shape(*args)
  File "/usr/local/lib/python3.7/site-packages/numpy/lib/stride_tricks.py", line 191, in _broadcast_shape
    b = np.broadcast(*args[:32])
ValueError: shape mismatch: objects cannot be broadcast to a single shape

  File "main.py", line 28, in <module>
    three_pict.cm()
  File "/Users/username/Desktop/picture/three_pict.py", line 30, in cm
    x_car=x1*math.cos(y1*math.pi/16)
TypeError: only size-1 arrays can be converted to Python scalars

該当のソースコード

def cm():
    print('物理量を入力(v,b):')
    br=input()
    print('読み取るディレクトリを入力')
    df=input()
    dac_name1=str(df)+'/000?_'+str(br)+'x_rank=*.dac'
    dac_name2=str(df)+'/000?_'+str(br)+'y_rank=*.dac'
    dac_name3=str(df)+'/000?_'+str(br)+'z_rank=*.dac'
    data1, x1, y1, z1 = dac_read(dac_name1,dimension=3)
    data2, x2, y2, z2 = dac_read(dac_name2,dimension=3)
    data3, x3, y3, z3 = dac_read(dac_name3,dimension=3)
    #データを4次元配列に変換
    l=data1.shape[0]
    print('総タイムステップは0から'+ str(l-1)+'です')
    print('現像したいタイムナンバーは？(全部の場合はallと入力):')
    ts=int(input())
    fig=plt.figure()
    ax = fig.gca(projection='3d')
    x_car=x1*math.cos(y1*math.pi/16)
    y_car=x1*math.sin(y1*math.pi/16)
    z_car=z1
    ax.quiver(x_car[::4], y_car[::4], z_car[::4], data1[ts,::4,:,::4], data2[ts,::4,:,::4], data3[ts,::4,:,::4], length=0.1, normalize=True)
　　#座標変換しないときはx_car→x1 y_car→y1　z_car→z1で行った 
    ax.plot_surface(x_car, y_car, z,color='g',alpha=0.5)
    fig.colorbar(p)
    plt.show()

試したこと

x_carとy_carを先に配列として定義してみましたがダメでした。

補足情報（FW/ツールのバージョンなど）

python3.7.5

Comments

[WathMorks]

x1.shape?
y1.shape?

[Kimi.S]

x1.shapeは64
y1.shapeは16です。

[WathMorks]

サイズが異なるため（要素ごとの）積
x_car=x1*math.cos(y1*math.pi/16)
y_car=x1*math.sin(y1*math.pi/16)
は定義できません。

[Kimi.S]

なるほど。x1をr変更し、y1をthetaにして、meshgridでサイズを整えて定義しました。
radius_matrix, theta_matrix = np.meshgrid(r,theta)
x = radius_matrix * np.cos(theta_matrix)
y = radius_matrix * np.sin(theta_matrix)
fig = plt.figure()
ax = fig.gca(projection='3d')
u =data1[ts,:,:,]*np.cos(theta_matrix)-data2[ts,:,:,]*np.sin(theta_matrix)
v =data1[ts,:,:,]*np.sin(theta_matrix)+data2[ts,:,:,]*np.cos(theta_matrix)
w =data3[ts,:,:,:]
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
ax.quiver(x, y, z, u, v, w, length=0.1, normalize=True)
fig.show()

その結果
Traceback (most recent call last):
File "three_pict2.py", line 41, in <module>
ax.quiver(x, y, z, u, v, w, length=0.1, normalize=True)
File "/usr/local/lib/python3.7/site-packages/mpl_toolkits/mplot3d/axes3d.py", line 2677, in quiver
bcast = np.broadcast_arrays(*input_args, *masks)
File "<__array_function__ internals>", line 6, in broadcast_arrays
File "/usr/local/lib/python3.7/site-packages/numpy/lib/stride_tricks.py", line 264, in broadcast_arrays
shape = _broadcast_shape(*args)
File "/usr/local/lib/python3.7/site-packages/numpy/lib/stride_tricks.py", line 191, in _broadcast_shape
b = np.broadcast(*args[:32])
ValueError: shape mismatch: objects cannot be broadcast to a single shape
というエラーが出ました。
これはzのサイズが異なるのでuvwが定義できなかったということでしょうか

[WathMorks]

data1, x1, y1, z1 = dac_read(dac_name1,dimension=3)
ここで読み込んだx1,y1について
「x1.shapeとy1.shapeは異なる」と想定されていますか？

Answer 1

dataの座標配列x1,y1,z1をmeshgridで3次元配列にしてから直交座標に変換したらうまく行きました。

theta_matrix, z_matrix, radius_matrix = np.meshgrid(y1, z1, x1)
x = radius_matrix * np.cos(theta_matrix)
y = radius_matrix * np.sin(theta_matrix)
z = z_matrix
print(x.shape)
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
ax.quiver(x[::10,:,::5], y[::10,:,::5], z[::10,:,::5], data[::,::,::], length=1, normalize=True)
plt.show()