How can I plot a numpy array of x, y, z in 3D surface plot?

Question

I have reviewed both of these threads, but am still struggling to make a 3D surface plot from a numpy array of x, y, z coordinates.

My array looks like this:

>>> points
array([[ 322697.1875    , 3663966.5       ,  -30000.        ],
       [ 325054.34375   , 3663966.5       ,  -30000.        ],
       [ 325054.34375   , 3665679.5       ,  -30000.        ],
       [ 322697.1875    , 3665679.5       ,  -30000.        ],
       [ 322697.1875    , 3663966.5       ,  -27703.12304688],
       [ 325054.34375   , 3663966.5       ,  -27703.15429688],
       [ 325054.34375   , 3665679.5       ,  -27703.70703125],
       [ 322697.1875    , 3665679.5       ,  -27703.67382812]])

ax.plot_surface accepts x, y, z points so I convert the above array into separate pieces below:

x = points[:, 0]
y = points[:, 1]
z = points[:, 2]

I then put it into a meshgrid for passing into ax.plot_surface():

import numpy as np

X, Y, Z = np.meshgrid(x, y, z)

And then try to plot:

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(16,10))
ax = plt.axes(projection = '3d')
ax.plot_surface(X, Y, Z, alpha=0.5)
plt.show()

When I run this I receive an error: rows, cols = Z.shape ValueError: too many values to unpack (expected 2).

I'm not sure where to go with this now, I don't expect the answer but a push in the correct direction would be great.

I would like the output to be similar in appearance to this but with my data:

UPDATE: If I do not include z in the meshgrid, but only x and y, I get this output when I run ax.plot_surface(X, Y, z, alpha=0.5):

This is really close, but I want all the sides to be filled in. Only one is showing as filled in. I've added the point coordinates to show the boundaries. I feel like it has something to do with the meshgrid that I'm creating. Here is the output of X, Y:

>>> X, Y = np.meshgrid(x, y)
(array([[322697.1875 , 325054.34375, 325054.34375, 322697.1875 ,
        322697.1875 , 325054.34375, 325054.34375, 322697.1875 ],
       [322697.1875 , 325054.34375, 325054.34375, 322697.1875 ,
        322697.1875 , 325054.34375, 325054.34375, 322697.1875 ],
       [322697.1875 , 325054.34375, 325054.34375, 322697.1875 ,
        322697.1875 , 325054.34375, 325054.34375, 322697.1875 ],
       [322697.1875 , 325054.34375, 325054.34375, 322697.1875 ,
        322697.1875 , 325054.34375, 325054.34375, 322697.1875 ],
       [322697.1875 , 325054.34375, 325054.34375, 322697.1875 ,
        322697.1875 , 325054.34375, 325054.34375, 322697.1875 ],
       [322697.1875 , 325054.34375, 325054.34375, 322697.1875 ,
        322697.1875 , 325054.34375, 325054.34375, 322697.1875 ],
       [322697.1875 , 325054.34375, 325054.34375, 322697.1875 ,
        322697.1875 , 325054.34375, 325054.34375, 322697.1875 ],
       [322697.1875 , 325054.34375, 325054.34375, 322697.1875 ,
        322697.1875 , 325054.34375, 325054.34375, 322697.1875 ]]), array([[3663966.5, 3663966.5, 3663966.5, 3663966.5, 3663966.5, 3663966.5,
        3663966.5, 3663966.5],
       [3663966.5, 3663966.5, 3663966.5, 3663966.5, 3663966.5, 3663966.5,
        3663966.5, 3663966.5],
       [3665679.5, 3665679.5, 3665679.5, 3665679.5, 3665679.5, 3665679.5,
        3665679.5, 3665679.5],
       [3665679.5, 3665679.5, 3665679.5, 3665679.5, 3665679.5, 3665679.5,
        3665679.5, 3665679.5],
       [3663966.5, 3663966.5, 3663966.5, 3663966.5, 3663966.5, 3663966.5,
        3663966.5, 3663966.5],
       [3663966.5, 3663966.5, 3663966.5, 3663966.5, 3663966.5, 3663966.5,
        3663966.5, 3663966.5],
       [3665679.5, 3665679.5, 3665679.5, 3665679.5, 3665679.5, 3665679.5,
        3665679.5, 3665679.5],
       [3665679.5, 3665679.5, 3665679.5, 3665679.5, 3665679.5, 3665679.5,
        3665679.5, 3665679.5]]))

If I just take x, y unique values I get an error thrown:

x = np.unique(x)
y = np.unique(y)

>>> x
array([322697.1875 , 325054.34375])
>>> y
array([3663966.5, 3665679.5])

X, Y = np.meshgrid(x, y)
>>> X, Y
(array([[322697.1875 , 325054.34375],
       [322697.1875 , 325054.34375]]), array([[3663966.5, 3663966.5],
       [3665679.5, 3665679.5]]))

>>> ax.plot_surface(X, Y, z, alpha=0.5)
Traceback (most recent call last):
  File "<pyshell#61>", line 1, in <module>
    ax.plot_surface(X, Y, z, alpha=0.5)
  File "/Users/NaN/anaconda/envs/py36/lib/python3.6/site-packages/mpl_toolkits/mplot3d/axes3d.py", line 1586, in plot_surface
    X, Y, Z = np.broadcast_arrays(X, Y, Z)
  File "/Users/NaN/anaconda/envs/py36/lib/python3.6/site-packages/numpy/lib/stride_tricks.py", line 259, in broadcast_arrays
    shape = _broadcast_shape(*args)
  File "/Users/NaN/anaconda/envs/py36/lib/python3.6/site-packages/numpy/lib/stride_tricks.py", line 193, in _broadcast_shape
    b = np.broadcast(*args[:32])
ValueError: shape mismatch: objects cannot be broadcast to a single shape

Answer 1

The arrays x, y, z need to be parametrized in two dimensions. One way of doing this is to use spherical coordinates as e.g. in Plot surfaces on a cube.

The remaining task is to distill the unique coordinates from the input data. I'm assuming here that there are only 2 distinct values per dimension.

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

def get_cube():   
    phi = np.arange(1,10,2)*np.pi/4
    Phi, Theta = np.meshgrid(phi, phi)

    x = np.cos(Phi)*np.sin(Theta)
    y = np.sin(Phi)*np.sin(Theta)
    z = np.cos(Theta)/np.sqrt(2)
    return x,y,z


points = np.array([[ 322697.1875    , 3663966.5       ,  -30000. ],
                   [ 325054.34375   , 3663966.5       ,  -30000. ],
                   [ 325054.34375   , 3665679.5       ,  -30000. ],
                   [ 322697.1875    , 3665679.5       ,  -30000. ],
                   [ 322697.1875    , 3663966.5       ,  -27703.12],
                   [ 325054.34375   , 3663966.5       ,  -27703.12],
                   [ 325054.34375   , 3665679.5       ,  -27703.12],
                   [ 322697.1875    , 3665679.5       ,  -27703.12]])

ux = np.unique(points[:,0])
uy = np.unique(points[:,1])
uz = np.unique(points[:,2])

x,y,z = get_cube()
offset = lambda X, o: o[0] + (X+.5)*np.diff(o)[0]


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

ax.plot_surface(offset(x, ux), offset(y, uy), offset(z, uz))

plt.show()