Scatter plotting 3D Numpy array using matplotlib

Question

I have a 3D numpy array of size (75, 150, 150) of numeric values which represents 75 layers, 150 grid cells in the x and 150 grid cells in the y directions. Is there a way to 3D plot this 3D array in x y, and z 3D map where the z dimension represents the layer of my model and x and y be the x and y dimensions? I have the following code but its giving me error: ValueError: shape mismatch: objects cannot be broadcast to a single shape.

array.shape = (75, 150, 150)

from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
n = 1687500

z=np.arange(76)
x=np.arange(151)
y=np.arange(151)

for i, j, k in x, y, z:
  y=array[k,i, :]
  x=array[k,:, j]
  colors = np.random.randint(0, 10, size=n)
  ax.scatter(z, x, y, c=colors, marker='o')
  plt.show()

Answer 1

You have an array of 75*150*150 (1687500) elements. In your first version of the scatter plot, you were using one column of 75 elements for Z and 2 columns of 150 values for X and Y. You are ignoring almost all information from the 3D array.

In the edited version, x,y and z each contain the full 3D array. All values will be the same.

Supposing the array contains 75 layers of grids of 150x150, you could plot it as follows. Note that it will be quite slow as the number of points is huge. Also, it will be hard to make sense of. np.meshgrid creates arrays of position values.

To filter out the zero values, you can replace them by np.nan. To color according to bins, you can use a BoundaryNorm.

from matplotlib import pyplot as plt
from matplotlib.colors import BoundaryNorm
import numpy as np

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

num_layers = 75
num_pnt = 150
z, x, y = np.meshgrid(np.arange(1, num_layers + 1), np.arange(num_pnt), np.arange(num_pnt), indexing='ij')

# create some random test data, suppose all values outside a cone are zero
array = np.random.rand(num_layers, num_pnt, num_pnt) ** 2
array[(x - num_pnt / 2) ** 2 + (y - num_pnt / 2) ** 2 > (num_layers - z) ** 2] = 0

array[array == 0] = np.nan  # replace zeros by NaN to make them invisible

bounds = [0, 0.0001, 0.001, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10, 1.00]
norm = BoundaryNorm(bounds, len(bounds) - 1)
cmap = plt.get_cmap('turbo', len(bounds) - 1)
scat = ax.scatter(x, y, z, c=array, marker='o', cmap=cmap, norm=norm)
cbar = plt.colorbar(scat, ax=ax, ticks=bounds, format='%.4f')
plt.show()

When there are many zeros, it helps to completely filter them away. You'll need to convert all arrays to 1D:

from matplotlib import pyplot as plt
from matplotlib.colors import BoundaryNorm
from matplotlib.ticker import FuncFormatter
import numpy as np

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

num_layers = 75
num_pnt = 150
z, x, y = np.meshgrid(np.arange(1, num_layers + 1), np.arange(num_pnt), np.arange(num_pnt), indexing='ij')

# create some random test data, suppose all values outside a cone are zero
array = np.random.rand(num_layers, num_pnt, num_pnt) ** 2
array[np.abs((x - num_pnt / 2) ** 2 + (y - num_pnt / 2) ** 2 - (num_layers - z) ** 2) > 5] = 0

# make the arrays 1D, so they are easier to filter
array = array.ravel()
filter = array != 0
x = x.ravel()[filter]
y = y.ravel()[filter]
z = z.ravel()[filter]
array = array[filter]

bounds = [0, 0.0001, 0.001, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10, 1.00]
norm = BoundaryNorm(bounds, len(bounds) - 1)
cmap = plt.get_cmap('turbo', len(bounds) - 1)
scat = ax.scatter(x, y, z, c=array, marker='o', cmap=cmap, norm=norm)
cbar = plt.colorbar(scat, ax=ax, ticks=bounds, format=FuncFormatter(lambda x, pos: f'{x * 100:3g} %'))
plt.show()