Plot density field depending on x, y, z coordinates with python

Question

I have extracted nodal fields from my geodynamic simulation, but I'm stuck on visualizing it. Here I take rock density rho [kg/m3] as an example. The data are represented on the nodal points of a 3D rectangular grid, with a resolution of 10 * 3 * 10 km. The lengths of my coordinates are X * Y * Z = 213 * 69 * 133 nodes, where X and Z are horizontal and Y is the vertical coordinate. The dimensions of rho are 69 * 213 * 133 (so y, x, z)

Ideally, I would like to draw some sort of box in (x,z,y) space and assign a color to it according to the value of rho. However, for now I'm fine with creating a color-coded scatter plot at the nodes. I've looked in many places, e.g. the matplotlib.mpl_toolkits.Axes3D documentation, and here, and almost all of them say something like

img = ax.scatter(x, z, y, c=rho, cmap=plt.hot())

Or versions using ax.plot_surface, or ax.plot_trisurf, and those explicitly require X, Y, Z to be of the same length. ax.scatter also only seems to work when x, y, z and rho have the same length, which seems unnecessary to me... However I try, I get errors like (on calling the scatter command):

  File "/usr/lib/python3/dist-packages/mpl_toolkits/mplot3d/axes3d.py", line 2354, in scatter
    xs, ys, zs = np.broadcast_arrays(
  File "<__array_function__ internals>", line 5, in broadcast_arrays
  File "/usr/lib/python3/dist-packages/numpy/lib/stride_tricks.py", line 264, in broadcast_arrays
    shape = _broadcast_shape(*args)
  File "/usr/lib/python3/dist-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

I looked up what that means, and supposedly it means that the dimensions don't match. But I checked in a debugger, and the dimensions are correct.

I tried flattening rho (same result) and I tried looping over each node and plotting each point separately (sample down below)


    rho, x, y, z = read_data('/path/to/my/data')
    fig = plt.figure()
    ax = fig.add_subplot(1,1,1, projection='3d')
    for i in range(len(x)):
        for j in range(len(z)):
            for k in range(len(y)):
                ax.scatter(x[i], z[j], y[k], c=ronew[k, i, j], cmap='hot')
    cbar = plt.colorbar(orientation='horizontal')
    cbar.set_label('Density [kg/m3]')
    plt.xlabel("X [km]")
    plt.ylabel("Z [km]")
    plt.zlabel("Depth [km]")
    plt.show()

I imagine I'm not the only one who needs to make these "4D" figures. Therefore, any help would be very much appreciated! To make helping easier, here's a sample data that produces the above copied error:

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import pyplot as plt

x = np.arange(0, 10)
y = np.arange(0, 6)
z = np.arange(0, 20)

rho = np.random.randint(1000, 3500, (10, 6, 20))

fig = plt.figure()
print(x.shape, y.shape, z.shape, rho.shape)
ax = fig.add_subplot(projection='3d')
ax.scatter(x, y, z, c=rho, cmap='hot')
plt.show()

You could use `np.meshgrid` to create a 3D grid out of the 3 1d arrays: `x3d, y3d, z3d = np.meshgrid(x, y, z)` and then `ax.scatter(x3d, y3d, z3d, c=rho, cmap='hot')` — JohanC, Aug 12 '21 at 15:25
Thanks for the suggestion. I tried it out, and then I get an error saying that `rho` is an invalid RGBA argument. — L. van Agtmaal, Aug 12 '21 at 15:31
Did you forget `c=` in `ax.scatter(..., c=rho, ...)`? Are you using the latest matplotlib 3.4.2? Do you also have that problem with the test code at the end of your question? — JohanC, Aug 12 '21 at 15:37
Note that `x`, `y`, and `z` do not explicitly define the coordinate system — they are actually the coordinates of *each point* in the field. That's why they have to be the same length. You might want to check out [`xarray`](http://xarray.pydata.org/en/stable/), which was made for data like yours. — Matt Hall, Aug 12 '21 at 23:06
@JohanC I did not forget the c, but apparently I am using matplotlib 3.1.4. I tried with the test code and that gives the same error. I'll see if I can update matplotlib — L. van Agtmaal, Aug 13 '21 at 07:09
@kwinkunks I tried using 3 for loops for x y z coordinates, then something like `ax.scatter(y[i], x[j], z[k], c=rho[i,j,k], cmap='hot')` but that produced the broadcast error as well — L. van Agtmaal, Aug 13 '21 at 07:12

score 0 · Answer 1 · answered Aug 13 '21 at 07:40

Okay, so thanks to @JohanC 's comment I found out that the 4D scatter is not possible with matplotlib version 3.1.4 (the reason is beyond me). After upgrading matplotlib to 3.4.3 I had no problems after calling xx, yy, zz = np.meshgrid(x, y, z) and plotting like ax.scatter(xx, zz, -yy, c=rho, cmap='hot')

I produced the following figure from the test code:

And this one from the real data file

Plot density field depending on x, y, z coordinates with python

1 Answers1