如何使用Python绘制任意平面内的向量场？

Question

我有一个三维速度向量场，在一个numpy数组中（zlength，ylength，xlength，3）。“3”包含速度分量（u，v，w）。在

我可以很容易地在正交的x-y，x-z，y-z平面上绘制向量场

X, Y = np.meshgrid(xvalues, yvalues)

xyfieldfig = plt.figure()
xyfieldax = xyfieldfig.add_subplot(111)
Q1 = xyfieldax.quiver(X, Y, velocity_field[zslice,:,:,0], velocity_field[zslice,:,:,1])

但是，我希望能够观察任意平面内的速度场。在

我试图将速度场投影到一个平面上，方法是：

^{pr2}$

然而，这（当然）仍然留给我一个具有相同形状的3d numpy数组：（zlength，ylength，xlength，3）。投影的_场现在包含每个（x，y，z）位置的速度矢量，这些矢量位于每个局部（x，y，z）位置的平面内。在

如何将速度场投影到一个平面上？或者，我现在如何沿着一个平面绘制投影的_场？在

提前谢谢！在

Answer 1

你很接近。丹尼尔F的建议是对的，你只需要知道如何做插值。这是一个有效的例子

from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
import numpy as np
import scipy.interpolate  


def norm(v,axis=0):
    return np.sqrt(np.sum(v**2,axis=axis))

#Original velocity field
xpoints = np.arange(-.2, .21, 0.05)
ypoints = np.arange(-.2, .21, 0.05)
zpoints = np.arange(-.2, .21, 0.05)

x, y, z = np.meshgrid(xpoints,ypoints,zpoints,indexing='ij')

#Simple example
#(u,v,w) are the components of your velocity field
u = x
v = y
w = z 

#Setup a template for the projection plane. z-axis will be rotated to point
#along the plane normal
planex, planey, planez  =
    np.meshgrid(np.arange(-.2,.2001,.1), 
                 np.arange(-.2,.2001,.1), [0.1], 
                 indexing='ij')

planeNormal = np.array([0.1,0.4,.4])
planeNormal /= norm(planeNormal)

#pick an arbirtrary vector for projection x-axis
u0 = np.array([-(planeNormal[2] + planeNormal[1])/planeNormal[0], 1, 1])
u1 = -np.cross(planeNormal,u0) 
u0 /= norm(u0)
u1 /= norm(u1)

#rotation matrix
rotation = np.array([u0,u1,planeNormal]).T


#Rotate plane to get projection vertices
rotatedVertices = rotation.dot( np.array( [planex.flatten(), planey.flatten(), planez.flatten()]) ).T


#Now you can interpolate gridded vector field to rotated vertices
uprime = scipy.interpolate.interpn( (xpoints,ypoints,zpoints), u, rotatedVertices, bounds_error=False )
vprime = scipy.interpolate.interpn( (xpoints,ypoints,zpoints), v, rotatedVertices, bounds_error=False )
wprime = scipy.interpolate.interpn( (xpoints,ypoints,zpoints), w, rotatedVertices, bounds_error=False )

#Projections
cosineMagnitudes = planeNormal.dot( np.array([uprime,vprime,wprime]) )

uProjected = uprime - planeNormal[0]*cosineMagnitudes
vProjected = vprime - planeNormal[1]*cosineMagnitudes
wProjected = wprime - planeNormal[2]*cosineMagnitudes

如果您想增加一些tensordot操作，可以减少行数。另外，这个或一些相近的变体在meshgrid中没有indexing='ij'也可以工作。在

原始字段：

^{pr2}$

投影场：

fig = plt.figure()
ax = fig.gca(projection='3d')

ax.quiver(rotatedVertices[:,0], rotatedVertices[:,1], rotatedVertices[:,2], 
          uprime, vprime,wprime, length=0.5, color='blue', label='Interpolation only')
ax.quiver(rotatedVertices[:,0], rotatedVertices[:,1], rotatedVertices[:,2], 
          uProjected, vProjected, wProjected, length=0.5, color='red', label='Interpolation + Projection')


plt.legend()