Wednesday, 11 August 2021

PYTHON Plotting 3D Vector Field

The Data

I have a vector field, which is 0 in all components except for the z component. I just have the data for one slice of this field. My goal is to show this slice in a 3D plot.

The slice:

enter image description here

import numpy as np
import matplotlib.pyplot as plt

# DATA FIELD (1,N,M)
dz = np.array([[[0.24884899, 0.24884899, 0.24884899, 0.24884899, 0.24884899,
         0.24884899],
        [0.248849  , 0.248849  , 0.248849  , 0.248849  , 0.248849  ,
         0.248849  ],
        [0.24885767, 0.24885513, 0.24885108, 0.24885113, 0.2488552 ,
         0.24885767],
        [0.2451304 , 0.24563262, 0.24642831, 0.24641793, 0.24561579,
         0.2451304 ],
        [0.0764377 , 0.12581053, 0.09866768, 0.10043774, 0.12461962,
         0.0764377 ],
        [0.03382106, 0.03394624, 0.03414449, 0.03414171, 0.03394174,
         0.03382106]]])

dx = np.zeros(np.shape(dz))
dy = np.zeros(np.shape(dz))

# DATA POINTS (N,)
X = np.array([0. , 0.2, 0.4, 0.6, 0.8, 1. ])
Y = X 
Z = np.array([-500., -360., -220.,  -80.,   60.,  200.])

To create the plot:

# COMPUTE LENGTH OF VECTORS
length = np.sqrt(dx[0]**2+ dz[0]**2)

# PLOT 2D Cross-Section of vector field
fig = plt.figure(dpi=300)
Q = plt.quiver(X, Z, dx[0], dz[0], length, units='xy' ,angles='xy', scale=0.005,  
       pivot = "tail", headaxislength = 5, headlength = 5, cmap='jet')
fig.colorbar(Q)
plt.gca().invert_yaxis()
plt.ylabel("Z")
plt.xlabel("X")
plt.title("2D Cross-Section")

The Problem:

So far so good, but when I try to plot it in 3D, it seems that I am doing something wrong, since

  1. the arrows are not visible anymore
  2. the plot seems to be rotated
# Create 3D Quiver Plot with color gradient
# Source: https://stackoverflow.com/questions/65254887/how-to-plot-with-matplotlib-a-3d-quiver-plot-with-color-gradient-for-length-giv

from mpl_toolkits.mplot3d import Axes3D # <--- This is important for 3d plotting 

def plot_3d_quiver(x, y, z, u, v, w):
    # COMPUTE LENGTH OF VECTOR -> MAGNITUDE
    c = np.sqrt(np.abs(v)**2 + np.abs(u)**2 + np.abs(w)**2)

    c = (c.ravel() - c.min())/c.ptp()
    # Repeat for each body line and two head lines
    c = np.concatenate((c, np.repeat(c, 2)))
    # Colormap
    c = plt.cm.jet(c)

    fig = plt.figure(dpi = 300)
    ax = fig.gca(projection = '3d')
    ax.quiver(x, y, z, u, v, w, colors = c, length = 1, arrow_length_ratio = 0.5)
    plt.gca().invert_zaxis()
    plt.show()


# Create Mesh !
xi, yi, zi = np.meshgrid(X, Y, Z, indexing = 'xy')

# slice mesh
skip_every = 1
skip_slice = 6
skip3D = (slice(None, None, skip_slice), slice(None, None, skip_every), slice(None, None, skip_every))

plot_3d_quiver(xi[skip3D], yi[skip3D], zi[skip3D], dx, dy, dz)

enter image description here

I am wondering if you could help me? Thanks a lot!

EDIT:

As pointed out by @Marc the arrows are there in the 3D plot, but apparently, the arrow size is relative to the size of the axis. See more here: https://github.com/matplotlib/matplotlib/issues/11746

So far there seems to be no other fix than to scale all axis to roughly the same values. So, if -as suggested by @Marc- I divide the z-axis by 1000, I can see the arrows again:

plot_3d_quiver(xi[skip3D], yi[skip3D], zi[skip3D]/1000, dx, dy, dz)

Then we get:

enter image description here

which clearly shows that the orientation is somehow off...

EDIT 2:

What I would like to achieve:

The data that I have is just for one cross-section (one plane) of the 3D space (-> dimensions (1,N,M)), hence if I would show it in 3D, it should look something like this:

enter image description here

FULL DATA

In the end, I would like to correctly plot the full vector field. You can download it here: https://filebin.net/csjvs7xde5lfbwek

# %% 
import pickle
import numpy as np
import matplotlib.pyplot as plt

# Import Full Data
with open('full_data.pickle', 'rb') as handle:
    full_data = pickle.load(handle)

# Axis
X = np.linspace(0,1,101)
Y = np.linspace(0,1,10)
Z = np.linspace(-500,200,101)

for cross_section in full_data["cross_sections"].keys():
    
    # extract field components in x, y, and z 
    dx,dy,dz = storage_dict["cross_sections"][cross_section]
    
    # Make them numpy imediatley
    dx = np.array(dx)
    dy = np.array(dy)
    dz = np.array(dz)
    
    # ------------------------
    # plot cross-sections z-x
    # -----------------------
    
    length=np.sqrt(sx**2+sz**2) 

    fig = plt.figure(dpi=300)
    # If U and V are 2-D arrays and X and Y are 1-D, 
    # and if len(X) and len(Y) match the column and row dimensions of U, 
    # then X and Y will be expanded with numpy.meshgrid().
    # To make the head a triangle, make headaxislength the same as headlength.
    Q = plt.quiver(X[skip1D], Z[skip1D], sx[skip2D], sz[skip2D],\
               length[skip2D], units='xy' ,angles='xy', scale=0.02,  
               pivot = "tail", headaxislength = 5, headlength = 5, cmap='jet')
    fig.colorbar(Q)
    plt.title("Cross-Section: " + str(cross_section))
    # Invert y axis to have same represntation than MC GRATING
    plt.gca().invert_yaxis()
    
    plt.pause(0.01)
  
# ------------------------
# plot full 3D Vector Field ?? 
# -----------------------


from PYTHON Plotting 3D Vector Field

No comments:

Post a Comment